Foundationপ্রথম নীতি থেকে
LEVEL 11লেসন ১০/১৭কঠিন১ ঘণ্টা ৫ মিনিট

SIMD Programming — হাতে-লেখা AVX2 আর NEON কোড

SIMD Programming — Real AVX2 and NEON Intrinsics

SIMD-র ইতিহাস আর ধারণা আগেই দেখা হয়েছে — এবার হাতে সত্যিকারের AVX2/NEON intrinsics লিখে, compiler-এর auto-vectorization-এর সাথে তুলনা করে, আর AoS/SoA-এর মতো বাস্তব data-layout সিদ্ধান্ত নিয়ে দেখা হবে।

এই লেসন শেষে আপনি পারবেন

  • একটা সম্পূর্ণ, compile-যোগ্য AVX2 intrinsics প্রোগ্রাম লিখে scalar array-add loop হাতে vectorize করতে পারবেন
  • হাতে-লেখা AVX2 assembly, scalar -O0, আর compiler-এর auto-vectorized -O3 output পাশাপাশি তুলনা করে পার্থক্য identify করতে পারবেন
  • একই operation-এর ARM NEON সংস্করণ লিখতে পারবেন এবং একই pattern কীভাবে ভিন্ন ISA-তে transfer হয় তা ব্যাখ্যা করতে পারবেন
  • একটা নির্দিষ্ট সমস্যায় Array-of-Structures (AoS) বনাম Structure-of-Arrays (SoA) বেছে নেওয়ার সিদ্ধান্ত SIMD applicability-র ভিত্তিতে নিতে পারবেন
  • Alignment-এর ইতিহাস ও বর্তমান বাস্তবতা, আর horizontal reduction কেন genuinely awkward তা কোড দিয়ে দেখাতে পারবেন

আগে যা বোঝা থাকা দরকার

আগে এটা বুঝি

computer-architecture/simd লেসনে SIMD-র পুরো গল্প বলা হয়েছে — MMX থেকে AVX-512 পর্যন্ত ইতিহাস, packed register-এর ধারণা, আর কোন কোড data-parallel তা চেনার পদ্ধতি। এই লেসন সেই গল্প আর repeat করবে না — সরাসরি কোড লেখায় যাবে।

assembly/reading-compiler-output.mdx লেসনে আমরা দেখেছিলাম -O0 বনাম -O2-এ একই function কতটা ভিন্ন assembly-তে compile হয়। সেই লেসনের শেষে একটা ইঙ্গিত ছিল — -O3-এ auto-vectorization যোগ হয়, কিন্তু সেই বিস্তারিত আলোচনা “SIMD লেসনে” পিছিয়ে রাখা হয়েছিল। আজ সেই প্রতিশ্রুতি পূরণ হবে — scalar, হাতে-লেখা AVX2, আর compiler-এর auto-vectorized output, তিনটাকে পাশাপাশি রেখে দেখব।

মূল ধারণা

একই কাজ, তিনটা পথ

আমাদের কাজ চিরপরিচিত — দুটো float array যোগ করে তৃতীয়টায় রাখা।

পথ ১ — Scalar (baseline):

void add_scalar(float *c, const float *a, const float *b, int n) {
    for (int i = 0; i \< n; i++) {
        c[i] = a[i] + b[i];
    }
}

পথ ২ — হাতে-লেখা AVX2 intrinsics:

#include \<immintrin.h>

void add_avx2(float *c, const float *a, const float *b, int n) {
    int i = 0;
    for (; i + 8 \<= n; i += 8) {
        __m256 va = _mm256_loadu_ps(&a[i]);   // ৮টা float memory থেকে register-এ
        __m256 vb = _mm256_loadu_ps(&b[i]);   // আরেকটা ৮টা float
        __m256 vc = _mm256_add_ps(va, vb);    // এক instruction-এ ৮টা addition
        _mm256_storeu_ps(&c[i], vc);          // ফলাফল ৮টা একসাথে ফেরত
    }
    for (; i \< n; i++) {                      // n যদি 8-এর গুণিতক না হয়, বাকিটা scalar-এ
        c[i] = a[i] + b[i];
    }
}

প্রতিটা intrinsic-এর signature সুনির্দিষ্ট:

IntrinsicSignatureকাজ
_mm256_loadu_ps__m256 (float const *)৮টা float, memory থেকে register-এ, unaligned নিরাপদ
_mm256_add_ps__m256 (__m256, __m256)দুটো register-এর ৮টা lane-ই একসাথে যোগ
_mm256_storeu_psvoid (float *, __m256)register থেকে memory-তে, unaligned নিরাপদ

__m256 একটা compiler-built-in type — একটা ২৫৬-bit YMM register-এর প্রতিনিধিত্ব করে (৮টা 32-bit float হিসেবে packed, computer-architecture/simd লেসনের সেই packed-register figure মনে করুন)। এটা কোনো সাধারণ C struct/array না — compiler এটাকে সরাসরি একটা hardware register-এ রাখার চেষ্টা করে, ঠিক int/float-এর মতোই একটা primitive type হিসেবে ব্যবহার করা হয়।

i + 8 \<= n শর্তটা লক্ষ্য করুন — main loop শুধু তখনই চলবে যখন কমপক্ষে ৮টা element বাকি আছে। শেষের tail (n % 8 টা বাকি element, যদি থাকে) আলাদা scalar loop-এ সামলানো হয় — এটা computer-architecture/simd লেসনের auto-vectorization আলোচনায় দেখা “alignment ও tail handling” সমস্যার হাতে-লেখা সমাধান।

ভেতরে কী ঘটছে

Assembly-তে আসল পার্থক্য — তিনটা output পাশাপাশি

reading-compiler-output.mdx-এর পদ্ধতি অনুসরণ করে তিনটা ভিন্ন compile করি:

gcc -O0 -S add_scalar.c -o scalar_o0.s          # baseline, কোনো optimization না
gcc -O3 -march=native -S add_scalar.c -o auto_o3.s   # compiler নিজে চেষ্টা করে
gcc -O2 -mavx2 -S add_avx2.c -o hand_avx2.s     # আমাদের হাতে-লেখা, শুধু register allocation optimize

Scalar -O0 (সংক্ষেপে, মূল loop body):

movl    -4(%rbp), %eax          ; i লোড
cltq
movss   (%rsi,%rax,4), %xmm0    ; a[i]
movl    -4(%rbp), %eax
cltq
addss   (%rdx,%rax,4), %xmm0    ; + b[i]  (register-memory add)
movl    -4(%rbp), %eax
cltq
movss   %xmm0, (%rdi,%rax,4)    ; c[i] = ফলাফল

একটা element প্রতি iteration, reading-compiler-output.mdx-এর সেই পরিচিত “প্রতিটা sub-step আলাদা instruction” প্যাটার্ন।

Auto-vectorized -O3 (grep vaddps-এ যা পাওয়া যাবে):

vmovups (%rsi,%rax,4), %ymm0     ; a-এর ৮টা float, একসাথে
vaddps  (%rdx,%rax,4), %ymm0, %ymm0  ; +b, memory operand সরাসরি ব্যবহার
vmovups %ymm0, (%rdi,%rax,4)     ; c-তে ৮টা একসাথে লেখা
addq    $32, %rax                ; পরের ৮টা float (32 byte) এগিয়ে যাওয়া

হাতে-লেখা AVX2 (add_avx2.c থেকে) প্রায় হুবহু একই instruction sequence তৈরি করে — এটাই প্রত্যাশিত, কারণ intrinsics নিজেই সরাসরি এই instruction-গুলোর একটা thin wrapper। পার্থক্য mostly সূক্ষ্ম: compiler auto-vectorized সংস্করণে vaddps-এর memory operand সরাসরি ব্যবহার করেছে (load + add এক instruction-এ ফিউজ করে, reading-compiler-output.mdx-এর সেই instruction fusion ধারণা এখানে vector instruction-এও প্রযোজ্য), যেখানে হাতে-লেখা কোডে _mm256_loadu_ps আলাদাভাবে লেখা থাকায় compiler-নির্ভর করে সেটাকে ফিউজ করবে কি না।

Alignment — ইতিহাস আর বর্তমান বাস্তবতা

_mm256_loadu_ps-এর u মানে unaligned — যেকোনো byte ঠিকানা থেকে load করতে পারে। এর বিপরীত _mm256_load_ps (কোনো u ছাড়া) aligned load — এটা দ্রুত হতে পারে, কিন্তু একটা কঠোর শর্তে: ঠিকানা অবশ্যই ৩২-byte boundary-তে aligned হতে হবে (ঠিকানা mod ৩২ = ০), নাহলে প্রোগ্রাম crash করে (#GP fault)।

float *a = aligned_alloc(32, n * sizeof(float));   // ৩২-byte aligned allocation
__m256 va = _mm256_load_ps(&a[0]);                 // নিরাপদ, aligned
__m256 vb = _mm256_load_ps(&a[3]);                 // CRASH! a[3] সাধারণত 32-এর গুণিতকে aligned না

AoS বনাম SoA — SIMD আদৌ প্রযোজ্য কি না তার পূর্বশর্ত

SIMD-র সবচেয়ে বাস্তব, ইন্টারভিউ-প্রাসঙ্গিক প্রশ্নটা হলো ডেটা layout-এর — কারণ _mm256_loadu_ps চায় ৮টা contiguous float, একটা struct-এর ভেতরের এলোমেলো ছড়ানো field না।

ধরুন একটা particle simulation — প্রতিটা particle-এর position আছে:

// Array of Structures (AoS) -- স্বাভাবিক, "object-oriented" instinct
typedef struct {
    float x, y, z;
    float vx, vy, vz;
} Particle;

Particle particles[N];

Memory-তে particles এভাবে সাজানো: x0,y0,z0,vx0,vy0,vz0, x1,y1,z1,vx1,... — প্রতিটা particle-এর সব field পাশাপাশি, কিন্তু একই field-এর বিভিন্ন particle একে অপরের থেকে ৬টা float দূরত্বে। যদি আমরা সব particle-এর x-কেই একসাথে vectorize করতে চাই (x[i] += vx[i] * dt), _mm256_loadu_ps(&particles[i].x) কাজ করবে না ঠিকভাবে — এটা x0, y0, z0, vx0, vy0, vz0, x1, y1-কে একসাথে একটা register-এ তুলে আনবে, যা fizically ৮টা float হলেও logically আটটা ভিন্ন জিনিস, একই field-এর ৮টা instance না। SIMD lane-ভিত্তিক অপারেশন (সব lane-এ একই কাজ) তাই এখানে অর্থহীন হয়ে যায়।

// Structure of Arrays (SoA) -- SIMD-friendly পুনর্বিন্যাস
typedef struct {
    float x[N], y[N], z[N];
    float vx[N], vy[N], vz[N];
} ParticleSoA;

ParticleSoA p;

এখন p.x[i], p.x[i+1], ..., p.x[i+7] সত্যিই contiguous, একই field-এর ৮টা instance — সরাসরি _mm256_loadu_ps(&p.x[i])-এ vectorize-যোগ্য, প্রতিটা lane একই ধরনের ডেটা (একটা x-coordinate) নিয়ে কাজ করছে।

AoS (memory-তে):
[x0 y0 z0 vx0 vy0 vz0][x1 y1 z1 vx1 vy1 vz1][x2 y2 z2 ...]
 ↑ একটা 8-wide load এখানে ধরলে মিশ্রিত, ভিন্ন-অর্থের ডেটা একসাথে আসবে

SoA (memory-তে):
[x0 x1 x2 x3 x4 x5 x6 x7 ...][y0 y1 y2 y3 ...][z0 z1 z2 ...]
 ↑ একটা 8-wide load এখানে ঠিক ৮টা x-value — সরাসরি vectorize-যোগ্য
AoS বনাম SoA — memory-তে বিন্যাস, আর SIMD-এর জন্য কোনটা কাজ করে।

Horizontal reduction — SIMD-র একটা সৎ সীমাবদ্ধতা

_mm256_add_ps স্বাভাবিকভাবেই lane-ভিত্তিক (element-wise) কাজ করে — প্রতিটা lane independent, কোনো lane আরেকটার সাথে মিশে না। কিন্তু কখনো কখনো আমাদের lane-এর মধ্যেই একটা মিশ্রণ দরকার হয় — যেমন একটা __m256-এর ভেতরের ৮টা float-এর যোগফল বের করা (reduction)। এটা SIMD-র natural per-lane চিন্তাভাবনার বিপরীত, তাই instruction-সংখ্যায় সত্যিই বেশি ব্যয়বহুল:

#include \<immintrin.h>

float hsum256_ps(__m256 v) {
    __m128 lo   = _mm256_castps256_ps128(v);      // নিচের ৪টা lane
    __m128 hi   = _mm256_extractf128_ps(v, 1);    // উপরের ৪টা lane
    lo = _mm_add_ps(lo, hi);                      // ৪টা partial sum
    __m128 shuf = _mm_movehdup_ps(lo);            // [1,1,3,3] pattern-এ broadcast
    __m128 sums = _mm_add_ps(lo, shuf);           // আরও গুটিয়ে আনা
    shuf = _mm_movehl_ps(shuf, sums);
    sums = _mm_add_ss(sums, shuf);                // শেষ স্কেলার যোগ
    return _mm_cvtss_f32(sums);                   // register থেকে সাধারণ float
}

৮টা independent addition করতে ছিল ১টা instruction (_mm256_add_ps)। ৮টা lane-কে একটা মাত্র সংখ্যায় গুটিয়ে আনতে লাগল ৬টা instruction — shuffle, extract, cast, একাধিক ধাপে add। এটাই “horizontal operation SIMD-তে awkward” বক্তব্যের সুনির্দিষ্ট প্রমাণ — SIMD hardware lane-গুলোকে independent রাখার জন্য optimize করা, lane-গুলোর মধ্যে ডেটা “পাশাপাশি সরানো” (shuffle) তার স্বাভাবিক শক্তির জায়গা না।

উদাহরণ

ARM NEON — একই প্যাটার্ন, ভিন্ন ISA

assembly/arm64-basics.mdx লেসনে আমরা ARM64-এর x0-x30 general-purpose register set দেখেছিলাম। NEON এর জন্য একটা সম্পূর্ণ আলাদা register file আছে — v0-v31, প্রতিটা ১২৮-bit (AVX２-এর অর্ধেক প্রস্থ, তাই এক লাফে ৪টা float, ৮টা না)।

#include \<arm_neon.h>

void add_neon(float *c, const float *a, const float *b, int n) {
    int i = 0;
    for (; i + 4 \<= n; i += 4) {
        float32x4_t va = vld1q_f32(&a[i]);   // ৪টা float memory থেকে
        float32x4_t vb = vld1q_f32(&b[i]);
        float32x4_t vc = vaddq_f32(va, vb);  // এক instruction-এ ৪টা addition
        vst1q_f32(&c[i], vc);                // ৪টা একসাথে ফেরত
    }
    for (; i \< n; i++) {
        c[i] = a[i] + b[i];
    }
}
AVX2 (x86-64)NEON (ARM64/AArch64)
Register type__m256float32x4_t
Register width২৫৬ bit১২৮ bit
float32 lane সংখ্যা৮৪
Load_mm256_loadu_psvld1q_f32
Add_mm256_add_psvaddq_f32
Store_mm256_storeu_psvst1q_f32
Header\<immintrin.h>\<arm_neon.h>
AVX2 বনাম NEON — একই shape, ভিন্ন width আর নাম।

প্যাটার্নটা হুবহু একই — load, op, store, tail handling — শুধু নাম আর width আলাদা। যে প্রোগ্রামার একবার এই “load-compute-store” কাঠামো বুঝেছেন, তিনি যেকোনো SIMD ISA-তে (AVX-512, RISC-V-এর Vector Extension) প্রায় একই ধরনের কোড লিখতে পারবেন — এটাই এই লেসনের একটা গুরুত্বপূর্ণ বার্তা: SIMD intrinsics শেখা মানে একটা নির্দিষ্ট instruction set মুখস্থ করা না, বরং এই shape-টা চেনা।

নিজে চালিয়ে দেখুন

EXPERIMENT

Scalar বনাম Auto-vectorized বনাম হাতে-লেখা AVX2 — তিনটাই মাপুন

Linux/macOS, gcc বা clang, AVX2-সক্ষম CPU· ২৫ মিনিট
// simd_compare.c
#include \<immintrin.h>
#include \<stdio.h>
#include \<stdlib.h>
#include \<time.h>

#define N (1 \<\< 26)   // ~৬৭ মিলিয়ন float

// এখানে restrict ইচ্ছাকৃতভাবে বাদ -- compiler-কে aliasing নিয়ে
// অনিশ্চিত রাখতে, auto-vectorization ব্যর্থ হওয়ার শর্ত তৈরি করতে
void add_scalar(float *c, float *a, float *b, int n) {
    for (int i = 0; i \< n; i++) c[i] = a[i] + b[i];
}

void add_avx2(float *c, float *a, float *b, int n) {
    int i = 0;
    for (; i + 8 \<= n; i += 8) {
        __m256 va = _mm256_loadu_ps(&a[i]);
        __m256 vb = _mm256_loadu_ps(&b[i]);
        _mm256_storeu_ps(&c[i], _mm256_add_ps(va, vb));
    }
    for (; i \< n; i++) c[i] = a[i] + b[i];
}

double bench(void (*fn)(float*, float*, float*, int),
             float *c, float *a, float *b, int n, int repeats) {
    struct timespec t0, t1;
    clock_gettime(CLOCK_MONOTONIC, &t0);
    for (int r = 0; r \< repeats; r++) fn(c, a, b, n);
    clock_gettime(CLOCK_MONOTONIC, &t1);
    return (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) / 1e9;
}

int main(void) {
    float *a = malloc(N * sizeof(float));
    float *b = malloc(N * sizeof(float));
    float *c = malloc(N * sizeof(float));
    for (int i = 0; i \< N; i++) { a[i] = 1.5f; b[i] = 2.5f; }

    double t_scalar = bench(add_scalar, c, a, b, N, 10);
    double t_avx2   = bench(add_avx2,   c, a, b, N, 10);

    printf("scalar (no restrict): %.4fs\n", t_scalar);
    printf("hand-written AVX2:    %.4fs  (%.2fx)\n",
           t_avx2, t_scalar / t_avx2);
    return 0;
}
# ১. auto-vectorization কী পেল তা প্রথমে যাচাই করুন
gcc -O3 -march=native -fopt-info-vec -c simd_compare.c -o /dev/null

# ২. পুরো প্রোগ্রাম চালান
gcc -O3 -march=native -mavx2 -o simd_compare simd_compare.c && ./simd_compare

সাধারণ ফলাফল — restrict ইচ্ছাকৃতভাবে বাদ দেওয়ায় compiler-এর জন্য add_scalar-কে নিরাপদে vectorize করা কঠিন (aliasing অনিশ্চয়তা, computer-architecture/simd লেসনের সেই আলোচনা):

simd_compare.c:12:5: missed: couldn't vectorize loop
simd_compare.c:12:5: missed: not vectorized: unsupported data-ref

scalar (no restrict): 0.298s
hand-written AVX2:    0.041s  (7.27x)

হাতে-লেখা AVX2 এখানে স্পষ্টভাবে জেতে — কারণ intrinsics ব্যবহার করলে aliasing-এর প্রশ্নই আসে না (প্রোগ্রামার নিজেই নিশ্চিত করছেন, compiler-কে অনুমতি নিতে হয় না)। restrict keyword যোগ করে add_scalar-কে আবার কম্পাইল করুন — তখন -fopt-info-vec সফল vectorization দেখাবে, আর scalar/AVX2-এর পার্থক্য অনেকটা কমে আসবে (কারণ এবার compiler নিজেই প্রায় একই AVX２ কোড generate করতে পারছে)। এই দুইটা রান একসাথে দেখেই এই experiment-এর মূল কথা প্রমাণ হয়: compiler সফল হলে হাতে লেখার প্রয়োজন কমে যায়, কিন্তু ব্যর্থ হলে (এখানে যেমন restrict ছাড়া) হাতে-লেখা কোডই একমাত্র নির্ভরযোগ্য বিকল্প।

এটা কী প্রমাণ করে

হাতে-লেখা intrinsics আর compiler-এর ভালো auto-vectorization প্রায় একই গতি দেয় যখন উভয়ই সফল হয় — কিন্তু compiler ব্যর্থ হলে (aliasing অনিশ্চয়তা) হাতে-লেখা কোডই একমাত্র নির্ভরযোগ্য পথ।

নিজে বানান

BUILD IT

AoS particle update-কে SoA-তে রূপান্তর করে vectorize করুন

C (AVX2 intrinsics) · ●●●●○
  1. একটা AoS Particle struct (x, y, z, vx, vy, vz) বানিয়ে N=1,000,000 particle-এর একটা array বানান, random initial value দিয়ে
  2. AoS সংস্করণে scalar position-update লিখুন: প্রতিটা particle-এর x/y/z += vx/vy/vz * dt, সময় মাপুন
  3. একই ডেটা SoA (আলাদা x[], y[], z[], vx[], vy[], vz[] array) হিসেবে পুনর্গঠন করুন
  4. SoA ডেটায় AVX2 দিয়ে vectorized position-update লিখুন (dt কে _mm256_set1_ps দিয়ে broadcast করে), সময় মাপুন
  5. দুই সংস্করণের ফলাফল element-by-element তুলনা করে সংখ্যাগতভাবে অভিন্ন প্রমাণ করুন, তারপর speedup রিপোর্ট করুন

মূল ধারণা: একই গাণিতিক কাজ (position += velocity * dt), কিন্তু layout বদলানোর ফলে SIMD প্রযোজ্য হয় কি না তা সম্পূর্ণ বদলে যায় — এই লেসনের AoS/SoA আলোচনার সবচেয়ে concrete প্রয়োগ।

#include \<immintrin.h>
#include \<stdio.h>
#include \<stdlib.h>
#include \<time.h>
#include \<math.h>

#define N (1 \<\< 20)   // ~10 লক্ষ particle

typedef struct { float x, y, z, vx, vy, vz; } Particle;

typedef struct {
    float *x, *y, *z, *vx, *vy, *vz;
} ParticleSoA;

void update_aos(Particle *p, int n, float dt) {
    for (int i = 0; i \< n; i++) {
        p[i].x += p[i].vx * dt;
        p[i].y += p[i].vy * dt;
        p[i].z += p[i].vz * dt;
    }
}

void update_soa_avx2(ParticleSoA *p, int n, float dt) {
    __m256 vdt = _mm256_set1_ps(dt);        // dt-কে ৮ বার broadcast
    int i = 0;
    for (; i + 8 \<= n; i += 8) {
        __m256 x  = _mm256_loadu_ps(&p->x[i]);
        __m256 vx = _mm256_loadu_ps(&p->vx[i]);
        _mm256_storeu_ps(&p->x[i], _mm256_fmadd_ps(vx, vdt, x));  // x + vx*dt, এক instruction-এ

        __m256 y  = _mm256_loadu_ps(&p->y[i]);
        __m256 vy = _mm256_loadu_ps(&p->vy[i]);
        _mm256_storeu_ps(&p->y[i], _mm256_fmadd_ps(vy, vdt, y));

        __m256 z  = _mm256_loadu_ps(&p->z[i]);
        __m256 vz = _mm256_loadu_ps(&p->vz[i]);
        _mm256_storeu_ps(&p->z[i], _mm256_fmadd_ps(vz, vdt, z));
    }
    for (; i \< n; i++) {
        p->x[i] += p->vx[i] * dt;
        p->y[i] += p->vy[i] * dt;
        p->z[i] += p->vz[i] * dt;
    }
}

int main(void) {
    float dt = 0.01f;

    Particle *aos = malloc(N * sizeof(Particle));
    ParticleSoA soa = {
        .x = malloc(N*sizeof(float)),  .y = malloc(N*sizeof(float)),
        .z = malloc(N*sizeof(float)),  .vx = malloc(N*sizeof(float)),
        .vy = malloc(N*sizeof(float)), .vz = malloc(N*sizeof(float)),
    };

    srand(1);
    for (int i = 0; i \< N; i++) {
        float x = (float)(rand() % 1000), y = (float)(rand() % 1000);
        float z = (float)(rand() % 1000);
        float vx = (float)(rand() % 10 - 5), vy = (float)(rand() % 10 - 5);
        float vz = (float)(rand() % 10 - 5);
        aos[i] = (Particle){x, y, z, vx, vy, vz};
        soa.x[i]=x; soa.y[i]=y; soa.z[i]=z;
        soa.vx[i]=vx; soa.vy[i]=vy; soa.vz[i]=vz;
    }

    struct timespec t0, t1;

    clock_gettime(CLOCK_MONOTONIC, &t0);
    for (int r = 0; r \< 100; r++) update_aos(aos, N, dt);
    clock_gettime(CLOCK_MONOTONIC, &t1);
    double t_aos = (t1.tv_sec-t0.tv_sec)+(t1.tv_nsec-t0.tv_nsec)/1e9;

    clock_gettime(CLOCK_MONOTONIC, &t0);
    for (int r = 0; r \< 100; r++) update_soa_avx2(&soa, N, dt);
    clock_gettime(CLOCK_MONOTONIC, &t1);
    double t_soa = (t1.tv_sec-t0.tv_sec)+(t1.tv_nsec-t0.tv_nsec)/1e9;

    printf("AoS scalar:      %.4fs\n", t_aos);
    printf("SoA AVX2:        %.4fs  (%.2fx)\n", t_soa, t_aos / t_soa);

    // Correctness -- ১০০ ধাপ পর x-এর মান তুলনা (floating-point tolerance সহ)
    int mismatches = 0;
    for (int i = 0; i \< N; i++) {
        if (fabsf(aos[i].x - soa.x[i]) > 1e-2f) mismatches++;
    }
    printf("mismatches: %d / %d\n", mismatches, N);
    return 0;
}
gcc -O2 -mavx2 -mfma -o particle_update particle_update.c -lm && ./particle_update

প্রত্যাশিত আউটপুট আকৃতি:

AoS scalar:      0.612s
SoA AVX2:        0.089s  (6.88x)
mismatches: 0 / 1048576

নিজে বাড়ান:

  1. AoS ডেটাকে AVX２ দিয়ে সরাসরি vectorize করার চেষ্টা করুন (কোনো SoA রূপান্তর ছাড়া) — কেন এটা স্বাভাবিকভাবে সম্ভব না, বা কেন একটা _mm256_loadu_ps কল ভুল ডেটা মিশিয়ে ফেলবে তা কোড লিখে directly অনুভব করুন
  2. _mm256_fmadd_ps-এর বদলে আলাদা _mm256_mul_ps + _mm256_add_ps দিয়ে লিখে তুলনা করুন — FMA কতটা পার্থক্য আনে
  3. float এর বদলে double দিয়ে পুরো experiment আবার করুন (_mm256_loadu_pd, ৪টা lane) — কেন speedup কম হয় তা lane-সংখ্যার হিসাব দিয়ে ব্যাখ্যা করুন

বাস্তব সিস্টেমে

যেখানে এই কোডটাই সরাসরি চলে

  • simdjson — একটা জনপ্রিয় open-source JSON parser যেটা পুরোপুরি hand-written AVX2/NEON intrinsics দিয়ে তৈরি, traditional parser-এর চেয়ে ৪-২৫ গুণ দ্রুত। এর কোডবেসেই ঠিক এই লেসনের load/compute/store প্যাটার্নের বাস্তব, production-grade উদাহরণ পাওয়া যায়।

  • ECS (Entity Component System) গেম architecture — আধুনিক গেম ইঞ্জিন (Unity DOTS, বহু custom engine) স্বেচ্ছায় SoA layout বেছে নেয় ঠিক এই লেসনের কারণেই — হাজার হাজার entity-র position আপডেট করা এই লেসনের build অংশের particle simulation-এরই বড় স্কেল সংস্করণ।

  • oneDNN / ML inference kernel — Intel-এর deep-learning kernel library-র সবচেয়ে ভেতরের loop hand-tuned AVX２/AVX-512 intrinsics, matrix multiply আর convolution-এর জন্য — এখানে horizontal reduction এড়ানোর কৌশল (accumulator-কে যতটা সম্ভব lane-ভিত্তিক রাখা, শেষেই একবার reduce করা) একটা মূল ডিজাইন নীতি।

  • NumPy-র C loop। Python-এর NumPy যখন a + b (দুটো array) চালায়, ভেতরে একটা C loop চলে যা compiler auto-vectorization-এর উপর নির্ভর করে (আর কিছু জায়গায় explicit intrinsics ব্যবহার করে) — এই কারণেই NumPy pure Python loop-এর চেয়ে বহুগুণ দ্রুত, ঠিক এই লেসনের experiment-এর মতোই একটা কারণে।

  • Database columnar engine (ClickHouse, DuckDB) — column-store database প্রাকৃতিকভাবেই SoA layout (প্রতিটা column একটা আলাদা, contiguous array) — যা তাদের “vectorized execution” কৌশলকে সম্ভব করে, আগের লেসনের prefetching আলোচনার সাথেও সরাসরি সংযুক্ত (contiguous column scan prefetcher-friendlyও বটে)।

  • Rust-এর std::simd ও C++23-এর std::simd — portable SIMD abstraction যা একই কোড লিখে AVX２ আর NEON দুটোতেই compile করতে দেয় — এই লেসনে হাতে যা #ifdef-এর মাধ্যমে করতে হতো (দুটো আলাদা ফাংশন), সেটাকে ভাষা-স্তরে সরল করার একটা প্রচেষ্টা।

  • libjpeg-turbo, x264/FFmpeg — computer-architecture/simd লেসনে এদের উল্লেখ ছিল conceptual স্তরে; এই লেসনের পর তাদের actual source code-এ গেলে ঠিক এই load-op-store প্যাটার্নই বারবার দেখা যাবে, শুধু বহুগুণ বেশি জটিল operation-এ (IDCT, motion estimation)।

যে ভুলগুলো সবাই করে

“Unaligned load এখনও সবসময় উল্লেখযোগ্যভাবে ধীর, তাই সবসময় aligned memory ব্যবহার করা উচিত।”

Nehalem (২০০৮)-এর পরবর্তী hardware-এ _mm256_loadu_ps আর _mm256_load_ps-এর মধ্যে পার্থক্য প্রায় নগণ্য যদি প্রকৃত ঠিকানা এমনিতেই aligned থাকে — hardware নিজে থেকেই দ্রুত পথ শনাক্ত করে। loadu ব্যবহার করাই বরং নিরাপদ ডিফল্ট, কারণ এটা কখনো fault করে না, আর প্রায় কোনো performance দাম ছাড়াই কাজ করে। load (কঠোর alignment দাবি করা) শুধু তখনই বেছে নিন যখন alignment guarantee করা সহজ এবং measurement সত্যিই একটা পার্থক্য দেখায় — অনুমান করে না।

“যেকোনো array of struct-কে সরাসরি SIMD দিয়ে vectorize করা যায়, শুধু intrinsics লিখলেই হয়।”

এই লেসনের AoS/SoA আলোচনা ঠিক এর বিপরীত দেখিয়েছে। SIMD load instruction contiguous memory থেকে একই ধরনের ডেটা আশা করে — AoS layout-এ একটা 8-wide load struct-এর বিভিন্ন field মিশিয়ে ফেলবে, যা logically অর্থহীন। Intrinsics লেখার আগেই ডেটা-লেআউট পরীক্ষা করা আবশ্যক — প্রয়োজনে SoA-তে রূপান্তর করেই vectorization সম্ভব হয় (অথবা AVX2/AVX-512-এর ধীরগতির gather/scatter instruction ব্যবহার করতে হয়, যা computer-architecture/simd লেসনে উল্লেখ করা হয়েছিল)।

“Horizontal reduction (যেমন sum) SIMD-তে সবসময় সমান কষ্টকর, ISA নির্বিশেষে।”

এই লেসনের NEON উদাহরণই এই দাবি ভুল প্রমাণ করে। AVX２-তে horizontal sum-এর জন্য ৬টা instruction লাগে (shuffle-ভিত্তিক কৌশল), কিন্তু AArch64 NEON-এ vaddvq_f32 নামের একটা dedicated instruction সরাসরি এই কাজ করে, এক call-এই। এটা প্রমাণ করে horizontal-reduction-এর “কষ্ট” কোনো universal SIMD-সীমা না, বরং একটা নির্দিষ্ট ISA-র hardware design সিদ্ধান্ত — তাই কর্মক্ষমতা-সংবেদনশীল কোডে সবসময় target ISA-র নির্দিষ্ট capability (intrinsics guide) দেখেই সিদ্ধান্ত নেওয়া উচিত, একটা ISA-র অভিজ্ঞতা থেকে সরাসরি সাধারণীকরণ করা উচিত না।

বুঝেছেন কি না দেখুন

1

একটা ফাংশনে restrict keyword ব্যবহার না করলে GCC -O3-এও প্রায়ই vectorize করতে ব্যর্থ হয়, যদিও হাতে-লেখা AVX2 কোড ঠিকই কাজ করে। এই পার্থক্যের কারণ সুনির্দিষ্টভাবে ব্যাখ্যা করুন।

প্রয়োগ

পার্থক্যটা কে দায়িত্ব নিচ্ছে correctness নিশ্চিত করার তার মধ্যে। restrict ছাড়া, compiler-কে ধরে নিতে হয় a, b, c pointer তিনটা memory-তে overlap করতেও পারে (যেমন c আসলে a-এরই একটা অংশ) — যদি সত্যিই overlap করে, ৮টা element একসাথে প্রসেস করা (vectorize করা) memory access-এর ক্রম বদলে দেবে, যা ভুল ফলাফল দিতে পারে। যেহেতু compiler নিশ্চিতভাবে জানে না overlap আছে কি না, এটা নিরাপদ দিকে ভুল করে — vectorize করে না।

হাতে-লেখা AVX２ কোডে, প্রোগ্রামার নিজেই (implicit ভাবে, ফাংশন কল করার সময়) নিশ্চিত করছেন যে a, b, c overlap করে না — কম্পাইলারকে কোনো অনুমান করতে হচ্ছে না, কারণ _mm256_loadu_ps/_mm256_storeu_ps ঠিক যা লেখা আছে তাই করবে, কোনো “কি যদি” বিশ্লেষণ ছাড়াই। এই “guarantee-এর দায়িত্ব compiler থেকে প্রোগ্রামারে সরে যাওয়া”-ই মূল কারণ — restrict যোগ করলে ঠিক এই একই guarantee compiler-কেও দেওয়া হয়, আর তখন auto-vectorization সাধারণত সফল হয়।

2

একটা টিম প্রতিটা struct Pixel { uint8_t r, g, b, a; }-এর r channel-এ brightness যোগ করতে চায়, লক্ষ লক্ষ pixel-এ। AoS (Pixel pixels[N]) layout-এই থেকে গিয়ে তারা _mm256_loadu_si256 দিয়ে সরাসরি ৩২ byte একসাথে load করে brightness যোগ করার চেষ্টা করল। এতে কী ভুল হবে, আর কেন?

যুক্তি

এই লোডটা যান্ত্রিকভাবে “কাজ করবে” (কোনো crash হবে না, ৩২ byte-ই load হবে), কিন্তু ফলাফল ভুল হবে — কারণ সেই ৩২ byte-এ শুধু r channel নেই, বরং interleaved r,g,b,a,r,g,b,a,... প্যাটার্নে মিশ্রিত আটটা সম্পূর্ণ pixel (৪ byte × ৮ = ৩২ byte)। যদি এই পুরো ৩২ byte-এ একটা saturating-add instruction চালানো হয় (ভুলবশত), g, b, a channel-ও brightness দিয়ে বদলে যাবে — যা উদ্দেশ্য ছিল না, শুধু r বদলানোর কথা ছিল।

এটা এই লেসনের AoS/SoA আলোচনার হুবহু একই সমস্যা, শুধু ভিন্ন domain-এ (image processing বনাম particle simulation)। সঠিক সমাধান হয় দুইটার একটা:

  1. De-interleave করে SoA-তে রূপান্তর — r[N], g[N], b[N], a[N] চারটা আলাদা array, তারপর শুধু r[]-এ সরাসরি vectorize করা।
  2. Shuffle instruction দিয়ে channel আলাদা করা — _mm256_shuffle_epi8-এর মতো instruction দিয়ে interleaved ডেটা থেকে শুধু r channel-গুলো একটা register-এ জড়ো করা, তারপর brightness যোগ করে আবার শাফল করে ফেরত রাখা — কাজ করে, কিন্তু অতিরিক্ত shuffle instruction-এর দাম আছে (এই লেসনের horizontal-reduction আলোচনার মতোই একটা “lane-এর মধ্যে ডেটা সরানো ব্যয়বহুল” নীতির আরেকটা রূপ)।

computer-architecture/simd লেসনের build অংশেই এই ঠিক এই “RGB channel interleaving SIMD-কে জটিল করে তোলে” পয়েন্টটা ইঙ্গিত করা হয়েছিল — এই প্রশ্নটাই সেই ইঙ্গিতের সম্পূর্ণ ব্যাখ্যা।

3

AVX2-তে _mm256_add_ps একটা instruction, কিন্তু hsum256_ps (horizontal sum) ছয়টা instruction লাগে। এই asymmetry hardware ডিজাইনের কোন মৌলিক নীতি থেকে আসে?

যুক্তি

মূল নীতি — SIMD hardware lane independence-এর জন্য optimize করা, lane-এর মধ্যে ডেটা চলাচলের জন্য না।

_mm256_add_ps-এ প্রতিটা লেনের ALU সম্পূর্ণ independently কাজ করে — লেন ৩-এর addition লেন ৫-এর কোনো ডেটার উপর নির্ভর করে না। এই independence-ই hardware-কে সব লেন একই cycle-এ, সমান্তরালে চালাতে দেয় — এটাই মূল কারণ কেন SIMD দ্রুত।

কিন্তু horizontal sum-এর জন্য ঠিক তার বিপরীত দরকার — লেন ০-এর ডেটা লেন ১-এর সাথে মিশতে হবে, তারপর সেই ফলাফল লেন ২-এর সাথে, ইত্যাদি। এর জন্য একটা আলাদা ধরনের hardware circuit লাগে — shuffle/permute network, যা লেনগুলোর মধ্যে ডেটা routing করে। এই routing hardware প্রতিটা সম্ভাব্য pattern-এর জন্য আলাদা instruction হিসেবে exposed (movehdup, movehl, extractf128, ইত্যাদি) কারণ “সব সম্ভাব্য শাফল প্যাটার্নের জন্য একটা single universal instruction” ডিজাইন করা হার্ডওয়্যার-ব্যয়বহুল এবং কম প্রয়োজনীয় (বেশিরভাগ SIMD কোড lane-independent অপারেশনই বেশি ব্যবহার করে)।

তাই এই asymmetry দুর্ঘটনাজনিত না — এটা একটা সচেতন design trade-off: “যা বেশি ব্যবহার হয় (lane-independent ops) তা সস্তা করো, যা কম ব্যবহার হয় (cross-lane ops) তার দাম মেনে নাও” — ঠিক এই module-এর বারবার আসা area-বনাম-frequency trade-off-এরই আরেকটা রূপ।

4

একটা codebase-এ ৫০টা জায়গায় Vector3 { float x, y, z; } struct ব্যবহার হয় (physics engine-এর মতো), সবসময় AoS হিসেবে। টিম চিন্তা করছে পুরো codebase SoA-তে রূপান্তর করবে কি না, শুধু SIMD-এর সুবিধা নিতে। কোন প্রশ্নগুলো আগে জিজ্ঞেস করা উচিত?

ডিজাইন

এই সিদ্ধান্ত নেওয়ার আগে অন্তত তিনটা প্রশ্নের উত্তর দরকার — কারণ SoA কোনো বিনামূল্যের upgrade না, একটা trade-off।

১. Access pattern আসলে কী প্রাধান্য পায়? যদি কোড বেশিরভাগ সময় “একটা vector নিয়ে সম্পূর্ণ কাজ” করে (normalize(v), dot(a, b) — যেখানে একটা single vector-এর x,y,z একসাথে দরকার), AoS-ই বরং ভালো — cache locality একটা vector-এর মধ্যেই থাকে। SoA শুধু তখনই জেতে যখন কোড “হাজার হাজার vector-এর একটা field নিয়ে একসাথে কাজ” করে (batch update, এই লেসনের build অংশের মতো)। ৫০টা জায়গার প্রতিটাতে এই pattern আলাদা হতে পারে — একটা blanket সিদ্ধান্ত ভুল হতে পারে।

২. Rewrite-এর খরচ কতটা? ৫০টা জায়গায় ব্যবহৃত একটা core type বদলানো মানে সব ৫০টা জায়গার কোড, test, আর সম্ভবত API পরিবর্তন — এটা একটা উল্লেখযোগ্য engineering প্রজেক্ট, শুধু “AVX২ ব্যবহার করতে চাই” বলার মতো সহজ না।

৩. আসল bottleneck কোথায়, profile করা হয়েছে কি? যদি physics update আসলে CPU সময়ের ৫% হয় (বাকি ৯৫% অন্য কোথাও — rendering, I/O), SIMD দিয়ে physics ৫× দ্রুত করলেও সামগ্রিক প্রোগ্রাম মাত্র সামান্য দ্রুত হবে (Amdahl-এর যুক্তি, পরের ভাগে বিস্তারিত আসবে) — তখন rewrite-এর খরচ সামগ্রিক লাভের তুলনায় ন্যায্য না-ও হতে পারে।

ব্যবহারিক মধ্যপন্থা: সম্পূর্ণ SoA না করে, শুধু hot path-এ (profile করে চিহ্নিত করা batch-update loop) স্থানীয়ভাবে SoA রূপান্তর করা — বাকি codebase AoS-ই থাকতে পারে, শুধু performance-critical অংশে conversion করে (loop শুরুতে AoS→SoA, শেষে ফেরত) বা সেই অংশটুকু সম্পূর্ণ SoA দিয়ে redesign করে। এভাবে বেশিরভাগ লাভ পাওয়া যায়, পুরো codebase-এর ঝুঁকি ছাড়াই।

5

_mm256_load_ps(&arr[3])-এ (aligned load) একটা প্রোগ্রাম crash করছে, কিন্তু _mm256_loadu_ps(&arr[3])-এ ঠিকমতো চলছে। arr একটা float arr[100] array, সাধারণ malloc/স্ট্যাক allocation দিয়ে তৈরি। কী ঘটছে?

প্রয়োগ

এটাই ঠিক এই লেসনের alignment আলোচনার একটা সরাসরি প্রয়োগ। _mm256_load_ps দাবি করে ঠিকানা ৩২-byte boundary-তে aligned হতে হবে — মানে ঠিকানা mod ৩২ = ০। সাধারণ malloc (কোনো বিশেষ alignment specify না করলে) সাধারণত শুধু ৮ বা ১৬-byte alignment নিশ্চিত করে (platform-নির্ভর), ৩২-byte না।

&arr[3] মানে arr-এর ঠিকানা থেকে 3 × 4 = 12 byte এগিয়ে। এমনকি যদি arr নিজে ৩২-byte aligned হয় (যেমন aligned_alloc(32, ...) দিয়ে তৈরি হতো), arr + 12 bytes অবশ্যই arr-এর alignment-এর সাথে সামঞ্জস্যপূর্ণ থাকবে না (12 mod 32 ≠ 0)। তাই _mm256_load_ps(&arr[3]) প্রায় নিশ্চিতভাবেই misaligned ঠিকানায় পড়বে — hardware সেটা মেনে নেয় না, একটা #GP (general protection) fault তোলে, প্রোগ্রাম crash করে।

_mm256_loadu_ps কোনো alignment দাবি করে না — একই ঠিকানায় (aligned হোক বা না হোক) নির্বিঘ্নে কাজ করে। এই কারণেই এই লেসনে loadu-কেই ডিফল্ট, নিরাপদ পছন্দ হিসেবে সুপারিশ করা হয়েছে — শুধু তখনই load বেছে নিন যখন আপনি নিশ্চিতভাবে জানেন প্রতিটা access-point (এই ক্ষেত্রে &arr[3], &arr[3+8], ইত্যাদি প্রতিটা loop iteration-এর ঠিকানা) ৩২-byte aligned থাকবে — যা এখানে সত্যি না, কারণ 3 কোনো ৮-এর গুণিতক না (৮টা float = ৩২ byte, তাই index ০, ৮, ১৬, … হলে align থাকত)।

এরপর কী

এই লেসন পর্যন্ত আমরা “একটা core-এর ভেতরে” যতটা parallelism সম্ভব তার প্রায় প্রতিটা রূপ দেখেছি — pipeline-এর ILP, out-of-order execution, আর এখন SIMD-এর data parallelism (একটা instruction, কয়েক লেন)। পরের ভাগে আমরা এই একই lane-ভিত্তিক ধারণাটাকে সম্পূর্ণ ভিন্ন এক স্কেলে নিয়ে যাব — GPU-তে, যেখানে ৮ বা ১৬টা লেন না, বরং একসাথে হাজার হাজার lane চলে (SIMT — Single Instruction, Multiple Threads model, SM ও warp-এর গঠন সহ)। এই লেসনের AoS/SoA আর contiguous-memory আলোচনা সেখানেও ফিরে আসবে, কারণ GPU-র memory access pattern-এর প্রয়োজনীয়তা SIMD-এর মতোই, শুধু আরও কঠোর মাত্রায়।

আরও পড়ুন

  • Intel Intrinsics Guide · প্রতিটা AVX2 intrinsic-এর signature, latency, throughput-এর interactive রেফারেন্স
  • ARM NEON Intrinsics Reference — ARM Limited · vld1q_f32, vaddq_f32, vaddvq_f32-সহ AArch64 NEON intrinsics-এর প্রামাণ্য তালিকা
  • Data-Oriented Design — Richard Fabian · AoS বনাম SoA আর হার্ডওয়্যার-সচেতন data layout-এর একটা প্রামাণ্য, ব্যবহারিক আলোচনা
  • Software Optimization Resources — Agner Fog · x86 instruction latency/throughput টেবিল — intrinsics-এর প্রকৃত hardware cost বোঝার জন্য অপরিহার্য