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)
コメント 0