Jump to content

BLAS: Level 1, 2, 3

From Haidhi Wiki

BLAS (Basic Linear Algebra Subprograms) adalah kosakata standar operasi aljabar linear padat (dense): perkalian vektor, matriks-vektor, matriks-matriks, dan kawan-kawannya. Hampir semua yang ada di atasnya — LAPACK, NumPy, PyTorch, HPL — pada akhirnya memanggil BLAS. Pembagian Level 1, 2, 3 awalnya cuma urutan sejarah, tapi kebetulan (dan ini bagian pentingnya) membelah dunia berdasarkan arithmetic intensity: Level 1–2 selamanya terkurung bandwidth memori, Level 3 satu-satunya yang bisa menyentuh puncak FLOPS mesin.

Prasyarat enak: GFLOPS dan Roofline model.

Tiga Level

Tahun Operasi Contoh Flops Data Intensitas Terikat oleh
Level 1 1979 vektor–vektor axpy, dot O(n) O(n) ~0.08 flop/byte bandwidth
Level 2 1988 matriks–vektor gemv, trsv O(n²) O(n²) ~0.25 flop/byte bandwidth
Level 3 1990 matriks–matriks gemm, trsm O(n³) O(n²) ~n/16 flop/byte compute

Logika intensitasnya:

  • Level 1 — daxpy (y ← αx + y): tiap elemen dibaca sekali, dipakai untuk 2 flop, selesai. 2n flop untuk ~24n byte lalu lintas → 0.083 flop/byte. Mau prosesor secepat apa pun, kecepatannya = kecepatan RAM.
  • Level 2 — dgemv (y ← αAx + βy): matriks A (8n² byte) dibaca sekali untuk 2n² flop → 0.25 flop/byte. Masih nasib yang sama.
  • Level 3 — dgemm (C ← αAB + βC): 2n³ flop atas data hanya O(n²) — tiap elemen A dan B dipakai n kali. Dengan blocking yang benar, data bisa parkir di cache dan dipakai berulang; intensitas naik sebanding n. Inilah satu-satunya operasi BLAS yang bisa lari di atap compute roofline.

Satu kalimat yang layak dihafal: "Level 3 bukan sekadar lebih banyak kerja — dia kerja yang bisa di-cache."

Membaca Nama Rutin

Nama BLAS kelihatan seperti kucing jalan di keyboard (dgemm, ztrsv) padahal sistematis: presisi + jenis matriks + operasi.

Bagian Arti
s / d / c / z float / double / complex / double complex
ge / sy / he / tr / gb general / symmetric / hermitian / triangular / general banded
mm / mv / sv / rk matrix-matrix / matrix-vector / solve / rank-k update

Jadi dgemm = double general matrix-matrix multiply, dan ztrsv = double-complex triangular solve dengan vektor. Rutin yang paling sering ketemu:

Level Rutin Rumus
1 daxpy y ← αx + y
1 ddot hasil ← xᵀy
1 dscal, dnrm2 x ← αx; ‖x‖₂
2 dgemv y ← αAx + βy
2 dger A ← αxyᵀ + A (rank-1 update)
2 dtrsv selesaikan Tx = b, T segitiga
3 dgemm C ← αAB + βC (op(A), op(B) bisa transpose)
3 dtrsm selesaikan TX = B untuk banyak kolom sekaligus
3 dsyrk C ← αAAᵀ + βC

Warisan Fortran: Column-Major & Leading Dimension

BLAS lahir dari Fortran, jadi API aslinya berasumsi matriks tersimpan column-major (kolom demi kolom di memori). C/C++ terbiasa row-major. Interface C-nya, CBLAS (<cblas.h>), memberi parameter pertama CblasRowMajor/CblasColMajor supaya tidak perlu mentranspose kepala sendiri.

Parameter yang paling sering bikin bingung: leading dimension (lda, ldb, ldc) = jarak (dalam elemen) antara awal baris berturut-turut (row-major) atau kolom berturut-turut (col-major) di memori. Untuk matriks utuh, lda = lebar alokasinya. Gunanya: mengoperasikan submatriks tanpa menyalin — 

// A_besar berukuran M x N (row-major). Kita mau blok m x k yang
// pojok kirinya di (bi, bj): cukup geser pointer + lda tetap N.
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
            m, n_blok, k,
            1.0, A_besar + bi * N + bj, N,   // lda = N, BUKAN k!
                 B, ldb, 0.0, C, ldc);

Salah mengisi lda tidak selalu crash — sering cuma hasil salah diam-diam. Kalau hasil gemm terlihat seperti sampah terstruktur, cek urutan major dan leading dimension dulu.

Mengukur Ketiga Level

Program lengkap: satu operasi wakil tiap level, ukur GFLOPS-nya. Compile: 

g++ -O2 -fopenmp -o blas_level blas_level.cpp -lcblas

#include <cstdio>
#include <cblas.h>
#include <omp.h>

// Ukur GFLOPS satu operasi wakil tiap level BLAS.
// Level 1: daxpy  (2n flops, ~24n byte)    -> intensitas ~0.08 flop/byte
// Level 2: dgemv  (2n^2 flops, ~8n^2 byte) -> intensitas ~0.25 flop/byte
// Level 3: dgemm  (2n^3 flops, ~32n^2 byte)-> intensitas ~n/16 flop/byte

int main() {
    // ---- Level 1: y = a*x + y ----
    {
        const long n = 1L << 26;                    // 2 array x 512 MiB
        double *x = new double[n], *y = new double[n];
        #pragma omp parallel for schedule(static)
        for (long i = 0; i < n; ++i) { x[i] = 1.0; y[i] = 2.0; }

        double best = 1e99;
        for (int u = 0; u < 5; ++u) {
            double t0 = omp_get_wtime();
            cblas_daxpy(n, 3.0, x, 1, y, 1);
            double t1 = omp_get_wtime();
            if (t1 - t0 < best) best = t1 - t0;
        }
        printf("L1 daxpy  n=%-9ld : %7.2f GFLOPS  (%.1f GB/s)\n",
               n, 2.0 * n / best / 1e9, 24.0 * n / best / 1e9);
        delete[] x; delete[] y;
    }

    // ---- Level 2: y = A*x ----
    {
        const long n = 8192;                        // A = 512 MiB
        double *A = new double[n * n], *x = new double[n], *y = new double[n];
        #pragma omp parallel for schedule(static)
        for (long i = 0; i < n * n; ++i) A[i] = 1e-3;
        for (long i = 0; i < n; ++i) { x[i] = 1.0; y[i] = 0.0; }

        double best = 1e99;
        for (int u = 0; u < 5; ++u) {
            double t0 = omp_get_wtime();
            cblas_dgemv(CblasRowMajor, CblasNoTrans, n, n,
                        1.0, A, n, x, 1, 0.0, y, 1);
            double t1 = omp_get_wtime();
            if (t1 - t0 < best) best = t1 - t0;
        }
        printf("L2 dgemv  n=%-9ld : %7.2f GFLOPS  (%.1f GB/s)\n",
               n, 2.0 * n * n / best / 1e9, 8.0 * n * n / best / 1e9);
        delete[] A; delete[] x; delete[] y;
    }

    // ---- Level 3: C = A*B ----
    {
        const long n = 1024;
        double *A = new double[n * n], *B = new double[n * n], *C = new double[n * n];
        for (long i = 0; i < n * n; ++i) { A[i] = 1e-3; B[i] = 1e-3; C[i] = 0.0; }

        double best = 1e99;
        for (int u = 0; u < 3; ++u) {
            double t0 = omp_get_wtime();
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n,
                        1.0, A, n, B, n, 0.0, C, n);
            double t1 = omp_get_wtime();
            if (t1 - t0 < best) best = t1 - t0;
        }
        printf("L3 dgemm  n=%-9ld : %7.2f GFLOPS\n", n, 2.0 * n * n * n / best / 1e9);
        delete[] A; delete[] B; delete[] C;
    }
    return 0;
}

Hasil nyata (i5-13450HX, paket Arch cblas = reference BLAS dari Netlib): 

L1 daxpy  n=67108864  :    2.47 GFLOPS  (29.6 GB/s)
L2 dgemv  n=8192      :    2.93 GFLOPS  (11.7 GB/s)
L3 dgemm  n=1024      :    3.77 GFLOPS

Bacaannya:

  • daxpy: 29.6 GB/s ≈ bandwidth memori mesin ini (bandingkan triad ~25 GB/s di OpenMP: Affinity & First-Touch). Level 1 sudah mentok di atap bandwidth — tidak ada implementasi yang bisa membuatnya jauh lebih cepat, fisika sudah bicara.
  • dgemv setali tiga uang, terikat laju pembacaan A.
  • dgemm: 3.77 GFLOPS. Nah ini skandal. Level 3 harusnya compute-bound — mesin ini sanggup ratusan GFLOPS — tapi cuma dapat 4-an. Kenapa? Lanjut ke bagian berikut.

BLAS Adalah Interface, Bukan Library

BLAS hanya mendefinisikan apa fungsinya; seberapa cepat sepenuhnya urusan implementasi:

Implementasi Catatan
Reference (Netlib) Fortran lurus, tanpa blocking/SIMD/threading. Tolok ukur kebenaran, bukan kecepatan. Di Arch ini paket blas/cblas.
OpenBLAS Open source, assembly kernel per arsitektur, multithread. Standar de-facto di Linux. (Arch: openblas / blas-openblas.)
Intel MKL Tercepat di CPU Intel; gratis dipakai.
BLIS Framework riset modern, basis AMD AOCL.
cuBLAS NVIDIA GPU — dipanggil dari host, kerjanya di device.

Bukti bahwa implementasi adalah segalanya — lomba dgemm n=1024, empat kandidat: 

// 1. Tiga loop polos, urutan ijk (inner loop melompati B selebar n)
void naive_ijk(const double *A, const double *B, double *C) {
    for (long i = 0; i < n; ++i)
        for (long j = 0; j < n; ++j) {
            double s = 0.0;
            for (long k = 0; k < n; ++k)
                s += A[i * n + k] * B[k * n + j];
            C[i * n + j] = s;
        }
}

// 2. Loop yang sama, urutan ikj (inner loop kontigu -> SIMD & cache jalan)
void reorder_ikj(const double *A, const double *B, double *C) {
    for (long i = 0; i < n; ++i)
        for (long k = 0; k < n; ++k) {
            const double a = A[i * n + k];
            for (long j = 0; j < n; ++j)
                C[i * n + j] += a * B[k * n + j];
        }
}

// 3. ikj + OpenMP di loop terluar
void omp_ikj(const double *A, const double *B, double *C) {
    #pragma omp parallel for schedule(static)
    for (long i = 0; i < n; ++i)
        for (long k = 0; k < n; ++k) {
            const double a = A[i * n + k];
            for (long j = 0; j < n; ++j)
                C[i * n + j] += a * B[k * n + j];
        }
}

// 4. Serahkan ke profesional
void blas_dgemm(const double *A, const double *B, double *C) {
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n,
                1.0, A, n, B, n, 0.0, C, n);
}

Hasil nyata (16 thread): 

n = 1024, threads = 16

naive ijk          :   3.294 s     0.65 GFLOPS
reorder ikj        :   0.436 s     4.92 GFLOPS
ikj + OpenMP       :   0.048 s    45.03 GFLOPS
cblas_dgemm (ref)  :   0.406 s     5.29 GFLOPS

Empat pelajaran dari satu tabel:

  1. Urutan loop = 7.6×. ijk → ikj mengubah akses B dari lompat-lompat selebar baris jadi kontigu: Cache line terpakai penuh dan compiler berani auto-vectorize (SIMD). Belum ada paralelisme sama sekali.
  2. OpenMP menambah 9× lagi di 16 thread.
  3. Kode 20 baris kita mengalahkan BLAS terinstall 8.5×. Bukan karena kita hebat — karena yang terinstall reference BLAS. Ini jebakan nyata: di banyak distro dan cluster, -lblas default menaut ke Netlib. Selalu cek implementasi apa yang sebenarnya ke-link (ldd ./program).
  4. 45 GFLOPS itu tetap kelas amatir. OpenBLAS/MKL di mesin yang sama tipikalnya ratusan GFLOPS (puncak teoretis mesin ini ~400+ GFLOPS DP) lewat blocking bertingkat mengikuti hierarki cache + microkernel assembly + packing. Teknik blocking-nya persis ide yang sama dengan CUDA: Matmul dengan Tiling — beda alamat, ide serupa.

Moral: jangan pernah menulis gemm sendiri untuk produksi. Tulis sekali untuk belajar (seperti di atas), lalu panggil library selamanya.

Kenapa Level 3 Itu Raja

Karena hanya Level 3 yang bisa lari di atap compute, seni menyusun algoritma numerik = memindahkan sebanyak mungkin kerja ke Level 3:

  • LAPACK menyusun ulang faktorisasi (Dekomposisi LU & Cholesky) jadi versi blocked: panel difaktorkan dengan Level 2 (porsi kecil), lalu sisanya di-update besar-besaran via dgemm/dtrsm (porsi raksasa). Rasio kerja Level 3-nya mendekati 100% saat matriks besar.
  • Deep learning duduk di atas gemm: layer fully-connected = gemm, convolution di-lower jadi gemm (im2col), attention = gemm bertumpuk. GPU tensor core pun pada dasarnya mesin gemm.
  • Pola praktis: kalau kode memanggil gemv berulang dengan matriks yang sama, tumpuk vektor-vektornya jadi matriks dan panggil gemm sekali. Kerjanya (flop) sama persis; lalu lintas memorinya runtuh karena A dibaca sekali, bukan k kali.

Catatan batas wilayah: semua ini dunia dense. Begitu matriksnya sparse (mayoritas nol — tipikal PDE, graf), format dan pustakanya beda total: Sparse Matrix: Format CSR & CSC dan SpMV — yang sayangnya kembali bernasib Level-2: memory-bound.

Jebakan Umum

  • Menaut reference BLAS tanpa sadar. Gejala: gemm "cuma" beberapa GFLOPS (terbukti di atas). Diagnosa: ldd ./program | grep blas. Obat di Arch: pasang openblas + link -lopenblas (atau paket blas-openblas yang menggantikan provider libblas.so sekalian).
  • Row-major vs column-major kebalik. Hasilnya transpose diam-diam. Pakai CBLAS dan konsisten CblasRowMajor, atau manfaatkan identitas (AB)ᵀ = BᵀAᵀ.
  • Perang thread. OpenBLAS punya thread pool sendiri. Memanggil BLAS dari dalam region OpenMP = oversubscription (16×16 thread). Aturan: BLAS multithread di luar region paralel, atau set OPENBLAS_NUM_THREADS=1 saat dipanggil dari dalam.
  • Matriks kecil. Overhead pemanggilan + setup mengalahkan kerjanya untuk n ≲ 50-an; kalau punya ribuan matriks kecil, cari batched BLAS (cublasDgemmBatched, MKL batch API), jangan loop.
  • Lupa β. dgemm menghitung C ← αAB + βC. β=0 menimpa C (dan tidak membaca isinya), β=1 mengakumulasi. Salah pilih = bug hasil yang halus.
  • Increment (incx) eksotis. daxpy dkk. menerima stride antar elemen. Berguna untuk kolom di matriks row-major, tapi stride besar = tiap elemen satu cache line = bandwidth efektif ambruk 8×.

Ringkasan

  • Level = kelas intensitas aritmetika: L1 (vektor, 0.08 flop/byte) dan L2 (matriks-vektor, 0.25) terkurung bandwidth; hanya L3 (matriks-matriks, ~n/16) yang bisa menyentuh puncak FLOPS.
  • Nama rutin = presisi + jenis matriks + operasi: dgemm = double, general, matrix-matrix.
  • BLAS itu interface. Reference Netlib untuk kebenaran; OpenBLAS/MKL/cuBLAS untuk kecepatan. Cek ldd, jangan berasumsi.
  • Terbukti di mesin sendiri: loop order 7.6×, +OpenMP 9×, dan reference BLAS kalah dari 20 baris kode — tapi implementasi serius masih ~8× di atasnya lagi.
  • Strategi algoritma numerik modern: dorong semua kerja ke Level 3.

Lihat Juga

Kategori:HPC Kategori:Metode Numerik

Retrieved from "https://wiki.haidhi.my.id/wiki/index.php?title=BLAS:_Level_1,_2,_3&oldid=798"