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とは以下のような行列

 
\begin{matrix}
a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{matrix}

を配列Aに次のように配置することである。

 
A[0] = a_{11} \\
A[1] = a_{21} \\
A[2] = a_{31} \\
A[3] = a_{12} \\
A[4] = a_{22} \\
A[5] = a_{32} \\
A[6] = a_{13} \\
A[7] = a_{23} \\
A[8] = a_{33}

Cにおける配列(C++における生配列)はメモリ上に連続的にデータを配置する。 つまり、各列 (column) をメモリ上に連続的に配置している。 Column major と対になる言葉として Row major というものがある。 これは同じ行列を次のように配置することである。

 
A[0] = a_{11} \\
A[1] = a_{12} \\
A[2] = a_{13} \\
A[3] = a_{21} \\
A[4] = a_{22} \\
A[5] = a_{23} \\
A[6] = a_{31} \\
A[7] = a_{32} \\
A[8] = a_{33}

つまり、各行 (row) をメモリ上に連続的に配置している。


下記のプログラムではランダムなエルミート行列を使っているので特に注意されていないが、 本来であれば Column/Row major はLAPACKドライバーを使う上では最も注意が必要な点と思われる。 (今回はエルミート行列なので転地しても固有値は変わらないはずである。)

C/C++ から FORTRAN でかかれたサブルーチンを呼び出すためには、FORTRAN でのサブルーチン名の最後にアンダースコア(_)をつけなければいけない。 つまり、zheev ドライバーを呼び出すためには zheev_ と記述する。

また、計算には上三角行列あるいは下三角行列のどちらかしか使わないので、対角化する行列を作る際は全部の要素を埋める必要はない。 つまり、エルミート対称性から対角要素を挟んで逆側の値は決まるので、原理的に計算に必要ないのである。 少しでもプログラムを高速化するためには重要である。

zheevドライバー

ドキュメント: https://netlib.org/lapack/explore-html/df/d9a/group__complex16_h_eeigen_gaf23fb5b3ae38072ef4890ba43d5cfea2.html

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列が Wi番目の要素に対応した固有ベクトルになっている。

WORK, RWORK は計算のために必要な作業用のメモリ領域であり、LWORKWORK の配列サイズである。 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ドライバー

ドキュメント: https://netlib.org/lapack/explore-html/df/d9a/group__complex16_h_eeigen_gaabef68a9c7b10df7aef8f4fec89fddbe.html

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列が Wi番目の要素に対応した固有ベクトルになっている。

zheevx ドライバーは与えられた行列 A の固有値の内の一部のみを計算する機能を持つ。 すべての固有値が必要ではない場合には、計算を打ち切って高速化することができる。 実際に見つかった固有値の個数は M に格納される。

  • RANGE'A' を指定した場合、すべての固有値を計算する。
  • RANGE'V' を指定した場合、(VL, VU] の半開区間に入る固有値のみを求める。
  • RANGE'I' を指定した場合 IL 番目から IU 番目の固有値のみを求める。

ABSTOL は固有値の許容計算誤差に関する引数だが、特に理由がなければ最良の精度を得るために 2*DLMCH('S') を与えるのがよい。

WORK, RWORK, IWORK は計算のために必要な作業用のメモリ領域であり、LWORKWORK の配列サイズである。 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ドライバー

ドキュメント: https://netlib.org/lapack/explore-html/df/d9a/group__complex16_h_eeigen_ga9b3e110476166e66f2f62fa1fba6344a.html

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列が Wi番目の要素に対応した固有ベクトルになっている。

WORK, RWORK, IWORK は計算のために必要な作業用のメモリ領域であり、LWORKWORK の、LRWORKRWORK の、LIWORKIWORK の配列サイズである。 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ドライバー

ドキュメント: https://netlib.org/lapack/explore-html/df/d9a/group__complex16_h_eeigen_ga60dd605c63d7183a4c289a4ab3df6df6.html

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列が Wi番目の要素に対応した固有ベクトルになっている。

zheevx ドライバーは与えられた行列 A の固有値の内の一部のみを計算する機能を持つ。 すべての固有値が必要ではない場合には、計算を打ち切って高速化することができる。 実際に見つかった固有値の個数は M に格納される。

  • RANGE'A' を指定した場合、すべての固有値を計算する。
  • RANGE'V' を指定した場合、(VL, VU] の半開区間に入る固有値のみを求める。
  • RANGE'I' を指定した場合 IL 番目から IU 番目の固有値のみを求める。

ABSTOL は固有値の許容計算誤差に関する引数である。特に理由がなければ DLMCH('S') を与えるのがよい。 v3.10.1時点ではまだそのような実装になっていないが、将来のバージョンではこの設定で最良の精度を得るようになる。

WORK, RWORK, IWORK は計算のために必要な作業用のメモリ領域であり、LWORKWORK の、LRWORKRWORK の、LIWORKIWORK の配列サイズである。 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およびドライバーの比較

使用されているBLASの比較をするとOpenBLASの成績が良い。 MKLも使っては見たものの、AMDのCPUを使っているので最適化度合いはイマイチか。

ドライバー毎の比較をすると zheevr の成績が良い。 zheevxzheevd はほぼ同等。 MKLを除いて zheev は他より2倍以上遅いようである。

さらに OpenBLAS を使用した場合についてマルチスレッドの効果を示したのが次の図である。 zheevd はマルチスレッド化の効果が大きく、4スレッドあたりで zheevrに逆転している。

スレッド数による効果(OpenBLAS)

全体的には zheevdzheevr の成績が良い。 マルチスレッド化の効果が大きければ zheevd を使うのがよいだろう。 ただし、zheevr は求める固有値の数を制限する機能があるなど柔軟なので、問題によってはよいだろう。 どちらにせよ zheevdzheevr の差は小さいのでどちらでも良い気がする。

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 が必要になるのだと思われる。