FreeFEM 4.6 で対称行列をPARDISOで解く [FreeFEM]
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
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
2020-11-26 20:05
コメント(0)
FreeFEM 4.6 でMUMPSで並列計算で求解する [FreeFEM]
FreeFEMのDocumentationのMUMPSに関する記述は (少なくともRelease4.6では) 古いモジュール "MUMPS_FreeFem" に関するものなので、最新のモジュール "MUMPS" の使い方はExample (examples/mpi/MUMPS.edp) を参照すること、という話。
varf lap(u,v) = ... として弱形式を記述して、matrix A = lap(Vh,Vh) 等として全体剛性行列を得る。このあと、load "MUMPS" を実行した上でset (A, solver=sparsesolver, master=xx) とすればソルバーとして直接法のMUMPSを使用するよう設定でき、x[] = A^-1 * b[] 等として求解できる。逐次版でもOpenMPを用いた並列化 (少なくともWindows版の"MUMPS_seq"ではそのようにビルドされているよう) でマルチコアで計算できるし、MPI版ならもちろん並列で求解できる。
ただMPI版で気になるのが、全体剛性行列はどのような値が用いられるのか。結論としては、set()の引数masterの値で動作を切り換えられる。
1) master=0の場合 (centralized)
係数行列はランク0のプロセスの行列Aの値を使う。他のプロセスの値は使用しない。このとき、MUMPSにはICNTL(18)=0 (the input matrix is centralized on the host) が与えられる。また、連立一次方程式の右辺bも同様に、ランク0のプロセスで与えた値を用いる。解xもランク0のプロセスのみが持つ。
つまり、ランク0が逐次版と全く同様に動作するように作成しておけば、求解部分のみを並列計算で行ってくれる。シンプル。
2) master=-1の場合 (distributed)
係数行列は全てのプロセスの (ローカルな) 行列Aの値の総和となる。このとき、MUMPSにはICNTL(18)=3 (user directly provides the distributed matrix, pattern and entries, ...) が与えられる。また、連立一次方程式の右辺bも同様に、全プロセスの値の合計を用いる (FreeFEMのプラグイン側でMPI_Reduceで求めてMUMPSに渡す)。解xも全てのプロセスにブロードキャストされる。
Example中のdistributed版では、change(Th, fregion=nuTriangle%mpisize) でregion番号をプロセス数を上限とする番号に付け直し、int2d(Th, mpirank) (...) としてregion番号が自分のランクと等しい要素のみをローカルな行列Aに計算している。FreeFEM側で各プロセスは自分の持ち分のみ行列要素を計算し、それを足し合わせる処理はMUMPSに任せていることになる。なるほど。
他にもMUMPSのset() のオプションにはstrategy, rinfoやinfoなどがあるようだけど、そちらは未調査です。
※以上、FreeFEM v4.6時点の話です。次にAPIが大変更されるまでは適用可能と思いますが。
参考：plugin/mpi/MUMPS.cpp
MUMPS Documentation
11月23日 17時10分
fartrip
varf lap(u,v) = ... として弱形式を記述して、matrix A = lap(Vh,Vh) 等として全体剛性行列を得る。このあと、load "MUMPS" を実行した上でset (A, solver=sparsesolver, master=xx) とすればソルバーとして直接法のMUMPSを使用するよう設定でき、x[] = A^-1 * b[] 等として求解できる。逐次版でもOpenMPを用いた並列化 (少なくともWindows版の"MUMPS_seq"ではそのようにビルドされているよう) でマルチコアで計算できるし、MPI版ならもちろん並列で求解できる。
ただMPI版で気になるのが、全体剛性行列はどのような値が用いられるのか。結論としては、set()の引数masterの値で動作を切り換えられる。
1) master=0の場合 (centralized)
係数行列はランク0のプロセスの行列Aの値を使う。他のプロセスの値は使用しない。このとき、MUMPSにはICNTL(18)=0 (the input matrix is centralized on the host) が与えられる。また、連立一次方程式の右辺bも同様に、ランク0のプロセスで与えた値を用いる。解xもランク0のプロセスのみが持つ。
つまり、ランク0が逐次版と全く同様に動作するように作成しておけば、求解部分のみを並列計算で行ってくれる。シンプル。
2) master=-1の場合 (distributed)
係数行列は全てのプロセスの (ローカルな) 行列Aの値の総和となる。このとき、MUMPSにはICNTL(18)=3 (user directly provides the distributed matrix, pattern and entries, ...) が与えられる。また、連立一次方程式の右辺bも同様に、全プロセスの値の合計を用いる (FreeFEMのプラグイン側でMPI_Reduceで求めてMUMPSに渡す)。解xも全てのプロセスにブロードキャストされる。
Example中のdistributed版では、change(Th, fregion=nuTriangle%mpisize) でregion番号をプロセス数を上限とする番号に付け直し、int2d(Th, mpirank) (...) としてregion番号が自分のランクと等しい要素のみをローカルな行列Aに計算している。FreeFEM側で各プロセスは自分の持ち分のみ行列要素を計算し、それを足し合わせる処理はMUMPSに任せていることになる。なるほど。
他にもMUMPSのset() のオプションにはstrategy, rinfoやinfoなどがあるようだけど、そちらは未調査です。
※以上、FreeFEM v4.6時点の話です。次にAPIが大変更されるまでは適用可能と思いますが。
参考：plugin/mpi/MUMPS.cpp
MUMPS Documentation
11月23日 17時10分
fartrip
2020-11-23 17:11
コメント(0)
