Foundationপ্রথম নীতি থেকে
LEVEL 11অ্যাডভান্সড~১২ ঘণ্টাC

Matrix Multiply Optimization

Matrix Multiply Optimization

একই matrix multiply চারবার লেখা — naive, loop-reordered, tiled, AVX2 SIMD, তারপর multi-threaded — প্রতি ধাপে GFLOPS মেপে roofline model-এর তত্ত্বকে হাতেকলমে প্রমাণ করা।

মাইলস্টোন

আগে যা পড়া দরকার

কেন এই প্রজেক্ট

Roofline model লেসনে আমরা দেখেছি একটা algorithm-এর সর্বোচ্চ সম্ভাব্য গতি নির্ভর করে সেটা compute-bound না memory-bound তার উপর, আর SIMD programming লেসনে দেখেছি একটা ভেক্টর ইনস্ট্রাকশন একসাথে ৮টা float-এ কাজ করতে পারে। কিন্তু এই দুইটা তত্ত্বকে একসাথে জোড়া লাগিয়ে নিজের চোখে GFLOPS-এর একটা ১০-৫০× লাফ দেখা — সেটা সম্পূর্ণ ভিন্ন অভিজ্ঞতা।

Matrix multiply এই পরীক্ষার জন্য নিখুঁত, কারণ এর গাণিতিক সংজ্ঞা এত সহজ (তিনটা nested loop) যে “সঠিক” কোড লিখতে ৫ মিনিট লাগে — অথচ এই একই গণনা cache-friendly ভাবে সাজালে, vectorize করলে, আর একাধিক কোরে ছড়িয়ে দিলে পারফরম্যান্স মাত্রার-পর-মাত্রা বদলে যায়, যদিও প্রতিবার ঠিক একই সংখ্যক গাণিতিক অপারেশন হচ্ছে। এই প্রজেক্টে আমরা একই C = A × B-কে চারবার লিখব:

  1. Naive — সরল ত্রিনেস্টেড লুপ, সঠিক কিন্তু cache-এর সাথে বিবাদে
  2. Loop-reordered (tiled) — cache hierarchy-কে মাথায় রেখে সাজানো
  3. SIMD — AVX2 FMA দিয়ে একসাথে ৮টা float
  4. Multi-threaded — একাধিক কোরে row ভাগ করে

আর প্রতিটা ধাপে GFLOPS (প্রতি সেকেন্ডে বিলিয়ন floating-point operation) মেপে দেখব কোন অপ্টিমাইজেশন কতটা লাভ দেয়, আর কেন।

GFLOPS মাপার পদ্ধতি

N×N দুইটা matrix গুণ করতে মোট floating-point অপারেশন লাগে 2·N³ — প্রতিটা আউটপুট এলিমেন্টে N-টা multiply আর N-টা add (মোট N² এলিমেন্ট, তাই 2·N²·N = 2N³)।

#include <time.h>

static double now_seconds(void) {
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return ts.tv_sec + ts.tv_nsec / 1e9;
}

// fn চালিয়ে GFLOPS রিপোর্ট করা — সব ভার্সনের জন্য একই পরিমাপ পদ্ধতি
static double measure_gflops(const char *label, void (*fn)(float*, float*, float*, int),
                              float *A, float *B, float *C, int n) {
    double t0 = now_seconds();
    fn(A, B, C, n);
    double t1 = now_seconds();
    double flops = 2.0 * (double)n * (double)n * (double)n;
    double gflops = flops / (t1 - t0) / 1e9;
    printf("%-20s  %8.3f s   %8.3f GFLOPS\n", label, t1 - t0, gflops);
    return gflops;
}

N = 512 দিয়ে দ্রুত সঠিকতা-যাচাই, তারপর N = 1024 বা 2048 দিয়ে আসল GFLOPS বেঞ্চমার্ক করুন — ছোট N-এ loop overhead ও timer resolution ফলাফল দূষিত করতে পারে।

ধাপে ধাপে

১. Naive — ijk অর্ডার

void matmul_naive(float *A, float *B, float *C, int n) {
    for (int i = 0; i \< n; i++)
        for (int j = 0; j \< n; j++) {
            float sum = 0.0f;
            for (int k = 0; k \< n; k++)
                sum += A[i*n + k] * B[k*n + j];
            C[i*n + j] = sum;
        }
}

সঠিক, কিন্তু cache-এর সাথে খারাপ ব্যবহার করে — কেন সেটা বুঝতে হলে innermost k-loop-এ মেমরি অ্যাক্সেস প্যাটার্ন দেখতে হবে। প্রতি k-এ A[i*n+k] sequential (ভালো), কিন্তু B[k*n+j]-এ k বাড়লে ঠিকানা n floats করে লাফায় — এটা column-wise access, প্রতিটা নতুন k-তে সম্ভবত সম্পূর্ণ নতুন cache line, আর সেই line-এর বাকি ৬৩টা বাইটের মধ্যে মাত্র ৪টা বাইট (একটা float) ব্যবহার হয়ে বাকিটা নষ্ট হয়। বড় n-এ (matrix cache-এর চেয়ে বড়) প্রতিটা B অ্যাক্সেসই কার্যত একটা cache miss।

২. ikj — শুধু loop-অর্ডার বদলে বিশাল পার্থক্য

void matmul_ikj(float *A, float *B, float *C, int n) {
    memset(C, 0, sizeof(float) * (size_t)n * n);
    for (int i = 0; i \< n; i++)
        for (int k = 0; k \< n; k++) {
            float a = A[i*n + k];                 // একবার লোড, পুরো j-loop জুড়ে reuse
            for (int j = 0; j \< n; j++)
                C[i*n + j] += a * B[k*n + j];
        }
}

এখানে একই তিনটা লুপ, শুধু ক্রম i,k,j। এখন innermost j-loop-এ দুইটাই — B[k*n+j] আর C[i*n+j] — j বাড়ার সাথে sequential (unit-stride)। একটা cache line load হলে তার পুরো ১৬টা float-ই (৬৪B/4B) পরপর ব্যবহার হয়ে যায়, কোনো অপচয় নেই। A[i*n+k] একটা রেজিস্টারে বসে পুরো j-loop জুড়ে reuse হয়, নতুন করে লোডই লাগে না।

গাণিতিকভাবে এটা ঠিক একই যোগফল — শুধু কোন ক্রমে যোগ-গুণ হচ্ছে তার পরিবর্তন (floating-point-এ যোগের ক্রম বদলালে rounding-এ নগণ্য পার্থক্য আসতে পারে, কিন্তু ফলাফল কার্যত অভিন্ন)। results_match()-এ একটা epsilon-tolerance দিয়ে এই দুইটা ভার্সন মিলিয়ে দেখুন:

int results_match(float *ref, float *test, int n, float eps) {
    for (int i = 0; i \< n * n; i++)
        if (fabsf(ref[i] - test[i]) > eps * fmaxf(1.0f, fabsf(ref[i])))
            return 0;
    return 1;
}

৩. Tiled/blocked — cache hierarchy প্রজেক্টের সংখ্যা কাজে লাগানো

ikj ইতিমধ্যে ভালো, কিন্তু বড় n-এ একটা B-এর সম্পূর্ণ row (n floats) আর C-এর row বারবার cache-এ থাকার জন্য প্রতিযোগিতা করে — যদি n এত বড় হয় যে একটা row-ও পুরো L1/L2-তে না আঁটে, reuse আবার নষ্ট হয়। সমাধান: matrix-কে ছোট tile-এ ভাগ করুন, এমন সাইজে যাতে একটা tile-এর A/B/C অংশ cache-এ ধরে যায় — Cache Effect Measurement প্রজেক্টে মাপা আপনার নিজের L1/L2 সাইজ এখানে সরাসরি কাজে লাগবে (৬৪KB L1-এ ~৩২×৩২ float tile আরামসে আঁটে)।

#define T 32   // tile size — নিজের মাপা L1 সাইজ অনুযায়ী টিউন করুন

void matmul_tiled(float *A, float *B, float *C, int n) {
    memset(C, 0, sizeof(float) * (size_t)n * n);
    for (int ii = 0; ii \< n; ii += T)
        for (int kk = 0; kk \< n; kk += T)
            for (int jj = 0; jj \< n; jj += T) {
                int i_max = ii + T \< n ? ii + T : n;
                int k_max = kk + T \< n ? kk + T : n;
                int j_max = jj + T \< n ? jj + T : n;
                for (int i = ii; i \< i_max; i++)
                    for (int k = kk; k \< k_max; k++) {
                        float a = A[i*n + k];
                        for (int j = jj; j \< j_max; j++)
                            C[i*n + j] += a * B[k*n + j];
                    }
            }
}

ভেতরের i,k,j মাইক্রো-কার্নেলটা ঠিক ধাপ ২-এর ikj-ই, শুধু এখন ii/kk/jj-এর তিনটা বাইরের লুপ পুরো কাজটাকে ছোট T×T ব্লকে ভাগ করে দেয়। একটা (ii,kk) ব্লকের জন্য A-এর T×T অংশ, আর একটা (kk,jj) ব্লকের জন্য B-এর T×T অংশ — দুটোই cache-এ থাকা অবস্থায় পুরো ব্যবহার হয়ে যায়, তারপর পরের tile-এ সরে যায়। এটাই cache hierarchy লেসনের “capacity” আর “reuse before eviction” ধারণার সরাসরি বাস্তবায়ন।

৪. SIMD — AVX2 FMA দিয়ে innermost loop

Tiling দিয়ে cache-friendly হলো; এবার SIMD দিয়ে একই কার্নেলের প্রতিটা ইনস্ট্রাকশন ৮টা float একসাথে প্রসেস করাব।

#include <immintrin.h>

void matmul_simd_tiled(float *A, float *B, float *C, int n) {
    memset(C, 0, sizeof(float) * (size_t)n * n);
    for (int ii = 0; ii \< n; ii += T)
      for (int kk = 0; kk \< n; kk += T)
        for (int jj = 0; jj \< n; jj += T) {
          int i_max = ii + T \< n ? ii + T : n;
          int k_max = kk + T \< n ? kk + T : n;
          int j_max = jj + T \< n ? jj + T : n;
          for (int i = ii; i \< i_max; i++)
            for (int k = kk; k \< k_max; k++) {
              __m256 a = _mm256_set1_ps(A[i*n + k]);   // A[i][k] আটটা লেনেই broadcast
              int j = jj;
              for (; j + 8 <= j_max; j += 8) {
                  __m256 b = _mm256_loadu_ps(&B[k*n + j]);
                  __m256 c = _mm256_loadu_ps(&C[i*n + j]);
                  c = _mm256_fmadd_ps(a, b, c);         // c += a*b, ৮টা float একসাথে
                  _mm256_storeu_ps(&C[i*n + j], c);
              }
              for (; j \< j_max; j++)                    // tile সীমানায় বাকি অংশ স্কেলারে
                  C[i*n + j] += A[i*n + k] * B[k*n + j];
          }
        }
}

_mm256_fmadd_ps একটা fused multiply-add — a*b + c একটা মাত্র ইনস্ট্রাকশনে, শুধু দ্রুতই না, একটা মধ্যবর্তী rounding ধাপও বাদ দেয় (তাই মাঝেমধ্যে naive ভার্সনের চেয়ে সামান্য বেশি নিখুঁত!)। j_max % 8 != 0 হলে শেষ কয়েকটা এলিমেন্ট স্কেলার লুপে পড়ে যায় — এই remainder-হ্যান্ডলিং বাদ দিলে buffer-এর বাইরে read/write হয়ে যাবে, তাই এটা optional না, বাধ্যতামূলক। কম্পাইল করুন gcc -O2 -mavx2 -mfma -o matmul matmul.c।

৫. Multi-threaded — row-parallel map pattern

Matrix multiply-এর row-গুলো একে অপরের থেকে সম্পূর্ণ স্বাধীন — C-এর row i গণনা করতে শুধু A-এর row i আর সম্পূর্ণ B লাগে, অন্য কোনো row C[j]-এর উপর নির্ভরতা নেই। এটাই parallel patterns লেসনের map pattern-এর নিখুঁত উদাহরণ — একই ফাংশন স্বাধীনভাবে প্রতিটা row-এ প্রয়োগ করা, কোনো coordination ছাড়াই।

#include <pthread.h>

typedef struct {
    float *A, *B, *C;
    int n, row_start, row_end;
} thread_arg_t;

static void *worker(void *argp) {
    thread_arg_t *arg = (thread_arg_t *)argp;
    int n = arg->n;
    for (int i = arg->row_start; i \< arg->row_end; i++) {
        for (int j = 0; j \< n; j++) arg->C[i*n + j] = 0.0f;
        for (int k = 0; k \< n; k++) {
            __m256 a = _mm256_set1_ps(arg->A[i*n + k]);
            int j = 0;
            for (; j + 8 <= n; j += 8) {
                __m256 b = _mm256_loadu_ps(&arg->B[k*n + j]);
                __m256 c = _mm256_loadu_ps(&arg->C[i*n + j]);
                c = _mm256_fmadd_ps(a, b, c);
                _mm256_storeu_ps(&arg->C[i*n + j], c);
            }
            for (; j \< n; j++)
                arg->C[i*n + j] += arg->A[i*n + k] * arg->B[k*n + j];
        }
    }
    return NULL;
}

void matmul_threaded(float *A, float *B, float *C, int n, int num_threads) {
    pthread_t threads[num_threads];
    thread_arg_t args[num_threads];
    int rows_per_thread = (n + num_threads - 1) / num_threads;

    for (int t = 0; t \< num_threads; t++) {
        int start = t * rows_per_thread;
        int end = start + rows_per_thread \< n ? start + rows_per_thread : n;
        args[t] = (thread_arg_t){ .A = A, .B = B, .C = C, .n = n,
                                   .row_start = start, .row_end = end };
        pthread_create(&threads[t], NULL, worker, &args[t]);
    }
    for (int t = 0; t \< num_threads; t++)
        pthread_join(threads[t], NULL);
}

কেন এটা data-race-মুক্ত: প্রতিটা থ্রেড C-এর সম্পূর্ণ আলাদা row-রেঞ্জে লেখে, কোনো দুইটা থ্রেড কখনো একই মেমরি লোকেশনে লেখে না (A আর B শুধু read হয়, কখনো write না)। তাই কোনো lock বা atomic লাগে না — এটাই embarrassingly parallel সমস্যার সংজ্ঞা। gcc -O2 -mavx2 -mfma -pthread -o matmul matmul.c দিয়ে কম্পাইল করুন। থ্রেড সংখ্যা sysconf(_SC_NPROCESSORS_ONLN) দিয়ে নিজের কোর-সংখ্যা অনুযায়ী ঠিক করুন।

নিজের CPU-র theoretical peak হিসাব করা

Roofline model লেসনের সূত্র সরাসরি এখানে কাজে লাগবে:

Peak GFLOPS = cores × clock_GHz × FLOPs_per_cycle_per_core

AVX2 + একটা FMA পোর্টে: ৮ float/vector × ২ (FMA = mul+add) = ১৬ FLOPs/সাইকেল। দুইটা FMA পোর্ট থাকা কোরে (অনেক desktop CPU-তে সাধারণ): ৩২ FLOPs/সাইকেল। উদাহরণ: ৮ কোর, ৩ GHz, ৩২ FLOPs/সাইকেল → 8 × 3 × 32 = 768 GFLOPS single-precision তাত্ত্বিক সর্বোচ্চ (সব কোর মিলিয়ে)। lscpu/cpuid-এ avx2, fma ফ্ল্যাগ আছে কি না আর FMA পোর্ট সংখ্যা (vendor-এর optimization manual) দেখে নিজের সংখ্যা বসান।

ফলাফল টেবিল — একটা টেমপ্লেট (N=1024, ৮-কোর মেশিন)

সংস্করণGFLOPS (উদাহরণ)Peak-এর %naive-এর তুলনায় speedup
Naive (ijk)~0.6~0.1%1×
ikj reorder~2.5~0.3%~4×
Tiled~5.0~0.7%~8×
Tiled + SIMD (AVX2/FMA)~20.0~2.6%~33×
Tiled + SIMD + threaded (৮ কোর)~130.0~17%~215×

এই সংখ্যাগুলো উদাহরণ — হার্ডওয়্যার-ভেদে বাস্তব মান অনেক আলাদা হবে। যেটা সত্য থাকবে প্রায় সব হার্ডওয়্যারে: naive ঘন ঘন theoretical peak-এর মাত্র ১-৫% পায় (memory-bound, cache miss-এ আটকে), আর পূর্ণ-অপ্টিমাইজড সংস্করণ naive-এর তুলনায় সাধারণত ১০-৫০×+ (কখনো তারও বেশি) দ্রুত হয়, যদিও কোনো সংস্করণই সচরাচর ১০০% peak ছুঁতে পারে না — মেমরি bandwidth, loop overhead, আর tile-boundary অদক্ষতা কিছু না কিছু রেখে দেয়। এটাই roofline model-এর মূল কথা: প্রতিটা algorithm-এর একটা বাস্তবসম্মত সিলিং আছে, আর সেই সিলিং-এর কতটা কাছে পৌঁছানো গেল সেটাই optimization-এর আসল মাপকাঠি, শুধু “দ্রুততর হলো কি না” না।

নিজেকে চ্যালেঞ্জ করুন

  1. Tile size sweep — T = 8, 16, 32, 64, 128 দিয়ে চালিয়ে GFLOPS বনাম tile-size গ্রাফ বানান; কোন সাইজে সবচেয়ে ভালো ফল পাচ্ছেন, আর সেটা আপনার মাপা L1 সাইজের সাথে কতটা সামঞ্জস্যপূর্ণ?
  2. Register blocking — innermost কার্নেলে একবারে ২×২ বা ৪×৪ আউটপুট ব্লক গণনা করুন (একই লোড করা A/B মান কয়েকবার reuse করে memory traffic আরও কমানো)
  3. Thread-count sweep — ১, ২, ৪, ৮, ১৬ থ্রেডে চালিয়ে Amdahl/Gustafson লেসনের speedup curve নিজের ডেটায় প্লট করুন — কোথায় diminishing returns শুরু হয়?
  4. Transpose trick — naive-এর B[k*n+j] কলাম-অ্যাক্সেস সমস্যাটা B-কে আগে থেকে transpose করে (Bt[j*n+k]) সমাধান করা যায় কি না পরীক্ষা করুন — transpose করার খরচ কি লাভের চেয়ে কম?
  5. AVX-512 বা NEON — যদি হার্ডওয়্যার সাপোর্ট করে, একই কার্নেল AVX-512 (১৬ float/vector) বা ARM NEON (৪ float/vector)-এ পোর্ট করে তুলনা করুন

এটা যেখানে গিয়ে মিশবে

এখানে যা শিখলেনপরে কোথায় লাগবে
ikj loop-reorder-এর cache-access যুক্তিএই মডিউলেরই Cache Effect Measurement প্রজেক্ট থেকে পাওয়া নিজের L1/L2 সংখ্যা এখানে সরাসরি tile-size বাছাইয়ে কাজে লাগে
SIMD vectorization (AVX2/FMA)Level 11-এর GPU Architecture ও SIMT লেসন — একই “একসাথে অনেক লেনে গণনা” ধারণা, কিন্তু হাজার হাজার লেনে
Row-parallel map patternLevel 11-এর Amdahl ও Gustafson লেসন — thread-count বাড়ালে speedup কোথায় থেমে যায় তার তত্ত্ব
GFLOPS বনাম theoretical peak মাপাLevel 11-এর Profiling (perf, flamegraph, PMU counters) লেসন — কোথায় সময় যাচ্ছে তা hardware counter দিয়ে সরাসরি দেখা
Compute-bound vs memory-bound classificationLevel 12 (Cloud)-এ CPU-bound workload-এর জন্য container CPU limit/pinning ঠিক করা — ভুল classification মানেই ভুল resource allocation