KAL
Menu

BLAS dan LAPACK

Memakai dgemm dan rutin BLAS serta LAPACK dari kode Fortran modern dengan OpenBLAS.

Daftar isi

BLAS (Basic Linear Algebra Subprograms) dan LAPACK (Linear Algebra PACKage) adalah dua pustaka numerik yang menjadi fondasi hampir semua perangkat lunak komputasi ilmiah, dan implementasi referensinya ditulis dalam Fortran 77. OpenBLAS adalah implementasi modern yang menyediakan rutin BLAS sekaligus LAPACK dalam satu pustaka. Topik ini membahas pemanggilan dgemm untuk perkalian matriks dan dgesv untuk penyelesaian sistem persamaan linear dari kode Fortran modern, melalui antarmuka bind(c) yang diarahkan pada simbol pustaka tersebut. Kedua program disertakan, dikompilasi dengan OpenBLAS, lalu dijalankan, dan keluarannya diverifikasi terhadap perhitungan manual.

dgemm dan antarmuka pustaka BLAS#

Rutin dgemm menghitung C = alpha*A*B + beta*C untuk matriks real double precision, dengan opsi transposisi pada A dan B melalui argumen karakternya. Simbol yang diekspor OpenBLAS mengikuti konvensi antarmuka Fortran lama, yaitu nama rutin diberi tambahan garis bawah, sehingga atribut antarmukanya adalah bind(c, name="dgemm_"). Argumen karakter dinyatakan character(kind=c_char) tanpa atribut value, karena pihak C menerimanya sebagai char*. Argumen matriks dinyatakan berbentuk eksplisit seperti a(lda, *), dengan lda sebagai dimensi utama (leading dimension) yang pada kasus terpadat sama dengan jumlah baris matriksnya.

BLAS tersusun bertingkat menurut bentuk operasinya. Level 1 mencakup operasi vektor seperti perkalian titik, level 2 mencakup operasi matriks-vektor, dan level 3 mencakup operasi matriks-matriks seperti dgemm. Nama GEMM merupakan singkatan general matrix multiply, dan rutin inilah mesin utama aljabar linear numerik karena LAPACK sendiri memanggilnya pada tahap pemfaktoran matriks. Biaya dgemm tumbuh dengan orde O(n3)O(n^3) untuk matriks persegi berukuran n, sedangkan jumlah datanya hanya tumbuh dengan orde O(n2)O(n^2), sehingga pemanfaatan cache yang dilakukan OpenBLAS memberi perbesaran performa yang besar pada ukuran yang besar.

Program berikut mengalikan dua matriks 2 x 2. Perhatikan urutan nilai pada reshape, yang mengisi matriks menurut kolom sebagaimana tata letak memori Fortran.

fortran demo-blas.f90 
program demo_blas
   use, intrinsic :: iso_c_binding
   implicit none
   interface
      subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc) &
           bind(c, name="dgemm_")
         import :: c_char, c_int, c_double
         character(kind=c_char), intent(in) :: transa, transb
         integer(c_int), intent(in) :: m, n, k, lda, ldb, ldc
         real(c_double), intent(in) :: alpha, beta
         real(c_double), intent(in) :: a(lda, *)
         real(c_double), intent(in) :: b(ldb, *)
         real(c_double), intent(inout) :: c(ldc, *)
      end subroutine dgemm
   end interface
   real(c_double) :: a(2,2), b(2,2), c(2,2)
   real(c_double) :: alpha, beta
   a = reshape([1.0_c_double, 3.0_c_double, 2.0_c_double, 4.0_c_double], [2,2])
   b = reshape([5.0_c_double, 7.0_c_double, 6.0_c_double, 8.0_c_double], [2,2])
   alpha = 1.0_c_double
   beta = 0.0_c_double
   call dgemm("N", "N", 2_c_int, 2_c_int, 2_c_int, alpha, a, 2_c_int, &
        & b, 2_c_int, beta, c, 2_c_int)
   print *, "Matriks C hasil A x B:"
   print '(2F7.1)', c(1,:), c(2,:)
end program demo_blas
bash 
# jalur pustaka OpenBLAS pada lingkungan jalur ini
OBLIB=/home/kal/my-project/toolchain/fortran/usr/lib/x86_64-linux-gnu/openblas-pthread
# kompilasi dan sambungkan program dengan OpenBLAS
gfortran demo-blas.f90 -o demo-blas -L"$OBLIB" -lopenblas -Wl,-rpath,"$OBLIB"
# jalankan program hasil kompilasi
./demo-blas
text Keluaran 
 Matriks C hasil A x B:
   19.0   22.0
   43.0   50.0

Dengan tata letak kolom utama, deret [1, 3, 2, 4] membentuk matriks A yang baris pertamanya adalah 1 dan 2, sedangkan deret [5, 7, 6, 8] membentuk B yang baris pertamanya adalah 5 dan 6. Hasil kali keduanya adalah C dengan baris pertama 19 dan 22 serta baris kedua 43 dan 50, yang sama persis dengan keluaran. Argumen "N" menyatakan kedua matriks dipakai tanpa transposisi, dan beta bernilai nol sehingga C sebelumnya tidak diikutkan.

dgesv untuk sistem persamaan linear#

Rutin dgesv menyelesaikan sistem persamaan linear berikut.

Ax=b A x = b

Sistem diselesaikan dengan dekomposisi LU yang disertai pemivotan parsial, sehingga dgesv juga mengembalikan matriks faktor pada A dan vektor pivot pada ipiv. Tanda tangannya adalah dgesv(N, NRHS, A, LDA, IPIV, B, LDB, INFO), dengan NRHS sebagai jumlah ruas kanan, LDA dan LDB sebagai dimensi utama A dan B, serta info sebagai kode status yang bernilai nol apabila berhasil. Argumen ldb wajib disertakan meskipun jumlah ruas kanannya satu, karena jumlah argumen rutin Fortran tidak diperiksa kompiler terhadap antarmuka buatan sendiri.

fortran demo-lapack.f90 
program demo_lapack
   use, intrinsic :: iso_c_binding
   implicit none
   interface
      subroutine dgesv(n, nrhs, a, lda, ipiv, b, ldb, info) bind(c, name="dgesv_")
         import :: c_int, c_double
         integer(c_int), intent(in) :: n, nrhs, lda, ldb
         real(c_double), intent(inout) :: a(lda, *)
         real(c_double), intent(inout) :: b(ldb, *)
         integer(c_int), intent(out) :: ipiv(*)
         integer(c_int), intent(out) :: info
      end subroutine dgesv
   end interface
   real(c_double) :: a(2,2), b(2)
   integer(c_int) :: ipiv(2), info
   a = reshape([2.0_c_double, 1.0_c_double, 1.0_c_double, 3.0_c_double], [2,2])
   b = [5.0_c_double, 10.0_c_double]
   call dgesv(2_c_int, 1_c_int, a, 2_c_int, ipiv, b, 2_c_int, info)
   print *, "solusi x:", b
   print *, "info    :", info
end program demo_lapack
bash 
# jalur pustaka OpenBLAS pada lingkungan jalur ini
OBLIB=/home/kal/my-project/toolchain/fortran/usr/lib/x86_64-linux-gnu/openblas-pthread
# kompilasi dan sambungkan program dengan OpenBLAS
gfortran demo-lapack.f90 -o demo-lapack -L"$OBLIB" -lopenblas -Wl,-rpath,"$OBLIB"
# jalankan program hasil kompilasi
./demo-lapack
text Keluaran 
 solusi x:   1.0000000000000000        3.0000000000000000     
 info    :           0

Sistem yang diselesaikan adalah 2x + y = 5 dan x + 3y = 10, yang solusinya x = 1 dan y = 3, sama dengan keluaran. Nilai info yang nol menegaskan tidak ada pembagian dengan nol maupun baris yang bergantung linear, sedangkan nilai positif menyatakan faktor U pada posisi tersebut tepat nol, sehingga matriksnya singular dan proses dihentikan. Vektor ipiv mencatat baris yang ditukar pada setiap langkah pemivotan, dan catatan itu diperlukan kembali apabila faktorisasi yang sama dihitung ulang secara manual. Setelah rutin selesai, A berisi faktor LU, sehingga penyelesaian ruas kanan lain dari matriks yang sama dapat dilanjutkan dengan rutin dgetrs tanpa memfaktorkan ulang.

Kesalahan umum#

  • Antarmuka yang tidak lengkap. Argumen dgesv berjumlah delapan, dan satu argumen yang terlewat seperti ldb tidak terdeteksi kompiler karena antarmukanya buatan pengguna. Program tetap tersambung, tetapi rutin membaca alamat yang salah saat berjalan dan program berhenti dengan segfault. Tanda tangan rutin selalu dicocokkan dengan halaman dokumentasi resminya sebelum antarmuka ditulis.
  • Lupa garis bawah pada nama simbol. Antarmuka dengan name="dgesv" tanpa garis bawah menghasilkan error undefined reference pada tahap penyambungan, karena simbol yang diekspor OpenBLAS adalah dgesv_.
  • Urutan reshape yang diabaikan. Matriks yang diisi menurut baris tanpa sadar membentuk sistem persamaan yang berbeda, dan rutin tetap mengembalikan solusi yang salah tanpa pesan error. Perilaku tata letak ini dibahas rinci pada topik tentang array lintas bahasa di modul ini.

Ringkasan#

  • OpenBLAS menyediakan rutin BLAS dan LAPACK dengan simbol bergaya Fortran lama, sehingga antarmukanya ditulis bind(c, name="dgemm_") beserta garis bawahnya.
  • dgemm mengalikan matriks dengan dimensi utama yang eksplisit, sedangkan dgesv menyelesaikan Ax = b melalui dekomposisi LU dan melaporkan statusnya pada info.
  • Kekurangan satu argumen pada antarmuka buatan sendiri tidak terdeteksi kompiler dan berujung pada kerusakan memori saat program berjalan.

Referensi#