C++ 開発で LAPACK を使ったエルミート行列対角化をする
はじめに
C++ 開発で LAPACK を使う場合のための覚書として記す。 ソースコードとコンパイルオプションを明示して、最低限動かすために必要な情報を残す。 また、REFERENCE BLAS と OpenBLAS の速度比較も行う。
BLAS と LAPACK
BLAS は Basic Linear Algebra Subprograms の略で線形代数に必要な基礎的な演算を高速に行う数値計算ライブラリである。 これには内積やノルムの計算、ベクトルや行列の乗算などが含まれている。
BLAS にはいくつかのバリエーションがあり、Netlib で公開されているものが公式の BLAS となる。 これは特に REFERENCE BLAS とも呼ばれている。 つまり、すべての BLAS の「お手本」である。 ドキュメントが整っており、どのバリエーションを使う場合も Netlib のドキュメントを確認すればよい。 REFERENCE BLAS は高速な計算アルゴリズムを示すためのもので、コードは可読性に重点を置かれている。 これに対して、より速度に重点をおいて種々の最適化を施された BLAS を最適化 BLAS (Optimized BLAS) という。 最適化 BLAS には現在有名なものに MKL (Intel Math Kernel Library)、OpenBLAS などがある。 (より正確に言えば MKL は BLAS 以外にも LAPACK 相当の機能も含んでおり、更に多くの数値計算ルーチンを包含している。)
LAPACK は Linear Algebra PACKage の略でより一般的な線形代数の問題を解くための数値計算ライブラリである。 これには連立1次方程式、固有値問題、特異値分解をはじめとして、それらを解くために必要となる多くの補助ルーチンが含まれる。
LAPACK は BLAS のルーチンを使って計算を行うために必ずこれが必要になる。 もちろん REFERENCE BLAS でもよいが、最適化 BLAS を使うことで LAPACK の計算も高速になる。
エルミート行列対角化のためのLAPACKドライバー
LAPACK は 1000 を超えるルーチンからなり、それらを組み合わせて線形代数の諸問題を数値的に計算するためのライブラリである。 多くのルーチンを組み合わせて一般的な問題(連立1次方程式、固有値問題、特異値分解など)を解くためのルーチンは特別にドライバー(driver routine)と呼ばれている。
エルミート行列の対角化には zheev
系のドライバーを用いる。
「系」というのには理由があり、同じ問題でも実装によっていくつかの選択肢が存在するためである。
zheev
- Simple driver
- 最も単純で初期から存在するドライバー
zheevx
- Expert driver
- Simple driver とアルゴリズムは同じだが、固有値を指定の範囲・個数まで求めて計算を打ち切るための仕組みが提供されている
- 固有値が全部必要ない場合には計算が早く終了することが期待できる
- 基本的には Simple driver の上位互換と考えて良い
zheevd
- Divide and conquer アルゴリズムによって固有値問題を解くドライバー
zheevr
- Relatively Robust Representation アルゴリズムによって固有値問題を解くドライバー
- 比較的新しく、v3.10.1 時点ではまだ実装は完全ではないようである
- 具体的に言うと引数 ABSTOL のドキュメントに少し書いてある
- しかし、十分に使えるようではある
- 高速かつ Expert driver と同じく一部の固有値のみを計算することができる
サンプルソースコードでは 1000x1000 のランダムなエルミート行列を対角化するプログラムを示す。 LAPACK はもともと FORTRAN 言語でかかれており、その呼び出し規則に従わなければならない。 つまり、行列は Column major で連続したメモリ上に配置されていなければならない。 C/C++ のプログラムでは1000*1000=1000000要素の1次元配列としてメモリを確保し、Column majorで数値を格納する。
コンピュータのメモリ空間は1次元的な配置になっているので、 2次元の行列をメモリ上に展開するには何らかの規則が必要である。 Column majorとは以下のような行列
を配列Aに次のように配置することである。
Cにおける配列(C++における生配列)はメモリ上に連続的にデータを配置する。 つまり、各列 (column) をメモリ上に連続的に配置している。 Column major と対になる言葉として Row major というものがある。 これは同じ行列を次のように配置することである。
つまり、各行 (row) をメモリ上に連続的に配置している。
下記のプログラムではランダムなエルミート行列を使っているので特に注意されていないが、 本来であれば Column/Row major はLAPACKドライバーを使う上では最も注意が必要な点と思われる。 (今回はエルミート行列なので転地しても固有値は変わらないはずである。)
C/C++ から FORTRAN でかかれたサブルーチンを呼び出すためには、FORTRAN でのサブルーチン名の最後にアンダースコア(_
)をつけなければいけない。
つまり、zheev
ドライバーを呼び出すためには zheev_
と記述する。
また、計算には上三角行列あるいは下三角行列のどちらかしか使わないので、対角化する行列を作る際は全部の要素を埋める必要はない。 つまり、エルミート対称性から対角要素を挟んで逆側の値は決まるので、原理的に計算に必要ないのである。 少しでもプログラムを高速化するためには重要である。
zheevドライバー
void zheev_(const char& JOBZ, const char& UPLO, const int& N, complex* A, const int& LDA, double* W, complex* WORK, const int& LWORK, double* RWORK, int& INFO);
引数 A
に対角化する行列を与えると W
に与えた配列に計算された固有値が格納される。
配列 A
は対角化計算の過程で破壊的変更を受ける。
つまり計算が終了した時点で A
は与えた配列と同一の内容ではないので注意。
JOBZ
に 'V'
を指定した場合は A
の第i列が W
のi番目の要素に対応した固有ベクトルになっている。
WORK
, RWORK
は計算のために必要な作業用のメモリ領域であり、LWORK
は WORK
の配列サイズである。
RWORK
は NxN の行列を対角化する場合、3N-2 のサイズ決め打ちでよい。
WORK
の最適なサイズは LWORK
に -1 を与えて zheev
を呼び出すと計算でき、WORK
の最初の要素に格納される。
つまり、固有値を得るためには 1. LWORK
の決定 → 2. 固有値の計算というように2回 zheev
を呼び出す。
それぞれ最適サイズより大きくても問題ない。
最適サイズより小さい場合、計算速度が悪化することが予想される。
何度も固有値問題を解く場合、WORK
, RWORK
はサイズが十分な限り使い回せる。
すなわち、ブロック対角化のように何回か対角化をする場合、最も大きい行列で WORK
, RWORK
を用意しておけば何度も確保し直す必要はない。
計算終了後に INFO
に計算の終了状態に関する情報が返される。
INFO = 0
なら正常終了。
INFO < 0
なら INFO
番目の引数に問題があり、異常終了。
INFO > 0
なら計算が収束しなかったことを意味する。
zheevxドライバー
void zheevx_(const char& JOBZ, const char& RANGE, const char& UPLO, const int& N, complex* A, const int& LDA, const double& VL, const double& VU, const int& IL, const int& IU, const double& ABSTOL, int& M, double* W, complex* Z, const int& LDZ, complex* WORK, const int& LWORK, double* RWORK, int* IWORK, int* IFAIL, int& INFO);
引数 A
に対角化する行列を与えると W
に与えた配列に計算された固有値が格納される。
配列 A
は対角化計算の過程で破壊的変更を受ける。
つまり計算が終了した時点で A
は与えた配列と同一の内容ではないので注意。
JOBZ
に 'V'
を指定した場合は Z
の第i列が W
のi番目の要素に対応した固有ベクトルになっている。
zheevx
ドライバーは与えられた行列 A
の固有値の内の一部のみを計算する機能を持つ。
すべての固有値が必要ではない場合には、計算を打ち切って高速化することができる。
実際に見つかった固有値の個数は M
に格納される。
RANGE
に'A'
を指定した場合、すべての固有値を計算する。RANGE
に'V'
を指定した場合、(VL
,VU
] の半開区間に入る固有値のみを求める。RANGE
に'I'
を指定した場合IL
番目からIU
番目の固有値のみを求める。
ABSTOL
は固有値の許容計算誤差に関する引数だが、特に理由がなければ最良の精度を得るために 2*DLMCH('S') を与えるのがよい。
WORK
, RWORK
, IWORK
は計算のために必要な作業用のメモリ領域であり、LWORK
は WORK
の配列サイズである。
RWORK
は NxN の行列を対角化する場合、7*N のサイズ決め打ちでよい。
IWORK
は NxN の行列を対角化する場合、5*N のサイズ決め打ちでよい。
WORK
の最適なサイズは LWORK
に -1 を与えて zheevx
を呼び出すと計算でき、WORK
の最初の要素に格納される。
つまり、固有値を得るためには 1. LWORK
の決定 → 2. 固有値の計算というように2回 zheevx
を呼び出す。
それぞれ最適サイズより大きくても問題ない。
最適サイズより小さい場合、計算速度が悪化することが予想される。
何度も固有値問題を解く場合、WORK
, RWORK
, IWORK
はサイズが十分な限り使い回せる。
すなわち、ブロック対角化のように何回か対角化をする場合、最も大きい行列で WORK
, RWORK
, IWORK
を用意しておけば何度も確保し直す必要はない。
計算終了後に INFO
に計算の終了状態に関する情報が返される。
INFO = 0
なら正常終了。
INFO < 0
なら INFO
番目の引数に問題があり、異常終了。
INFO > 0
なら計算が収束しなかったことを意味する。
zheevdドライバー
void zheevd_(const char& JOBZ, const char& UPLO, const int& N, complex* A, const int& LDA, double* W, complex* WORK, const int& LWORK, double* RWORK, const int& LRWORK, int* IWORK, const int& LIWORK, int& INFO);
引数 A
に対角化する行列を与えると W
に与えた配列に計算された固有値が格納される。
配列 A
は対角化計算の過程で破壊的変更を受ける。
つまり計算が終了した時点で A
は与えた配列と同一の内容ではないので注意。
JOBZ
に 'V'
を指定した場合は A
の第i列が W
のi番目の要素に対応した固有ベクトルになっている。
WORK
, RWORK
, IWORK
は計算のために必要な作業用のメモリ領域であり、LWORK
は WORK
の、LRWORK
は RWORK
の、LIWORK
は IWORK
の配列サイズである。
WORK
, RWORK
, IWORK
の最適なサイズは LWORK
, LRWORK
, LIWORK
のどれかに -1 を与えて zheevd
を呼び出すと計算でき、WORK
, RWORK
, IWORK
それぞれの最初の要素に格納される。
つまり、固有値を得るためには 1. LWORK
, LRWORK
, LIWORK
の決定 → 2. 固有値の計算というように2回 zheevd
を呼び出す。
それぞれ最適サイズより大きくても問題ない。
最適サイズより小さい場合、計算速度が悪化することが予想される。
何度も固有値問題を解く場合、WORK
, RWORK
, IWORK
はサイズが十分な限り使い回せる。
すなわち、ブロック対角化のように何回か対角化をする場合、最も大きい行列で WORK
, RWORK
, IWORK
を用意しておけば何度も確保し直す必要はない。
計算終了後に INFO
に計算の終了状態に関する情報が返される。
INFO = 0
なら正常終了。
INFO < 0
なら INFO
番目の引数に問題があり、異常終了。
INFO > 0
なら計算が収束しなかったことを意味する。
zheevrドライバー
void zheevr_(const char& JOBZ, const char& RANGE, const char& UPLO, const int& N, complex* A, const int& LDA, const double& VL, const double& VU, const int& IL, const int& IU, const double& ABSTOL, int& M, double* W, complex* Z, const int& LDZ, int* ISUPPZ, complex* WORK, const int& LWORK, double* RWORK, const int& LRWORK, int* IWORK, const int& LIWORK, int& INFO);
引数 A
に対角化する行列を与えると W
に与えた配列に計算された固有値が格納される。
配列 A
は対角化計算の過程で破壊的変更を受ける。
つまり計算が終了した時点で A
は与えた配列と同一の内容ではないので注意。
JOBZ
に 'V'
を指定した場合は Z
の第i列が W
のi番目の要素に対応した固有ベクトルになっている。
zheevx
ドライバーは与えられた行列 A
の固有値の内の一部のみを計算する機能を持つ。
すべての固有値が必要ではない場合には、計算を打ち切って高速化することができる。
実際に見つかった固有値の個数は M
に格納される。
RANGE
に'A'
を指定した場合、すべての固有値を計算する。RANGE
に'V'
を指定した場合、(VL
,VU
] の半開区間に入る固有値のみを求める。RANGE
に'I'
を指定した場合IL
番目からIU
番目の固有値のみを求める。
ABSTOL
は固有値の許容計算誤差に関する引数である。特に理由がなければ DLMCH('S') を与えるのがよい。
v3.10.1時点ではまだそのような実装になっていないが、将来のバージョンではこの設定で最良の精度を得るようになる。
WORK
, RWORK
, IWORK
は計算のために必要な作業用のメモリ領域であり、LWORK
は WORK
の、LRWORK
は RWORK
の、LIWORK
は IWORK
の配列サイズである。
WORK
, RWORK
, IWORK
の最適なサイズは LWORK
, LRWORK
, LIWORK
のどれかに -1 を与えて zheevr
を呼び出すと計算でき、WORK
, RWORK
, IWORK
それぞれの最初の要素に格納される。
つまり、固有値を得るためには 1. LWORK
, LRWORK
, LIWORK
の決定 → 2. 固有値の計算というように2回 zheevr
を呼び出す。
それぞれ最適サイズより大きくても問題ない。
最適サイズより小さい場合、計算速度が悪化することが予想される。
何度も固有値問題を解く場合、WORK
, RWORK
, IWORK
はサイズが十分な限り使い回せる。
すなわち、ブロック対角化のように何回か対角化をする場合、最も大きい行列で WORK
, RWORK
, IWORK
を用意しておけば何度も確保し直す必要はない。
計算終了後に INFO
に計算の終了状態に関する情報が返される。
INFO = 0
なら正常終了。
INFO < 0
なら INFO
番目の引数に問題があり、異常終了。
INFO > 0
なら計算が収束しなかったことを意味する。
コンパイル
Windows
Windows では MSYS2 を使うのが楽。
MSYSの MinGW 64-bit 環境を使う想定で、コンパイルには gcc の c++ コンパイラ、g++
を使う。
LAPACK、REFERENCE BLAS のインストールは以下のコマンドでできる。
REFERENCE BLAS は LAPACK のパッケージに同梱されている。
$ pacman -S mingw-w64-x86_64-lapack
OpenBLAS のインストールは以下のコマンドでできる。 LAPACK のドライバー(後述)も openblas のバッケージに含まれている。
$ pacman -S mingw-w64-x86_64-openblas
REFERENCE BLAS か OpenBLAS かについては基本的に OpenBLAS を使えばよい。
コンパイルオプションとして -lopenblas
を指定すると OpenBLAS がリンクされる。
つまり、例えば次のようにすると実行ファイル a.out
が生成される。
$ g++ example_zheev.cpp -o a.out -lopenblas
REFERENCE LAPACKが openblas に静的リンクされているみたいなので、-llapack
はいらない。むしろ指定順序によってはつけるとLAPACKに静的リンクされている REFERENCE BLAS が優先されてしまうのでつけてはいけない。
あまり、使う理由はないが REFERENCE BLAS を使いたい場合は -llapack
を代わりに使う。REFERENCE BLAS が静的リンクされているようなので、-lblas
は(MSYSのパッケージを使う場合は)必要ない。
$ g++ example_zheev.cpp -o a.out -llapack
ちなみにコンパイルされた実行ファイルがどちらに依存しているかは ldd
コマンドで確認することができる。
$ ldd a.out
Ubuntu
Ubuntu では apt
コマンドを使って、REFERENCE BLAS、OpenBLAS をインストールできる。
LAPACK (REFERENCE BLAS) は次のコマンド、
$ apt install liblapack-dev
OpenBLAS は次のコマンドを実行するとよい。
$ apt install libopenblas-dev
また、MKL は Intel 公式の手順 に従うことで容易に使用可能となる。
複数の lapack 互換ライブラリがインストールされている場合、-llapack
フラグでリンクされるライブラリを update-alternatives
コマンドを使って選択することができる。
$ update-alternatives --get-selections | grep liblapack.so
とすると現在の状態を調べることができる。複数表示される場合は適当なソースを実際にコンパイルして ldd
コマンドで確認するのが確実と思われる。これはもっと良い方法がありそうだが…
$ g++ example_zheev.cpp -o a.out -llapack $ ldd a.out (略) liblapack.so.3 => /lib/x86_64-linux-gnu/liblapack.so.3 (0x00007fbab32f3000) (略)
liblapack.so.3
がリンクされるようである。update-alternatives
コマンドで該当する選択肢を探す。
$ update-alternatives --get-selections | grep liblapack.so.3 liblapack.so.3-x86_64-linux-gnu manual /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3
この場合、liblapack.so.3-x86_64-linux-gnu
が識別子のようである。次のコマンドを実行すると -llapack
で使用されるライブラリを選択することができる。
$ sudo update-alternatives --config liblapack.so.3-x86_64-linux-gnu
あるいは -llapack
フラグの代わりに -lopenblas
、-lmkl_rt
を使うと必ず OpenBLAS、MKL が使われる。
こちらの方が簡単かもしれない。
- OpenBLAS を使う例
$ g++ example_zheev.cpp -o a.out -lopenblas
パフォーマンス比較
使用した環境は以下の通り。
- Ubuntu 20.04
- gcc 9.4.0 (Ubuntu 9.4.0-1ubuntu1~20.04.1)
- AMD Ryzen7 3700X (8 core, 16 threads)
- LAPACK v3.9.0
- REFERENCE BLAS v3.9.0
- OpenBlas v0.3.8
1スレッドに制限した場合の結果は以下。
使用されているBLASの比較をするとOpenBLASの成績が良い。 MKLも使っては見たものの、AMDのCPUを使っているので最適化度合いはイマイチか。
ドライバー毎の比較をすると zheevr
の成績が良い。
zheevx
と zheevd
はほぼ同等。
MKLを除いて zheev
は他より2倍以上遅いようである。
さらに OpenBLAS を使用した場合についてマルチスレッドの効果を示したのが次の図である。
zheevd
はマルチスレッド化の効果が大きく、4スレッドあたりで zheevr
に逆転している。
全体的には zheevd
と zheevr
の成績が良い。
マルチスレッド化の効果が大きければ zheevd
を使うのがよいだろう。
ただし、zheevr
は求める固有値の数を制限する機能があるなど柔軟なので、問題によってはよいだろう。
どちらにせよ zheevd
と zheevr
の差は小さいのでどちらでも良い気がする。
REFERENCE BLAS/1thread/zheev の場合に比べて OpenBLAS/8thread/zheevd は20倍近い高速化をしている。
CLAPACK と LAPACKE
今回は FORTRAN のサブルーチンを C++ から呼び出す方式をとったが、他のLAPACKを使うためのインターフェースとして CLAPACK と LAPACKE というものがある。
LAPACKはもともと FORTRAN で書かれたプログラムなのだが、CLAPACKはこれを f2c でC言語に自動変換したもの。 正直なところ、普通にCからFORTRANのサブルーチンを呼べるのであまり何が便利なのかわかっていない。
LAPACKE はもともと MKL で使われていた LAPACK をラップしたインターフェースらしい。
現在では本家LAPACKに逆輸入されて公式に同梱・サポートされている。
行列の引数に対しては Column/Row major を指定することができ、必要に応じて再配置してくれるようである。
lapacke_{ドライバー名}
とlapacke_{ドライバー名}_work
の2つの系統の関数が提供されており、前者は作業用のメモリ領域を自動で確保してくれるので便利。
後者は Column/Row major の指定以外はLAPACKのネイティブインターフェースを踏襲している。
lapacke_{ドライバー名}
は便利であるが、究極的に速度を求める場合は結局、一度確保した作業用メモリを使いまわしたりする必要があるので結局、lapacke_{ドライバー名}_work
が必要になるのだと思われる。