View source for BLAS: Level 1, 2, 3

'''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, [[Benchmark: HPL & HPCG|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 =

{| class="wikitable"
! !! Tahun !! Operasi !! Contoh !! Flops !! Data !! Intensitas !! Terikat oleh
|-
| '''Level 1''' || 1979 || vektor–vektor || <code>axpy</code>, <code>dot</code> || O(n) || O(n) || ~0.08 flop/byte || bandwidth
|-
| '''Level 2''' || 1988 || matriks–vektor || <code>gemv</code>, <code>trsv</code> || O(n²) || O(n²) || ~0.25 flop/byte || bandwidth
|-
| '''Level 3''' || 1990 || matriks–matriks || <code>gemm</code>, <code>trsm</code> || O(n³) || O(n²) || ~n/16 flop/byte || '''compute'''
|}

Logika intensitasnya:
* '''Level 1''' — <code>daxpy</code> (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''' — <code>dgemv</code> (y ← αAx + βy): matriks A (8n² byte) dibaca sekali untuk 2n² flop → 0.25 flop/byte. Masih nasib yang sama.
* '''Level 3''' — <code>dgemm</code> (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 model|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 (<code>dgemm</code>, <code>ztrsv</code>) padahal sistematis: '''presisi + jenis matriks + operasi'''.

{| class="wikitable"
! Bagian !! Arti
|-
| <code>s</code> / <code>d</code> / <code>c</code> / <code>z</code> || float / double / complex / double complex
|-
| <code>ge</code> / <code>sy</code> / <code>he</code> / <code>tr</code> / <code>gb</code> || general / symmetric / hermitian / triangular / general banded
|-
| <code>mm</code> / <code>mv</code> / <code>sv</code> / <code>rk</code> || matrix-matrix / matrix-vector / solve / rank-k update
|}

Jadi <code>dgemm</code> = '''d'''ouble '''ge'''neral '''m'''atrix-'''m'''atrix multiply, dan <code>ztrsv</code> = double-complex triangular solve dengan vektor. Rutin yang paling sering ketemu:

{| class="wikitable"
! Level !! Rutin !! Rumus
|-
| 1 || <code>daxpy</code> || y ← αx + y
|-
| 1 || <code>ddot</code> || hasil ← xᵀy
|-
| 1 || <code>dscal</code>, <code>dnrm2</code> || x ← αx; ‖x‖₂
|-
| 2 || <code>dgemv</code> || y ← αAx + βy
|-
| 2 || <code>dger</code> || A ← αxyᵀ + A (rank-1 update)
|-
| 2 || <code>dtrsv</code> || selesaikan Tx = b, T segitiga
|-
| 3 || <code>dgemm</code> || C ← αAB + βC (op(A), op(B) bisa transpose)
|-
| 3 || <code>dtrsm</code> || selesaikan TX = B untuk banyak kolom sekaligus
|-
| 3 || <code>dsyrk</code> || 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''' (<code><cblas.h></code>), memberi parameter pertama <code>CblasRowMajor</code>/<code>CblasColMajor</code> supaya tidak perlu mentranspose kepala sendiri.

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

<syntaxhighlight lang="cpp">
// 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);
</syntaxhighlight>

Salah mengisi <code>lda</code> tidak selalu crash — sering cuma hasil salah diam-diam. Kalau hasil <code>gemm</code> 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:

<pre>
g++ -O2 -fopenmp -o blas_level blas_level.cpp -lcblas
</pre>

<syntaxhighlight lang="cpp">
#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;
}
</syntaxhighlight>

Hasil nyata (i5-13450HX, paket Arch <code>cblas</code> = '''reference BLAS dari Netlib'''):

<pre>
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
</pre>

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:

{| class="wikitable"
! Implementasi !! Catatan
|-
| '''Reference (Netlib)''' || Fortran lurus, tanpa blocking/SIMD/threading. Tolok ukur '''kebenaran''', bukan kecepatan. Di Arch ini paket <code>blas</code>/<code>cblas</code>.
|-
| '''OpenBLAS''' || Open source, assembly kernel per arsitektur, multithread. Standar de-facto di Linux. (Arch: <code>openblas</code> / <code>blas-openblas</code>.)
|-
| '''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 <code>dgemm</code> n=1024, empat kandidat:

<syntaxhighlight lang="cpp">
// 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);
}
</syntaxhighlight>

Hasil nyata (16 thread):

<pre>
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
</pre>

Empat pelajaran dari satu tabel:
# '''Urutan loop = 7.6×.''' <code>ijk</code> → <code>ikj</code> mengubah akses B dari lompat-lompat selebar baris jadi kontigu: [[Cache line]] terpakai penuh dan compiler berani auto-vectorize ([[SIMD: Auto-Vectorization & Laporan Compiler|SIMD]]). Belum ada paralelisme sama sekali.
# '''OpenMP menambah 9×''' lagi di 16 thread.
# '''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, <code>-lblas</code> default menaut ke Netlib. '''Selalu cek implementasi apa yang sebenarnya ke-link''' (<code>ldd ./program</code>).
# '''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 <code>gemm</code> 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 <code>dgemm</code>/<code>dtrsm</code> (porsi raksasa). Rasio kerja Level 3-nya mendekati 100% saat matriks besar.
* Deep learning duduk di atas <code>gemm</code>: 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 <code>gemv</code> berulang dengan matriks yang sama, tumpuk vektor-vektornya jadi matriks dan panggil <code>gemm</code> 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: <code>ldd ./program | grep blas</code>. Obat di Arch: pasang <code>openblas</code> + link <code>-lopenblas</code> (atau paket <code>blas-openblas</code> yang menggantikan provider <code>libblas.so</code> sekalian).
* '''Row-major vs column-major kebalik.''' Hasilnya transpose diam-diam. Pakai CBLAS dan konsisten <code>CblasRowMajor</code>, atau manfaatkan identitas (AB)ᵀ = BᵀAᵀ.
* '''Perang thread.''' OpenBLAS punya thread pool sendiri. Memanggil BLAS dari dalam region [[OpenMP: Parallel For & Scheduling|OpenMP]] = oversubscription (16×16 thread). Aturan: BLAS multithread ''di luar'' region paralel, atau set <code>OPENBLAS_NUM_THREADS=1</code> saat dipanggil dari dalam.
* '''Matriks kecil.''' Overhead pemanggilan + setup mengalahkan kerjanya untuk n ≲ 50-an; kalau punya ribuan matriks kecil, cari ''batched BLAS'' (<code>cublasDgemmBatched</code>, MKL batch API), jangan loop.
* '''Lupa β.''' <code>dgemm</code> menghitung C ← αAB + '''β'''C. β=0 menimpa C (dan tidak membaca isinya), β=1 mengakumulasi. Salah pilih = bug hasil yang halus.
* '''Increment (incx) eksotis.''' <code>daxpy</code> 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: <code>dgemm</code> = double, general, matrix-matrix.
* BLAS itu interface. Reference Netlib untuk kebenaran; OpenBLAS/MKL/cuBLAS untuk kecepatan. Cek <code>ldd</code>, 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 =

* [[GFLOPS]] dan [[Roofline model]] — kerangka membaca angka-angka di atas
* [[Cache line]] dan [[SIMD]] — kenapa urutan loop berdampak 7.6×
* [[Dekomposisi LU & Cholesky]] — konsumen terbesar Level 3
* [[Sparse Matrix: Format CSR & CSC]], [[SpMV]] — dunia sebelah, saat matriks kebanyakan nol
* [[CUDA: Matmul dengan Tiling]] — blocking yang sama, di GPU
* [[Benchmark: HPL & HPCG]] — HPL pada dasarnya kontes dgemm sedunia
* [[OpenMP: Parallel For & Scheduling]] — dipakai di versi 45 GFLOPS di atas

[[Kategori:HPC]]
[[Kategori:Metode Numerik]]