SIMD আর vector — একসাথে অনেক ডেটা, SSE-AVX থেকে RISC-V Vector

১ · একসাথে অনেক ডেটা — data-level parallelism আর Flynn-এর শ্রেণি

ঈদের আগে একটি মিষ্টির দোকানের কথা ভাবুন। এক হাজার বাক্সে মিষ্টি ভরতে হবে, প্রতিটি বাক্সে একই কাজ — বাক্স খোলা, ছয়টি রসগোল্লা রাখা, ঢাকনা বন্ধ। একজন কর্মচারী একটি একটি করে বাক্স ভরলে সময় লাগবে অনেক। কিন্তু কাউন্টারে যদি আটটি বাক্স পাশাপাশি সাজিয়ে একটি চওড়া চামচ দিয়ে একবারে আটটি বাক্সে মিষ্টি দেওয়া যায়, তাহলে একই হাতের নড়াচড়ায় আট গুণ কাজ। শর্ত একটিই — প্রতিটি বাক্সে একই কাজ, আর এক বাক্সের কাজ অন্য বাক্সের ওপর নির্ভর করে না।

কম্পিউটারে এমন কাজ প্রচুর: একটি ছবির ২০ লক্ষ pixel-এর উজ্জ্বলতা বাড়ানো, একটি গানের প্রতিটি sample-এর volume কমানো, দুটি লম্বা array যোগ করা, AI model-এর লাখ লাখ গুণ-যোগ। প্রতিটি উপাদানে একই হিসাব, আর উপাদানগুলো একে অপরের থেকে স্বাধীন। এই সুযোগকে বলে data-level parallelism (DLP). আগের পর্বগুলোতে আমরা অন্য ধরনের সমান্তরালতা দেখেছি — pipeline-এ এক instruction-এর বিভিন্ন ধাপ একসাথে (পর্ব ৩৫), আর পর্ব ৩৯-এ প্রতি cycle-এ কয়েকটি স্বাধীন instruction একসাথে (instruction-level parallelism). DLP-র কৌশল ভিন্ন: instruction একটিই, কিন্তু সে অনেকগুলো ডেটায় কাজ করে।

Flynn-এর চার শ্রেণি

১৯৬৬ সালে কম্পিউটার বিজ্ঞানী Michael Flynn কম্পিউটারকে দুটি সরল প্রশ্নে ভাগ করার প্রস্তাব দেন: একসাথে কয়টি instruction-এর ধারা (stream) চলে, আর কয়টি ডেটার ধারা? উত্তর “এক” বা “অনেক” — তাই চারটি ঘর। এই ভাগকে বলে Flynn's taxonomy:

শ্রেণিগুলো পরিষ্কার মনে হলেও আজকের chip একসাথে কয়েকটি শ্রেণিতে পড়ে: একটি ৮-core ফোনের CPU সামগ্রিকভাবে MIMD, কিন্তু প্রতিটি core-এর ভেতরে SIMD unit আছে, আর প্রতিটি core নিজে SISD-এর মতো instruction চালায়। Flynn-এর ভাগ তাই একটি চিন্তার কাঠামো, কঠিন সীমারেখা নয়। পর্ব ৫৪-এ সমান্তরাল আর্কিটেকচার আলোচনায় এই শ্রেণিগুলোতে আবার ফিরব, NUMA আর interconnect-সহ। এই পর্বের কাজ শুধু SIMD ঘরটি — একটি core-এর ভেতরে একটি instruction কীভাবে অনেক ডেটা ছোঁয়।

২ · SIMD register আর lane — packed data, packed যোগ আর saturating হিসাব

একটি সাধারণ register-এ একটি সংখ্যা থাকে। একটি SIMD বা vector register আসলে একটি চওড়া bit-এর সারি — যেমন x86-এর XMM register 128 bit চওড়া। এই 128 bit-কে কীভাবে “পড়া” হবে, তা ঠিক করে instruction:

সূত্রটি সরল: lane সংখ্যা = register-এর চওড়া / উপাদানের চওড়া। পাশাপাশি গুঁজে রাখা এই মানগুলোকে বলে packed data. Lane-গুলো একে অপরের থেকে আলাদা — lane 0-র যোগে carry হলে তা lane 1-এ যায় না। পর্ব ১১-এর adder মনে করুন: একটি 128-bit adder-এর carry-শৃঙ্খল প্রতি 32 bit পরপর কেটে দিলেই চারটি স্বাধীন 32-bit adder পাওয়া যায়। Hardware-এ packed integer যোগ ঠিক এভাবেই সস্তায় বানানো হয় — একই তার, শুধু lane-এর সীমানায় carry আটকানোর একটি control signal.

এখন একটি packed float যোগ (x86-এ addps — add packed single) lane ধরে ধরে দেখুন। বোঝার সুবিধার জন্য stepper প্রতিটি lane আলাদা ধাপে দেখায়, কিন্তু মনে রাখবেন hardware চারটি lane একই সময়ে হিসাব করে:

Wrap-around বনাম saturating — ছবির উজ্জ্বলতা বাড়ানো

একটি grayscale ছবির প্রতিটি pixel একটি 8-bit unsigned সংখ্যা: 0 মানে কালো, 255 মানে সাদা। পুরো ছবি উজ্জ্বল করতে প্রতিটি pixel-এ 60 যোগ করি। সাধারণ 8-bit যোগে 200 + 60 = 260, যা 8 bit-এ ধরে না — রাখা হয় 260 − 256 = 4. এটি wrap-around (পর্ব ০৫)। ফল ভয়াবহ: প্রায় সাদা pixel হঠাৎ প্রায় কালো হয়ে যায়, ছবির উজ্জ্বল অংশে কালো দাগ পড়ে।

ছবির জন্য সঠিক আচরণ হলো saturating arithmetic — সীমা ছাড়ালে ফল সর্বোচ্চ মানে (255) আটকে থাকবে। Scalar কোডে এর জন্য প্রতিটি pixel-এ একটি তুলনা আর শর্ত লাগে। SIMD ISA-গুলো তাই আলাদা instruction দেয়: x86-এর paddusb (packed add unsigned saturate bytes), ARM NEON-এর uqadd, RISC-V Vector-এর vsaddu. একটি instruction-এ ১৬টি pixel, কোনো branch ছাড়া।

নিচের simulator একটি 128-bit register-কে JavaScript-এর typed array দিয়ে মডেল করে — একই ১৬ byte-এর ওপর Float32Array বসালে ৪টি lane, Uint8Array বসালে ১৬টি। Packed যোগ-গুণ, wrap বনাম saturate, একটি shuffle আর একটি horizontal যোগফল lane ধরে ধরে ছাপে। Run চাপুন:

// A 128-bit register = 16 bytes; view it as 4 x f32 or 16 x u8 lanes
const reg = () => new ArrayBuffer(16);
const show = (name, view) => console.log(name.padEnd(7) + "[" + Array.from(view).join(", ") + "]");
const a = reg(), b = reg(), c = reg();
const fa = new Float32Array(a), fb = new Float32Array(b), fc = new Float32Array(c);
fa.set([1.5, 2, 3, 4]); fb.set([10, 20, 30, 40]);
for (let i = 0; i < 4; i++) fc[i] = fa[i] + fb[i];          // addps: one op, 4 lanes
show("a", fa); show("b", fb); show("a+b", fc);
for (let i = 0; i < 4; i++) fc[i] = fa[i] * fb[i];          // mulps
show("a*b", fc);
// 16 x u8 pixels, brighten by +60
const p = new Uint8Array(16), wrap = new Uint8Array(16), sat = new Uint8Array(16);
p.set([0, 50, 100, 150, 180, 195, 196, 200, 210, 220, 230, 240, 250, 253, 254, 255]);
for (let i = 0; i < 16; i++) {
  wrap[i] = p[i] + 60;                                       // paddb: wraps mod 256
  sat[i] = Math.min(p[i] + 60, 255);                         // paddusb: saturates
}
show("pixel", p); show("wrap", wrap); show("sat", sat);
// shuffle: reverse 4 lanes (like _mm_shuffle_ps(x, x, 0x1B))
const sh = new Float32Array(4), ctl = [3, 2, 1, 0];
for (let i = 0; i < 4; i++) sh[i] = fa[ctl[i]];
show("rev(a)", sh);
// horizontal sum of 8 lanes as a tree: 8 -> 4 -> 2 -> 1
let v = [3, 1, 4, 1, 5, 9, 2, 6], step = 0;
while (v.length > 1) {
  const h = v.length / 2, next = [];
  for (let i = 0; i < h; i++) next.push(v[i] + v[i + h]);  // add upper half onto lower half
  v = next; step++;
  console.log("tree step " + step + ": [" + v.join(", ") + "]");
}
console.log("sum = " + v[0] + " in " + step + " vector adds (scalar needs 7)");

a      [1.5, 2, 3, 4]
b      [10, 20, 30, 40]
a+b    [11.5, 22, 33, 44]
a*b    [15, 40, 90, 160]
pixel  [0, 50, 100, 150, 180, 195, 196, 200, 210, 220, 230, 240, 250, 253, 254, 255]
wrap   [60, 110, 160, 210, 240, 255, 0, 4, 14, 24, 34, 44, 54, 57, 58, 59]
sat    [60, 110, 160, 210, 240, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255]
rev(a) [4, 3, 2, 1.5]
tree step 1: [8, 10, 6, 7]
tree step 2: [14, 17]
tree step 3: [31]
sum = 31 in 3 vector adds (scalar needs 7)

৩ · Vertical বনাম horizontal কাজ, shuffle আর reduction tree

এতক্ষণ যা দেখলাম — দুটি register-এর একই নম্বরের lane মিলিয়ে কাজ — তাকে বলে vertical কাজ। কাগজে দুই register ওপর-নিচে লিখলে প্রতিটি হিসাব খাড়া কলাম বরাবর হয়। SIMD hardware এতেই সবচেয়ে দক্ষ, কারণ প্রতিটি lane-এর তার সোজা নিজের adder-এ যায়।

কিন্তু অনেক সময় একই register-এর ভেতরের lane-গুলোকে একে অপরের সাথে মেলাতে হয় — যেমন সব lane যোগ করে একটি যোগফল (dot product-এর শেষ ধাপ), বা সব lane-এর মধ্যে সর্বোচ্চ মান। একে বলে horizontal operation. এর জন্য ডেটাকে lane-এর সীমানা পেরিয়ে সরাতে হয়, আর সেই কাজ করে shuffle (বা permute) instruction — যা register-এর lane-গুলোকে যেকোনো নতুন ক্রমে সাজাতে পারে: উল্টে দেওয়া, জোড়া অদলবদল, একটি মান সব lane-এ ছড়ানো (broadcast). x86-এ shufps, pshufd, pshufb; ARM-এ tbl, zip, uzp; RISC-V-এ vrgather, vslidedown.

ঠিক এই জোড়ায় জোড়ায় কৌশলের নাম reduction tree. প্রতিটি ধাপে register-এর ওপরের অর্ধেক shuffle করে নিচের অর্ধেকের ওপর এনে একটি vertical যোগ — ফলে প্রতিটি ধাপে বাকি মানের সংখ্যা অর্ধেক হয়:

শেষ কথাটি গুরুত্বপূর্ণ। একটি লম্বা array-এর যোগফলে দক্ষ SIMD কোড লুপের ভেতরে শুধু vertical যোগ করে — ৮টি আলাদা “আংশিক যোগফল” lane-এ জমায় — আর লুপ শেষে একবার reduction tree চালায়। Horizontal কাজ প্রতি ধাপে করলে লাভের বেশিরভাগ shuffle-এ খরচ হয়ে যায়।

৪ · x86-এর SIMD — MMX থেকে AVX-512, আর ARM NEON

Desktop প্রসেসরে SIMD ঢুকেছিল মাল্টিমিডিয়ার হাত ধরে — নব্বইয়ের দশকে কম্পিউটারে ভিডিও, গান আর 3D গেম চালানোর চাপ থেকে। প্রায় তিন দশকে x86-এ একের পর এক extension যোগ হয়েছে। বছরগুলো প্রথম chip বাজারে আসার মোটামুটি সময়:

লক্ষ করুন, প্রতিটি নতুন চওড়ার জন্য নতুন register নাম, নতুন instruction encoding আর নতুন intrinsic লাগে — একই “দুটি array যোগ” কাজের জন্য x86-এ এখন paddd, vpaddd xmm, vpaddd ymm, vpaddd zmm — চার রূপ। পুরোনো প্রোগ্রাম নতুন চওড়ার সুবিধা পায় না যতক্ষণ না আবার compile বা নতুন করে লেখা হয়। এই ক্রমাগত বাড়তে থাকা ISA-কে অনেকে বলেন “ISA churn” — অংশ ৬-এ দেখব vector ISA কীভাবে এ সমস্যা এড়ায়।

ARM NEON

ARM-এর SIMD-এর নাম NEON (আনুষ্ঠানিক নাম Advanced SIMD). এটি সবসময় 128 bit চওড়া; 64-bit ARM-এ (AArch64) ৩২টি V register (v0–v31). প্রায় প্রতিটি আধুনিক Android ফোন আর Apple-এর ARM-ভিত্তিক Mac-এ NEON আছে — আপনার ফোনে ক্যামেরার ছবি প্রসেস করা বা ভিডিও চালানোর সময় সম্ভবত NEON কাজ করছে। NEON-এর instruction-এর নামে উপাদানের টাইপ লেখা থাকে: add v0.4s, v1.4s, v2.4s মানে চারটি 32-bit (s = single word) lane যোগ, uqadd v0.16b, … মানে ১৬টি byte-এর saturating যোগ।

নিজের হাতে — SSE intrinsics দিয়ে C

C থেকে SIMD instruction ব্যবহারের সবচেয়ে সরাসরি পথ হলো intrinsic — দেখতে function-এর মতো, কিন্তু compiler প্রতিটিকে একটি নির্দিষ্ট instruction-এ বদলায়। x86-এর নামের নিয়ম: _mm = 128-bit (_mm256 = 256, _mm512 = 512), তারপর কাজের নাম, তারপর টাইপ — ps = packed single float, pd = packed double, epi8/epi32 = packed signed integer, epu8 = packed unsigned byte. তাই _mm_add_ps → addps, _mm_adds_epu8 → paddusb.

x86-64-এ SSE2 প্রতিটি CPU-তে থাকতেই হবে, তাই নিচের কোড যেকোনো x86-64 machine-এ চলবে। প্রোগ্রামটি দুটি ৪০৯৬-উপাদানের float array একবার scalar লুপে, একবার SSE দিয়ে (প্রতি ঘোরায় ৪টি) যোগ করে, ফল মিলিয়ে দেখে আর সময়ের অনুপাত ছাপে:

#define _POSIX_C_SOURCE 199309L
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <emmintrin.h>                    /* SSE/SSE2: baseline on every x86-64 CPU */

#define N 4096
#define REPS 20000
static float a[N], b[N], c1[N], c2[N];

static double now(void) {
    struct timespec t; clock_gettime(CLOCK_MONOTONIC, &t);
    return (double)t.tv_sec + (double)t.tv_nsec / 1e9;
}
int main(void) {
    for (int i = 0; i < N; i++) { a[i] = (float)i * 0.5f; b[i] = (float)(N - i); }
    double t0 = now();
    for (int r = 0; r < REPS; r++)
        for (int i = 0; i < N; i++) c1[i] = a[i] + b[i];          /* 1 float per add */
    double t1 = now();
    for (int r = 0; r < REPS; r++)
        for (int i = 0; i < N; i += 4) {                         /* 4 floats per add */
            __m128 va = _mm_loadu_ps(&a[i]), vb = _mm_loadu_ps(&b[i]);
            _mm_storeu_ps(&c2[i], _mm_add_ps(va, vb));
        }
    double t2 = now();
    printf("identical results: %s\n", memcmp(c1, c2, sizeof c1) == 0 ? "yes" : "NO");
    printf("c[0]=%.1f c[1]=%.1f c[4095]=%.1f\n", c1[0], c1[1], c1[N - 1]);
    printf("scalar loop iterations: %d, SSE loop iterations: %d\n", N, N / 4);
    printf("time ratio scalar/SSE: %.1f (varies per run)\n", (t1 - t0) / (t2 - t1));
    printf("avx2 supported here: %s\n", __builtin_cpu_supports("avx2") ? "yes" : "no");
    return 0;
}

identical results: yes
c[0]=4096.0 c[1]=4095.5 c[4095]=2048.5
scalar loop iterations: 4096, SSE loop iterations: 1024
time ratio scalar/SSE: 1.3 (varies per run)
avx2 supported here: yes

এবার অংশ ২-এর pixel উদাহরণটি আসল hardware-এ। _mm_set1_epi8(60) ১৬টি lane-এ 60 বসায় (broadcast), আর _mm_adds_epu8 একটি instruction-এ ১৬টি saturating যোগ করে। ফলটি scalar “তুলনা করে আটকানো” কোডের সাথে মিলিয়ে দেখা হয়:

#include <stdio.h>
#include <stdint.h>
#include <emmintrin.h>                     /* SSE2: _mm_adds_epu8 = PADDUSB */

int main(void) {
    uint8_t px[16] = {0, 50, 100, 150, 180, 195, 196, 200,
                      210, 220, 230, 240, 250, 253, 254, 255};
    uint8_t scal[16], simd[16], wrap[16];
    for (int i = 0; i < 16; i++) {                   /* scalar: compare and clamp */
        unsigned s = (unsigned)px[i] + 60u;
        scal[i] = (uint8_t)(s > 255u ? 255u : s);
        wrap[i] = (uint8_t)s;                        /* plain add keeps only 8 bits */
    }
    __m128i v = _mm_loadu_si128((const __m128i *)px);
    __m128i k = _mm_set1_epi8(60);                   /* 60 in all 16 lanes */
    _mm_storeu_si128((__m128i *)simd, _mm_adds_epu8(v, k));   /* 16 lanes, 1 instr */
    int same = 1;
    for (int i = 0; i < 16; i++) same &= (scal[i] == simd[i]);
    printf("in  :"); for (int i = 0; i < 16; i++) printf(" %3u", px[i]);   printf("\n");
    printf("wrap:"); for (int i = 0; i < 16; i++) printf(" %3u", wrap[i]); printf("\n");
    printf("sat :"); for (int i = 0; i < 16; i++) printf(" %3u", simd[i]); printf("\n");
    printf("scalar clamp == _mm_adds_epu8: %s\n", same ? "yes" : "NO");
    return 0;
}

in  :   0  50 100 150 180 195 196 200 210 220 230 240 250 253 254 255
wrap:  60 110 160 210 240 255   0   4  14  24  34  44  54  57  58  59
sat :  60 110 160 210 240 255 255 255 255 255 255 255 255 255 255 255
scalar clamp == _mm_adds_epu8: yes

৫ · ক্লাসিক vector machine — Cray-1, chaining, stride, gather-scatter আর mask

SIMD-এর ধারণা desktop-এর চেয়ে অনেক পুরোনো। সত্তরের দশকের supercomputer-গুলো ছিল vector processor — যাদের instruction পুরো একটি সংখ্যার সারি (vector) নিয়ে কাজ করত।

Cray-1-এর মতো vector machine-এর কয়েকটি ধারণা আজকের RISC-V Vector-এ প্রায় হুবহু ফিরে এসেছে, তাই সেগুলো ভালো করে জানা দরকার:

• Vector length register (VL). Register-এ ৬৪টি জায়গা, কিন্তু কাজ হয়তো ৫০টি উপাদানে। প্রোগ্রাম VL-এ ৫০ লিখে দিলে পরের সব vector instruction শুধু প্রথম ৫০টি উপাদানে কাজ করে। ৬৪-এর বেশি উপাদান হলে লুপ ৬৪-র টুকরায় ভাগ হয় — এই কৌশলের নাম strip-mining (খনিতে স্তরে স্তরে খোঁড়ার মতো)। এটাই vector length-এর ধারণা।
• Pipelined functional unit, একটি একটি করে উপাদান। Packed SIMD-এর মতো ৬৪টি adder পাশাপাশি থাকত না; একটি গভীর pipelined adder প্রতি cycle-এ একটি নতুন উপাদান নিত। একটি instruction ৬৪ cycle ধরে চলত, কিন্তু fetch-decode হতো একবার।
• Chaining. V3 = V1 × V2 তারপর V5 = V3 + V4 — দ্বিতীয় instruction-কে পুরো V3-এর জন্য ৬৪ cycle অপেক্ষা করতে হয় না; গুণের প্রথম উপাদান বের হতেই যোগ শুরু। এটি পর্ব ৩৬-এর forwarding-এর vector রূপ, এর ফলে গুণ আর যোগ একসাথে চলে, প্রতি cycle-এ দুটি floating point কাজ হতো।
• Strided access. Row-major matrix-এর একটি column পড়তে প্রতি উপাদানে একটি পুরো row সমান লাফ দিতে হয়। Vector load-এ একটি stride দিলেই হতো — একটি instruction-এ “প্রতি N-তম উপাদান”।
• Gather / scatter. Sparse matrix-এ উপাদান এলোমেলো জায়গায়। একটি index-vector ধরে মেমরির নানা জায়গা থেকে জোগাড় (gather) বা ছড়িয়ে লেখা (scatter) — পরের প্রজন্মের vector machine-এ আসে, আর আজ AVX2/AVX-512, SVE আর RVV-তে আছে। তবে এটি সাধারণ টানা load-এর চেয়ে অনেক ধীর, কারণ প্রতিটি উপাদান আলাদা cache line-এ থাকতে পারে।
• Vector mask. “যেখানে a[i] > 0, সেখানেই b[i] = a[i] × 2” — প্রতি lane-এ একটি bit বলে দেয় সেই lane-এ ফল লেখা হবে কি না। Cray-1-এ এর জন্য একটি vector mask register ছিল। এটি পর্ব ২৯-এর predication-এর vector রূপ: branch-এর বদলে দুই দিকের হিসাব করে mask দিয়ে বাছাই।

আশি-নব্বইয়ের দশকে সস্তা microprocessor-এর অনেকগুলো একসাথে জুড়ে বানানো MIMD machine (cluster) বিশেষায়িত vector supercomputer-কে বাজার থেকে প্রায় সরিয়ে দেয়। কিন্তু ধারণাগুলো হারায়নি — প্রথমে ছোট আকারে desktop-এর packed SIMD-এ, পরে GPU-তে, আর এখন আবার প্রায় পূর্ণ রূপে ARM SVE আর RISC-V Vector-এ।

৬ · RISC-V Vector — vsetvli, VLEN, SEW, LMUL আর strip-mining

RISC-V-এর vector extension (RVV, ISA-র নামে অক্ষর “V”)-এর version 1.0 ২০২১ সালের শেষ দিকে অনুমোদিত (ratified) হয়। এর মূল সিদ্ধান্ত Cray-ধারার: register কত চওড়া, তা ISA বলে না — chip ঠিক করে। প্রোগ্রাম চলার সময় জেনে নেয়। চারটি নাম মনে রাখলেই পুরো নকশা বোঝা যায়:

লুপের প্রতিটি ধাপের শুরুতে vsetvli instruction-কে বলা হয় “আমার আরও AVL (application vector length) সংখ্যক উপাদান বাকি, উপাদান 32-bit” — সে vl ঠিক করে একটি register-এ ফেরত দেয়। প্রচলিত বাস্তবায়নে vl = min(AVL, VLMAX). নিচে C-এর for (i = 0; i < n; i++) c[i] = a[i] + b[i]; (32-bit int)-এর একটি হাতে লেখা RVV লুপ — -O2 -march=rv64gcv-এ GCC বা Clang যে ধরনের লুপ বানায় তার কাছাকাছি:

# উদাহরণ (চালানো হয়নি) — RVV 1.0, vector add: c[i] = a[i] + b[i], 32-bit int
# a0 = n, a1 = &a[0], a2 = &b[0], a3 = &c[0]
vadd:
    vsetvli t0, a0, e32, m1, ta, ma   # t0 = vl = min(n, VLMAX); SEW=32, LMUL=1
    vle32.v v0, (a1)                  # vl ints from a  -> v0
    vle32.v v1, (a2)                  # vl ints from b  -> v1
    vadd.vv v2, v0, v1                # vl lane-wise adds
    vse32.v v2, (a3)                  # v2 -> vl ints of c
    sub     a0, a0, t0                # n -= vl
    slli    t1, t0, 2                 # bytes done = vl * 4
    add     a1, a1, t1                # advance the three pointers
    add     a2, a2, t1
    add     a3, a3, t1
    bnez    a0, vadd                  # elements left? loop again
    ret

প্রতিটি লাইনের অর্থ:

• vsetvli t0, a0, e32, m1, ta, ma — বাকি উপাদান a0, উপাদান 32-bit (e32), একটি register-এর দল (m1). ta (tail agnostic) মানে vl-এর পরের lane-গুলোতে কী থাকবে তা নিয়ে প্রোগ্রামের মাথাব্যথা নেই, ma (mask agnostic) mask-বন্ধ lane-এর জন্য একই কথা — এতে hardware-এর কাজ সহজ হয়। ফল vl যায় t0-তে।
• vle32.v — unit-stride vector load: মেমরি থেকে পরপর vlটি 32-bit উপাদান। (Strided load vlse32.v আর gather vluxei32.v-ও আছে — Cray-এর সেই ধারণাগুলো।)
• vadd.vv — vector-vector যোগ, শুধু প্রথম vlটি উপাদানে। (.vx রূপে একটি scalar register সব lane-এ যোগ হয়, .vi-তে একটি ছোট immediate.)
• vse32.v — vector store, ঠিক vlটি উপাদান; তার পরের মেমরি ছোঁয়া হয় না।
• sub, slli, তিনটি add, bnez — সাধারণ RV64I scalar instruction (পর্ব ২৬): বাকি সংখ্যা কমানো, pointer এগোনো (vl × ৪ byte), আর বাকি থাকলে আবার লুপ।

একই কাজের scalar লুপে প্রতি উপাদানে ৯টি instruction লাগে (দুটি lw, add, sw, তিনটি pointer-এর addi, counter-এর addi, bnez). Saxpy (y[i] += a × x[i])-এর float রূপেও কাঠামো একই — শুধু vadd.vv-এর জায়গায় vfmacc.vf v1, fa0, v0 (প্রতি lane-এ v1 += fa0 × v0).

এবার n = ১৩ আর VLEN = 128 (VLMAX = ৪) ধরে লুপটি ধাপে ধাপে চালাই। শেষ ধাপে vl কীভাবে বদলায়, খেয়াল করুন:

নিচের model কয়েকটি n আর VLEN/SEW/LMUL-এর জন্য প্রতিটি vsetvli-এর vl, ঘোরার সংখ্যা আর মোট instruction (প্রতি ঘোরায় ১১টি) বনাম scalar (প্রতি উপাদানে ৯টি) হিসাব করে। প্রথম দুটি লাইনে একই binary দুটি ভিন্ন chip-এ চলছে — খেয়াল করুন:

// RVV strip-mining model of c[i] = a[i] + b[i] (32-bit ints)
// VLMAX = VLEN * LMUL / SEW; vsetvli returns vl = min(AVL, VLMAX) (common choice)
function stripMine(n, VLEN, SEW, LMUL) {
  const VLMAX = VLEN * LMUL / SEW;
  console.log(`n=${n} VLEN=${VLEN} SEW=${SEW} LMUL=${LMUL} -> VLMAX=${VLMAX}`);
  let avl = n, iter = 0, done = 0;
  const total = Math.ceil(n / VLMAX);
  while (avl > 0) {
    const vl = Math.min(avl, VLMAX);           // what vsetvli writes to t0
    iter++;
    if (total <= 5 || iter <= 2 || iter === total)
      console.log(`  iter ${iter}: avl=${avl} -> vl=${vl}, elements ${done}..${done + vl - 1}`);
    else if (iter === 3) console.log("  ...");
    done += vl; avl -= vl;
  }
  const vec = iter * 11;    // vsetvli, 2 vle32, vadd, vse32, sub, slli, 3 add, bnez
  const scalar = n * 9;     // 2 lw, add, sw, 3 addi (pointers), addi (count), bnez
  console.log(`  ${iter} iterations: ~${vec} instr vs ~${scalar} scalar instr` +
              ` (${(scalar / vec).toFixed(1)}x fewer)`);
}
stripMine(13, 128, 32, 1);     // same binary...
stripMine(13, 256, 32, 1);     // ...on a CPU with wider registers
stripMine(1000, 128, 32, 1);
stripMine(1000, 128, 32, 4);   // LMUL=4: a group of 4 registers
stripMine(1000, 512, 32, 8);

n=13 VLEN=128 SEW=32 LMUL=1 -> VLMAX=4
  iter 1: avl=13 -> vl=4, elements 0..3
  iter 2: avl=9 -> vl=4, elements 4..7
  iter 3: avl=5 -> vl=4, elements 8..11
  iter 4: avl=1 -> vl=1, elements 12..12
  4 iterations: ~44 instr vs ~117 scalar instr (2.7x fewer)
n=13 VLEN=256 SEW=32 LMUL=1 -> VLMAX=8
  iter 1: avl=13 -> vl=8, elements 0..7
  iter 2: avl=5 -> vl=5, elements 8..12
  2 iterations: ~22 instr vs ~117 scalar instr (5.3x fewer)
n=1000 VLEN=128 SEW=32 LMUL=1 -> VLMAX=4
  iter 1: avl=1000 -> vl=4, elements 0..3
  iter 2: avl=996 -> vl=4, elements 4..7
  ...
  iter 250: avl=4 -> vl=4, elements 996..999
  250 iterations: ~2750 instr vs ~9000 scalar instr (3.3x fewer)
n=1000 VLEN=128 SEW=32 LMUL=4 -> VLMAX=16
  iter 1: avl=1000 -> vl=16, elements 0..15
  iter 2: avl=984 -> vl=16, elements 16..31
  ...
  iter 63: avl=8 -> vl=8, elements 992..999
  63 iterations: ~693 instr vs ~9000 scalar instr (13.0x fewer)
n=1000 VLEN=512 SEW=32 LMUL=8 -> VLMAX=128
  iter 1: avl=1000 -> vl=128, elements 0..127
  iter 2: avl=872 -> vl=128, elements 128..255
  ...
  iter 8: avl=104 -> vl=104, elements 896..999
  8 iterations: ~88 instr vs ~9000 scalar instr (102.3x fewer)

Packed SIMD বনাম vector-length agnostic নকশা

উপরের কোডকে বলে vector-length agnostic (VLA) — সে জানে না register কত চওড়া। ARM-ও একই পথে গেছে: SVE (২০১৬ সালে ঘোষিত) আর Armv9-এর SVE2-তে register-এর চওড়া chip ভেদে 128 থেকে 2048 bit (128-এর গুণিতকে) হতে পারে। SVE tail সামলায় predicate register দিয়ে: whilelo p0.s, x_i, x_n instruction প্রতিটি lane-এর জন্য “i + lane < n কি না” হিসাব করে একটি mask বানায়, আর সব load-store-যোগ সেই mask মেনে চলে। Fujitsu-র A64FX chip (জাপানের Fugaku supercomputer, যেটি ২০২০ সালে বিশ্বের দ্রুততম তালিকার শীর্ষে ছিল) 512-bit SVE ব্যবহার করে বলে জানা যায়।

৭ · Compiler আর SIMD — auto-vectorization, intrinsics আর portable পথ

বেশিরভাগ প্রোগ্রামার কখনো addps নিজে লেখেন না। আধুনিক compiler একটি সাধারণ C লুপকে নিজেই SIMD instruction-এ বদলাতে পারে — একে বলে auto-vectorization. GCC-তে এটি -O3-এ পুরোপুরি চালু, আর GCC 12 থেকে -O2-তেও একটি সাবধানী রূপে; Clang -O2-তেই করে। কিন্তু compiler-কে প্রমাণ করতে হয় যে রূপান্তরটি প্রোগ্রামের অর্থ বদলাবে না। যেসব শর্ত না মিললে compiler পিছিয়ে যায়:

1. Trip count আগে থেকে গোনা যায়। for (i = 0; i < n; i++) ভালো — লুপ শুরুর সময়েই জানা কতবার ঘুরবে। while (*p != 0) বা ভেতরে break থাকলে কঠিন (নতুন কিছু GCC সংস্করণ এমন কিছু “early-exit” লুপও vectorize করতে শিখছে, তবে সব ক্ষেত্রে নয়)।
2. Aliasing নেই। void add(float *x, float *y, float *z, int n)-এ compiler জানে না z আর x একই মেমরিতে ওভারল্যাপ করে কি না (memory aliasing). যদি z = x + 1 হয়, তাহলে এক ঘোরার লেখা পরের ঘোরার পড়াকে বদলায় — একসাথে চারটি পড়া ভুল ফল দেবে। Compiler তখন হয় চলার সময় ওভারল্যাপ যাচাইয়ের একটি পরীক্ষা বসায় (দুটি সংস্করণ, “loop versioning”), নয়তো vectorize-ই করে না। প্রোগ্রামার C99-এর restrict দিয়ে প্রতিশ্রুতি দিতে পারেন: float *restrict z — “এই pointer-এর মেমরি আর কেউ ছোঁবে না”। প্রতিশ্রুতি মিথ্যা হলে আচরণ undefined, তাই সাবধানে।
3. লুপের এক ঘোরা থেকে পরের ঘোরায় নির্ভরতা নেই। a[i] = a[i-1] + b[i] (prefix sum)-এ প্রতিটি উপাদান আগেরটির ফলের অপেক্ষায় — এটি loop-carried dependence, সরাসরি vectorize হয় না (বিশেষ শিফট-যোগ কৌশল লাগে)। a[i] = a[i+4] + 1-এর মতো দূরের নির্ভরতা হলে ৪টি lane পর্যন্ত নিরাপদ — compiler দূরত্ব হিসাব করে।
4. Float-এর reduction-এ ক্রম বদলানোর অনুমতি। s += x[i] integer-এ সহজে vectorize হয়। Float-এ হয় না, কারণ অংশ ৩-এর মতো ক্রম বদলালে ফল বদলাতে পারে, আর C মান কঠোর ক্রম চায়। -ffast-math (বা শুধু -fassociative-math) দিলে compiler ক্রম বদলাতে পারে। কিন্তু -ffast-math আরও অনেক কিছু ধরে নেয় — যেমন NaN আর infinity কখনো আসবে না — ফলে isnan()-এর পরীক্ষা পর্যন্ত মুছে যেতে পারে। নিরাপদ বিকল্প: শুধু ওই লুপে #pragma omp simd reduction(+:s) (OpenMP), যা স্পষ্টভাবে বলে এই যোগফলের ক্রম বদলানো চলবে।
5. Alignment আর tail. পুরোনো SSE-তে aligned load (movaps) unaligned-এর চেয়ে অনেক দ্রুত ছিল; আধুনিক x86-এ পার্থক্য সাধারণত সামান্য, বিশেষত cache line-এর সীমানা না পেরোলে। Compiler শুরুতে কয়েকটি উপাদান আলাদা করে (peeling) বাকিটা aligned করে, আর শেষে বাকি থাকা উপাদানের জন্য একটি scalar “epilogue” লুপ বা (AVX-512, SVE, RVV-তে) mask বসায়।
6. লুপের ভেতরে function call বা জটিল শাখা নেই। সরল if হলে compiler দুই দিক হিসাব করে mask দিয়ে মেলাতে পারে (if-conversion); কিন্তু অজানা function call বা গভীর শাখা থাকলে সাধারণত হাল ছেড়ে দেয়।

Compiler কী করল তা অনুমান না করে জিজ্ঞেস করুন: GCC-তে -fopt-info-vec-optimized (সফল) আর -fopt-info-vec-missed (কেন পারেনি), Clang-এ -Rpass=loop-vectorize আর -Rpass-missed=loop-vectorize. Compiler Explorer-এ assembly দেখুন — addps/vaddps-এর মতো packed (ps) instruction বা ymm/zmm register দেখলে vectorize হয়েছে; শুধু xmm দেখে নয়, কারণ scalar float-এর addss-ও xmm ব্যবহার করে (পর্ব ৩১)। পরিচিত একটি নমুনা — GCC -O3-এ x86-64-এর জন্য float array যোগের মূল লুপ মোটামুটি এমন দেখায়:

# উদাহরণ — gcc -O3 -masm=intel (x86-64, SSE2, Intel syntax), c[i] = a[i] + b[i] এর ভেতরের লুপ; সংস্করণভেদে কিছুটা বদলায়
.L3:
    movups  xmm0, XMMWORD PTR [rsi+rax]     # 4 floats of b
    movups  xmm1, XMMWORD PTR [rdi+rax]     # 4 floats of a
    addps   xmm0, xmm1                      # 4 adds
    movups  XMMWORD PTR [rdx+rax], xmm0     # 4 floats to c
    add     rax, 16
    cmp     rax, rcx
    jne     .L3

Intrinsics বনাম portable পথ

Compiler না পারলে তিনটি পথ খোলা থাকে:

• Intrinsics (_mm_add_ps, NEON-এর vaddq_f32, RVV-এর __riscv_vadd_vv_i32m1) — পূর্ণ নিয়ন্ত্রণ, কিন্তু কোড এক ISA-তে বাঁধা; x86, ARM আর RISC-V-এর জন্য তিনবার লিখতে হয়, আর AVX-512 থাকলে-না-থাকলে চলার সময় বেছে নেওয়ার (runtime dispatch) কোডও লাগে।
• Portable SIMD library বা ভাষার সুবিধা — GCC/Clang-এর vector extension (typedef float v4sf __attribute__((vector_size(16)));), OpenMP-এর #pragma omp simd, C++-এর std::experimental::simd (C++26-এ std::simd হিসেবে মানে আসার কথা), Google-এর Highway library বা Intel-এর ISPC. একবার লিখে বহু ISA-তে compile.
• তৈরি library — BLAS, FFmpeg, OpenCV, libjpeg-turbo, glibc-এর memchr/strlen: বিশেষজ্ঞরা প্রতিটি ISA-র জন্য হাতে টিউন করে রেখেছেন। বাস্তবে সবচেয়ে বড় লাভ প্রায়ই এখান থেকেই আসে।

৮ · কখন SIMD কাজে লাগে, কখন লাগে না — আর Amdahl-এর ছাদ

যেখানে SIMD ঝলমল করে

• ছবি আর ভিডিও: প্রতিটি pixel-এ একই হিসাব — উজ্জ্বলতা, রং বদল, ছোট-বড় করা, video codec-এর block হিসাব। ভিডিও decoder-এ হাতে লেখা SIMD assembly খুব প্রচলিত।
• অডিও: sample-এর ওপর filter, volume, mixing — 16-bit বা float lane.
• AI আর বৈজ্ঞানিক হিসাব: matrix গুণ আর dot product মূলত গুণ-যোগ (FMA). ছোট টাইপ (int8, bfloat16) ব্যবহার করলে একই register-এ আরও বেশি lane, তাই আধুনিক CPU-তে AI-এর জন্য বিশেষ dot-product instruction যোগ হয়েছে।
• Cryptography আর checksum: ChaCha20-এর মতো cipher আর অনেক hash/CRC হিসাব SIMD-এ দ্রুত; x86-এর AES-NI instruction-ও XMM register ব্যবহার করে।
• String-এর কাজ: memchr (একটি byte খোঁজা) — ১৬টি byte একসাথে তুলনা (pcmpeqb), ফলকে একটি 16-bit bitmask-এ আনা (pmovmskb), আর mask শূন্য না হলে প্রথম 1-bit-এর অবস্থান (__builtin_ctz) — একবারে ১৬টি অক্ষর পরীক্ষা, কোনো প্রতি-byte branch ছাড়া। simdjson-এর মতো library এই কৌশলে JSON parsing-কে খুব দ্রুত করেছে বলে পরিচিত।

# উদাহরণ (চালানো হয়নি) — SSE2 দিয়ে memchr-এর মূল ধারণা, 16 byte একসাথে
__m128i key  = _mm_set1_epi8((char)c);                 // c in all 16 lanes
__m128i blk  = _mm_loadu_si128((const __m128i *)p);    // 16 bytes of text
int     mask = _mm_movemask_epi8(_mm_cmpeq_epi8(blk, key)); // bit k = 1 if p[k] == c
if (mask != 0) return p + __builtin_ctz(mask);         // first match
// else p += 16 and repeat; the last < 16 bytes need care (tail)

যেখানে SIMD প্রায় অকেজো

• Branch-ভরা কোড: প্রতিটি উপাদানে ভিন্ন পথ (parser-এর state machine, compiler নিজে, ব্যবসায়িক নিয়ম)। Mask দিয়ে দুই দিক হিসাব করা যায়, কিন্তু শাখা অনেক হলে বেশিরভাগ lane ফেলে দেওয়া কাজ করে।
• Pointer chasing: linked list বা tree ধরে হাঁটা — পরের ঠিকানা আগের load-এর ফলের ওপর নির্ভর, তাই একসাথে আটটি পড়ার উপায় নেই।
• ছোট n: ৩টি উপাদানের জন্য setup আর tail-এর খরচ লাভকে ছাড়িয়ে যায়।
• মেমরি-সীমাবদ্ধ (memory-bound) লুপ: দুটি বিশাল array যোগে প্রতি উপাদানে হিসাব মাত্র একটি, কিন্তু মেমরি থেকে ১২ byte আনা-নেওয়া। Cache-এ না ধরলে DRAM-এর bandwidth-ই সীমা; SIMD instruction কমায়, কিন্তু ডেটা আনার গতি বাড়ায় না। Cache-বান্ধব লুপের ক্রম আর blocking পর্ব ৪৮-এ।

Amdahl-এর ছাদ

পর্ব ০৩-এর Amdahl's law এখানে হুবহু খাটে। প্রোগ্রামের সময়ের f অংশ vectorize হলো, আর সেই অংশ s গুণ দ্রুত হলো:

ধরা যাক একটি ছবি-প্রসেসিং প্রোগ্রামের ৮০% সময় pixel লুপে (f = ০.৮)। ৪ lane-এর SSE সেই অংশকে আদর্শভাবে ৪ গুণ দ্রুত করলে মোট speedup = ১ / (০.২ + ০.৮/৪) = ১ / ০.৪ = ২.৫. ১৬ lane-এর AVX-512-এ ১ / (০.২ + ০.০৫) = ৪.০ — lane চার গুণ বেড়েছে, speedup মাত্র ১.৬ গুণ। অসীম lane-এও ছাদ ১ / ০.২ = ৫. আর বাস্তবে s কখনো পুরো lane-সংখ্যা হয় না — load-store, shuffle, tail আর মেমরির গতি মিলে সাধারণত কম। তাই বড় লাভের জন্য vectorize করা অংশ বাড়ানোই (f) আসল কাজ, শুধু চওড়া register নয়।

JavaScript আর ওয়েবে SIMD

এই পাতার simulator-গুলো কেন JavaScript-এ “lane মডেল” করে, আসল SIMD চালায় না? কারণ JavaScript-এ সরাসরি SIMD-এর কোনো API নেই। একসময় SIMD.js নামে একটি প্রস্তাব অনেক দূর এগিয়েছিল, কিন্তু ২০১৭-১৮ নাগাদ সেটি সরিয়ে নিয়ে মনোযোগ দেওয়া হয় WebAssembly SIMD-এর দিকে — একটি নির্দিষ্ট 128-bit v128 টাইপ, যা Chrome আর Firefox-এ ২০২১ সালে আর Safari-তে তার কিছু পরে চালু হয়। Browser চলার সময় সেটিকে CPU-র SSE বা NEON instruction-এ বদলায়। সাধারণ JavaScript-এর লুপ JIT compiler সাধারণত vectorize করে না বলেই জানা যায়, তাই আমাদের typed array-এর লুপ হলো শেখার মডেল, দ্রুততার কৌশল নয়।

৯ · AI-এর যুগে — SIMD কোড AI দিয়ে, যাচাই নিজে

“এই লুপটি AVX2 দিয়ে দ্রুত করে দিন” — AI এমন অনুরোধে সাথে সাথে intrinsics-ভরা কোড দেয়। দেখতে বিশ্বাসযোগ্য, কিন্তু SIMD কোডের ভুল সবচেয়ে চুপচাপ ভুল — compile হয়, চলে, শুধু কিছু উপাদানে ভুল ফল দেয়। এই বিষয়ে AI-এর সাধারণ ভুলগুলো:

• Tail ভুলে যাওয়া — n চারের গুণিতক ধরে নেওয়া, ফলে শেষের ১–৩টি উপাদান বাদ পড়ে বা array-এর বাইরে পড়া-লেখা হয়।
• Intrinsic-এর নাম বা argument-এর ক্রম বানানো — বিশেষত _mm_shuffle_ps-এর control byte, _mm_set_ps (উল্টো ক্রম!) বনাম _mm_setr_ps, signed বনাম unsigned saturating (epi8 বনাম epu8).
• CPU-তে নেই এমন extension ব্যবহার — AVX-512 কোড দিয়ে বলা “সব আধুনিক CPU-তে চলবে”, বা runtime যাচাই ছাড়া AVX2.
• Lane সংখ্যা, VLMAX বা strip-mining-এর ঘোরার সংখ্যা ভুল হিসাব; RVV-র পুরোনো 0.7 খসড়ার syntax (যেমন vsetvli-এ ta, ma ছাড়া, বা পুরোনো load নাম) আর 1.0 মিশিয়ে ফেলা।
• Float reduction-এ ক্রম বদলের প্রভাব উপেক্ষা করা, বা অন্ধভাবে -ffast-math সুপারিশ করা; আর chip-এর SIMD চওড়া বা বছর বানিয়ে বলা।

AI-এর দেওয়া SIMD কোড যাচাইয়ের চেকলিস্ট

• Scalar রূপের সাথে ফল মিলিয়েছি — বিশেষত ছোট n (0, 1, 3, 5) আর চওড়ার গুণিতক নয় এমন n-এ
• প্রতিটি intrinsic-এর নাম আর argument-এর ক্রম Intel Intrinsics Guide, Arm-এর NEON reference বা RISC-V Vector spec-এ মিলিয়েছি
• লক্ষ্য CPU-তে extension আছে কি না যাচাই করেছি (/proc/cpuinfo-এর flag, __builtin_cpu_supports) — নেই এমন instruction-এ প্রোগ্রাম crash করে
• Lane সংখ্যা, VLMAX আর ঘোরার সংখ্যা নিজে হিসাব করেছি (চওড়া / উপাদানের চওড়া; VLEN × LMUL / SEW)
• Optimisation flag ঠিক রেখে মেপেছি — -O0-এর সংখ্যা দিয়ে রায় দিইনি, আর compiler আগেই vectorize করেছে কি না -fopt-info-vec দিয়ে দেখেছি