FreeFEM 4.6 で対称行列をPARDISOで解く

FreeFEMはかねてから直接法ソルバーのIntel MKL PARDISOにプラグインで対応しているのですが、細かな対応具合はバージョンでまちまちでした。ただFreeFEM 4.6では行列の対称性としてStructually symmetric以外のPARDISOに入力可能なもの全てに対応していたので、備忘録としてメモします。

1. Linuxであれば、MKLを利用してビルドするとplugin/seq/PARDISO.cpp をビルドしてモジュールを作成してくれる。configureのオプションに--with-mkl=(mklのインストールパス)/linux/mkl/lib/intel64_lin と --enable-mkl-mlt をつければ良い。後者はOpenMP並列化でマルチコアで実行するために必要。Windowsであればプラグインの修正（自作）が必要？
2. edpスクリプト内でvarf lap(u,v)=... として弱形式を記述して、matrix A = lap(Vh,Vh) 等として全体剛性行列を得る。
3. load "PARDISO" を実行した上で、set(A, solver=sparsesolver, sym=xx, positive=xx) としてPARDISOを利用するよう指定した後、x[] = A^-1 * b[] などとして求解できる。

ここで重要なのが、set()の引数のsymとpositive。例によってDocumentationにはまだ記述が見つからず、Exampleのexamples/plugin/PARDISO.edp から使い方を見ました。行列の対称性により、下記の組み合せで指定するようです。

A. 実数非対称 sym=0 もしくはsym省略
B. 実数対称不定値 sym=1, positive=false もしくはpositive省略
C. 実数対称正定値 sym=1, positive=true
D. 複素数非対称 sym=0 もしくはsym省略
E. 複素数エルミート不定値 sym=1, positive=false もしくはpositive省略
F. 複素数エルミート正定値 sym=1, positive=true
G. 複素数対称非エルミート sym=2

特にsym=2の記述は別のサンプル(examples/ffddm/Helmholtz-3d-simple.edp)でようやく見つけました。

非対称(sym=0)のときは、行列Aをわざわざ明示的に作成しなくても、solve lap(u,v,solver=sparsesolver)=... などとしてもPARDISOで解くことができます。またproblem, solveのinit引数も有効で、例えばsolve lap(u,v,solver=sparsesolver, init=1)=... とすれば、最初の1回のみfactorization含め全て行い (phase 1-3)、以降は行列部分の再評価無しにsolveのみ (phase 3) を行ってくれます。

また、verbosityを5以上に上げればPARDISOに渡す引数のmsglvl=1となり、進捗や所要時間、使用コア数等の情報を印字してくれます。

速度はMUMPSとどっこいですが、メモリ使用量は気持ち少ないようです。MKLは無料で入手でき、MUMPSと違いMPIと関わらないことから、ビルドさえ成功すれば気軽に利用できる直接法ソルバーと思います。

11月26日20時3分
fartrip