nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

MKL11.3:逆行列計算ルーチンの罠

自分の今の環境では、MKL11.3を使ってfortran90でプログラムを弄っている。違う環境では、ここで言うことは当てはまらないかも知れない。

逆行列を計算したい。私の願いはそれだけであった。

一般の行列に対して対応可能なsubroutineは、GETRF (LU分解)とGETRI(それを踏まえて逆行列計算)である。頭文字のGはGeneralのGであろう。両方セットで使う。引数は、変換したい行列a(*)、変換するときの情報が入っているipiv (integer)、変換が上手く行ったか教えてくれるinfo (integer)、の三つのみである。aは実数でも倍精度実数でも複素数でも倍精度複素数でも、GETRFやGETRIは内部でそれぞれの型に対応した関数に渡して処理してくれる。
実際、GETRFやGETRIはLAPACK95のモジュールを使用しているのだが、これはInterfaceの集まりであり、計算そのものを実行しているのはF95_PRECISIONの中のCGETRF(単精度複素数用)、ZGETRF(倍精度複素数用)といった、型がキッチリ指定されたやつである。使うときには、LAPACK95のみをuseすると宣言してGETRF, GETRIをただ使えば十分である。

罠からは外れるが、使うときにはlapack95.modがあれば良く、f95_precision.modは無くて大丈夫。
ややこしい設定ファイルは、Windowsだと/Qmklとか/QMkl:sequentialとか、Linuxだと-mklで通ると書いてあるが、僕の環境では通らなかったので、自力でファイルを持って来たりした。使ったファイルだけ列挙すると、
libiomp5md.lib
mkl_core.lib
mkl_intel_ilp64.lib
mkl_intel_thread.lib
mkl_lapack95_ilp64.lib
である。libiomp5md.libは mkl/lib/intel64_win/ の中ではなく、compiler/lib/intel64_win/ に入っているため、探すのに苦労した気がする。
後、ilp64のファイルを使う場合には、コンパイルオプションに /4I8 を追加する必要がある。
本来は、Makefileをキチっと作って使い回すものなんでしょうけれども、Makefileはまだ勉強中であります。

んで、何が罠だったかというと、自分がやりたいのはエルミート行列の逆行列計算なので、エルミート行列用のHETRFとHETRIが使いたかった。頭文字のHはもちろんHermitianのHであろう。これまでは適当にGETRF,GETRIを使っていて、せっかくaがエルミートだから速く計算できる方が良かろうということであるのだが、引数の並びに問題があった。

call GETRF( a, [ipiv], [info] )
call GETRI( a, [ipiv], [info] )

に対して、

call HETRF( a, [uplo], [ipiv], [info] )
call HETRI( a, ipiv, [uplo], [info] )

である。uploは、Bunch-Kaufman diagonal pivoting methodにおいて上三角行列と下三角行列どちらを使うかというもので、'U'か'L'が入る。[...]はOPTIONAL引数で、省略可能である。
この引数の順番の違いで、コンパイルが全然通らなくてだいぶ困った。多分エラーメッセージが良くない。
「この汎用サブルーチン呼び出しと一致する特殊サブルーチンがありません。」
って言われたら、真っ先に設定ファイルが足りてないことを疑うのではないのだろうか?設定ファイルの多さと、GETRF,GETRIの流れから引数を疑うまでになかなか至らなかった。。。
もしMKLで詰まったら、大元のlapack.f90を見てみて下さい(というか、本来はそうするのが当たり前なんだろうけれども。。。)。