[Fortranのはじめ方] Part23: 行列計算とBLAS/LAPACKの連携

Fortranで高速な数値計算を実現するための標準ライブラリ活用法

科学技術計算において、行列計算は非常に重要な役割を果たします。例えば、連立一次方程式の求解、固有値問題、最小二乗問題など、多くの問題が行列演算に帰着します。Fortranはその成り立ちから数値計算を得意としており、これらの計算を効率的に行うための強力なライブラリが整備されています。

今回は、その中でも特に標準的で広く使われている線形計算ライブラリである BLAS (Basic Linear Algebra Subprograms) と LAPACK (Linear Algebra PACKage) について学び、Fortranプログラムからこれらを連携させて利用する方法を見ていきましょう。

BLASとLAPACKとは？ 🤔

BLAS (Basic Linear Algebra Subprograms)

BLASは、ベクトルや行列に関する基本的な演算（加算、スカラー倍、内積、行列ベクトル積、行列行列積など）を行うためのサブルーチン群の標準仕様です。
BLASは、演算の種類と対象によって3つのレベルに分けられています。

Level 1: ベクトル-ベクトル演算 (例: ddot – 倍精度実数ベクトルの内積, daxpy – 定数倍加算)
Level 2: 行列-ベクトル演算 (例: dgemv – 倍精度実数行列とベクトルの積)
Level 3: 行列-行列演算 (例: dgemm – 倍精度実数行列と行列の積)

特にLevel 3の演算は計算量が多いため、キャッシュメモリなどを効率的に利用した高速な実装が重要になります。

LAPACK (Linear Algebra PACKage)

LAPACKは、BLASを基礎として構築された、より高度な線形計算（連立一次方程式、線形最小二乗問題、固有値問題、特異値問題など）を行うためのサブルーチン群です。
BLASが基本的な演算の「部品」を提供するのに対し、LAPACKはそれらの部品を組み合わせて、より複雑な問題を解くための「完成品」に近い機能を提供します。

例: dgesv – 倍精度実数一般行列を用いた連立一次方程式 Ax=b の解法
例: dgeev – 倍精度実数一般行列の固有値・固有ベクトルの計算

LAPACKのルーチンは内部で効率的にBLASルーチンを呼び出すように設計されています。

なぜBLAS/LAPACKを使うのか？ ✨

自分で一から行列計算のコードを書くことも可能ですが、BLAS/LAPACKを利用することには大きなメリットがあります。

パフォーマンス: BLAS/LAPACKの実装（特にベンダー提供のものやOpenBLASなど）は、特定のCPUアーキテクチャ向けに高度に最適化されています。これにより、自分で書いたコードよりも大幅に高速な計算が期待できます。キャッシュ効率、ベクトル化、並列化などが考慮されています。
信頼性: 長年にわたり世界中の研究者や技術者によって利用され、テストされてきた実績があり、バグが少なく信頼性の高い計算結果を得られます。
標準化と移植性: インターフェース（サブルーチン名や引数）が標準化されているため、異なる環境（OS、コンパイラ、ハードウェア）でも同じように利用でき、コードの移植性が高まります。
開発効率: 実績のある高度なアルゴリズムを自分で実装する手間が省け、本来の目的である科学技術計算そのものに集中できます。

FortranからBLAS/LAPACKを使う方法 💻

FortranプログラムからBLAS/LAPACKのサブルーチンを利用するのは比較的簡単です。主に以下のステップで行います。

BLAS/LAPACKライブラリの準備: システムにBLAS/LAPACKライブラリがインストールされている必要があります。多くの場合、OSのパッケージマネージャ（Linuxのaptやyum、macOSのHomebrewなど）でインストールできます。高性能な実装（OpenBLAS, Intel MKLなど）を利用することも可能です。
プログラムでの呼び出し: Fortranコード内で、使用したいBLAS/LAPACKサブルーチンを CALL 文で呼び出します。多くの場合、引数は参照渡し（変数をそのまま渡す）です。サブルーチン名や引数の仕様は、BLAS/LAPACKのリファレンスで確認します。
コンパイルとリンク: プログラムをコンパイルする際に、使用するBLAS/LAPACKライブラリをリンクするようコンパイラに指示します。

簡単な例: ベクトルの内積 (BLAS Level 1: DDOT)

2つの倍精度実数ベクトル `x` と `y` の内積を計算する例です。


program dot_product_example
  implicit none
  integer, parameter :: n = 3
  double precision :: x(n), y(n), result
  external :: ddot ! BLASのddot関数を外部関数として宣言

  ! ベクトルに値を設定
  x = [1.0d0, 2.0d0, 3.0d0]
  y = [4.0d0, 5.0d0, 6.0d0]

  ! ddotを呼び出して内積を計算
  ! 引数: N (要素数), X (ベクトル1), INCX (Xの要素間隔), Y (ベクトル2), INCY (Yの要素間隔)
  result = ddot(n, x, 1, y, 1)

  ! 結果の表示
  write(*,*) 'Vector X:', x
  write(*,*) 'Vector Y:', y
  write(*,*) 'Dot Product (X . Y) =', result ! 1*4 + 2*5 + 3*6 = 4 + 10 + 18 = 32

end program dot_product_example

コンパイルとリンク（gfortranと標準的なBLASライブラリの場合）:


gfortran dot_product_example.f90 -o dot_example -lblas
./dot_example

-lblas オプションでBLASライブラリをリンクしています。（環境によっては -ldot など、または特定のライブラリパスの指定が必要な場合があります。）

例: 連立一次方程式の解法 (LAPACK: DGESV)

連立一次方程式 Ax = b を解く例です。ここで A は正方行列、x と b はベクトルです。dgesv サブルーチンを使います。


program solve_linear_equation
  implicit none
  integer, parameter :: n = 2 ! 行列のサイズ (N x N)
  integer, parameter :: nrhs = 1 ! 右辺ベクトル(b)の数
  double precision :: a(n, n) ! 係数行列 A
  double precision :: b(n, nrhs) ! 右辺ベクトル b (複数列も可能)
  integer :: ipiv(n) ! ピボット選択のための作業用配列
  integer :: info ! 実行結果ステータス (0なら成功)
  external :: dgesv ! LAPACKのdgesvサブルーチン

  ! 行列Aとベクトルbを設定
  ! | 2  1 | | x1 | = | 5 |
  ! | 1  3 | | x2 | = | 5 |
  ! 答えは x1=2, x2=1

  ! Fortranは列優先(column-major)なので注意して設定
  a(:, 1) = [2.0d0, 1.0d0]
  a(:, 2) = [1.0d0, 3.0d0]

  b(:, 1) = [5.0d0, 5.0d0]

  write(*, '(a)') "Solving A * x = b"
  write(*, '(a, 2f6.1)') "Matrix A (col1):", a(:, 1)
  write(*, '(11x, 2f6.1)') "         (col2):", a(:, 2)
  write(*, '(a, 2f6.1)') "Vector b:", b(:, 1)

  ! dgesvを呼び出して方程式を解く
  ! 引数: N (次元), NRHS (右辺数), A (係数行列, 解bで上書きされる), LDA (Aの第一次元),
  !       IPIV (作業配列), B (右辺ベクトル, 解xで上書きされる), LDB (Bの第一次元), INFO (結果)
  call dgesv(n, nrhs, a, n, ipiv, b, n, info)

  ! 結果の確認
  if (info == 0) then
    write(*, '(a, 2f6.1)') "Solution x:", b(:, 1)
  else if (info > 0) then
    write(*, *) "Error: Matrix is singular. INFO =", info
  else
    write(*, *) "Error: Illegal argument. INFO =", info
  end if

end program solve_linear_equation

コンパイルとリンク（gfortranと標準的なLAPACK, BLASライブラリの場合）:


gfortran solve_linear_equation.f90 -o solve_example -llapack -lblas
./solve_example

-llapack -lblas オプションでLAPACKと、それが依存するBLASライブラリをリンクしています。（順番が重要な場合があります。）

注意点:

BLAS/LAPACKのルーチンはFortran 77スタイルで書かれていることが多く、引数の型や渡し方に注意が必要です。
多くのルーチンでは、入力として与えた行列やベクトルが計算結果で上書きされることがあります。必要なら事前にコピーを取っておきましょう。
Fortranの配列は列優先 (Column-major) ですが、C言語など他の言語とのインターフェースでは行優先 (Row-major) が一般的なため、混在環境では注意が必要です。
サブルーチンの引数は非常に多い場合があります。リファレンスマニュアルをよく読んで、各引数の意味（特に作業用配列のサイズなど）を正確に理解することが重要です。

ライブラリの選択と入手

利用可能なBLAS/LAPACKライブラリはいくつかあります。

ライブラリ名	特徴	入手方法例
Netlib LAPACK/BLAS	リファレンス実装。移植性は高いが、最適化レベルは限定的。	公式サイトからソースコードをダウンロードしてビルド。一部OSではパッケージとしても提供。
OpenBLAS	オープンソース。多くのアーキテクチャ向けに最適化されており、高性能。	`sudo apt-get install libopenblas-dev` (Debian/Ubuntu), `brew install openblas` (macOS) など。
Intel Math Kernel Library (MKL)	Intel CPU向けに高度に最適化された商用（現在は無料利用可）ライブラリ。非常に高性能。	Intel oneAPI Base Toolkitの一部としてインストール。
ATLAS (Automatically Tuned Linear Algebra Software)	オープンソース。インストール時にシステムに合わせて自動的に最適化（チューニング）を行う。	`sudo apt-get install libatlas-base-dev` (Debian/Ubuntu) など。
AMD Optimizing CPU Libraries (AOCL)	AMD CPU向けに最適化されたライブラリ。BLAS, LAPACK, FFTなどを含む。	AMDのサイトからダウンロード。

どのライブラリを選択するかは、使用する環境（CPU、OS）や求める性能によって異なります。多くの場合、OpenBLASやIntel MKLが高性能な選択肢となります。コンパイル時のリンクオプションは、使用するライブラリによって異なりますので、各ライブラリのドキュメントを確認してください。

まとめ 💡

BLASとLAPACKは、Fortranで高速かつ信頼性の高い行列計算を行うための強力なツールです。これらの標準ライブラリを利用することで、複雑な線形計算アルゴリズムの実装から解放され、科学技術計算の本質的な部分に集中することができます。

最初は多くのサブルーチンや引数に戸惑うかもしれませんが、基本的な使い方をマスターすれば、その恩恵は計り知れません。特に大規模な問題を扱う際には、最適化されたBLAS/LAPACKライブラリの利用が計算時間短縮の鍵となります 🔑。

ぜひ、ご自身のFortranプログラムにBLAS/LAPACKを導入して、その性能を体験してみてください！

参考情報 📚

Netlib LAPACK: https://www.netlib.org/lapack/ (公式リファレンス、ユーザーガイドなど)
Netlib BLAS: https://www.netlib.org/blas/ (BLASの仕様など)
LAPACK Users’ Guide (Chapter 2): https://www.netlib.org/lapack/lug/node14.html (サブルーチンの命名規則など)
OpenBLAS: https://www.openblas.net/
Intel Math Kernel Library (MKL): https://www.intel.com/content/www/us/en/developer/tools/oneapi/mkl.html