概要

Parallel Sparse Direct Solver(PARDISO)は大規模な疎行列連立一次方程式を解く際に用いられるソルバーです．
物理シミュレーションのコードなどを書くとき，最終的に$Ax=b$の形の連立一次方程式を解くという問題に帰着されることはよくあります．いままでは秘伝のタレの奥深くでいい感じに謎の引数たちを渡してこの連立一次方程式を解いているということだけ把握していたのですが，自作のプログラムを実装する際にどのように使えば連立一次方程式を解けるのかわからなかったので，簡単な例題を解くことでPARDISOの使いかたを確認しました．

上記の目的から，本記事の内容はPARDISOの性能をフルに活用できるようになるための知識は提供しません．そのかわりに本記事を読むことで，PARDISOを使ったことのない人も連立一次方程式を解けるようになります．詳細かつ正確な解説は公式ドキュメントを参照してください．また，Panua TechnologyのホームページにはFortran以外のプログラミング言語(C, C++, Python, MATLAB)を用いたプログラムの実装例も掲載されています(https://www.panua.ch/pardiso/)．

問題設定

解く方程式は次式とします．
$$\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$$
ただし，$\boldsymbol{A}$は$n$次の実正方行列で，$\boldsymbol{x}$および$\boldsymbol{b}$はそれぞれ$n$次の列ベクトルとします．
この方程式を$\boldsymbol{x}$について解きます．

PARDISOを利用して解くことができるのは，$\boldsymbol{A}$が以下のいずれかに該当する場合です．

実構造対称
実対称かつ正定値
実対称かつ不定値
複素構造対称
複素エルミート正定値
複素エルミート不定値
複素対称
実非対称
複素非対称
ここで，構造対称(structually symmetric)行列とは，$a(i,j)\neq 0$ならば$a(j,i)\neq 0$が成り立つような行列です($a(i,j)=a(j,i)$とは限らない)．

$\boldsymbol{A}$が上記のどれに該当するかは，ソルバーを呼び出す前に指定しておく必要があります．

連立一次方程式をPARDISOを使って解く

PARDISOを利用して連立一次方程式を解くFortranプログラムを書いていきます．
この記事を書いている人はこちらの記事¹を参考に環境構築を行いました．

ここではとにかく連立一次方程式を解くことを目標とするため，必ずしも把握しなくてよいようなパラメータについては深く追いません．

また，最後にプログラム全体をまとめて再掲するので，先に全体を見たい方はそちらを参照してください．

Fortranではプログラム冒頭で変数宣言をする必要があります．今回のプログラム冒頭はこんな感じになります．

変数宣言

program main
   implicit none
   !PARDISO関連
   integer::pt(64)
   integer::maxfct
   integer::mnum
   integer::mtype               !matrix type
   integer::solver              !solver type
   integer::phase              
   integer,allocatable::perm(:)
   integer ::iparm(64)
   integer ::msglvl
   integer ::error
   double precision::dparm(64) 

   !行列とベクトル　行列はCSR形式で格納
   integer::n         !行列の次元
   integer::nnz       !行列の非零要素の数
   integer::nrhs
   integer,         allocatable::ia(:)!各行最初の非零要素の位置
   integer,         allocatable::ja(:)!非零要素の列インデックス
   double precision,allocatable::A(:) !行列の非零要素
   double precision,allocatable::b(:) !右辺ベクトル
   double precision,allocatable::x(:) !解ベクトル
   integer         ::i

PARDISOを利用する準備

PARDISOを利用して連立一次方程式を解くというときに，実際に方程式を解くのは pardisoサブルーチンですが，pardisoサブルーチンを呼び出す前に必要な準備がいくつかあります．ざっくりと分類すると，以下のようになります．

行列タイプの指定
使用するソルバーの指定
ptおよび制御パラメータの初期化

この分類には何か明確な基準があるわけではありません．
一つ目は解く問題によって設定する必要がある項目，二つ目は自由に選択できる項目，三つ目はデフォルトの値のままあまり気にしなくてよいもの，程度の分類です．

行列タイプの指定

冒頭で，PARDISOを利用することのできる$\boldsymbol{A}$の条件を列挙しました．プログラム中で

行列タイプの指定

    mtype = 11

のようにすることで，$\boldsymbol{A}$がどのような行列かを指定することができます．
mtypeの値と行列の性質の対応表を以下にまとめます．

mtype	対応する行列
1	実構造対称
2	実対称正定値
-2	実対称かつ不定値
3	複素構造対称
4	複素エルミートかつ正定値
-4	複素エルミートかつ不定値
6	複素対称
11	実非対称
13	複素非対称

ここでは，複数の例題を解く際にいちいち行列の性質を確認するのが面倒くさいためmtypeの値を11とします．

使用するソルバーの指定

続いて，ソルバーの指定をします．solverに代入する値によってどのように解を求めるかを指定することができます．

ソルバーの指定

    solver = 1

選べるソルバーは以下の二つです²．

solver	対応する解法
0	sparse direct solver
1	multi-recursive iterative solver

疎行列直接法ソルバと反復法ソルバについての詳しい説明は文献³を参照してください．ここではsolver=1を選択します．

`pt`および制御パラメータの初期化

pardisoinitサブルーチンはpt，iparmおよびdparmを初期化します．

pardisoinit

    call pardisoinit(pt, mtype, solver, iparm, dparm, error)
    if(error/=0)then
        if(error==-10) print *, "No license file pardiso.lic found"
        if(error==-11) print *, "License is expired"
        if(error==-12) print *, "Wrong username or hostname"
        stop
    endif
    iparm(3)=1
    msglvl  =1
    maxfct  =1
    mnum    =1

errorの値が0であればpardisoinitはうまく動作したことを意味します．

ptは内部データのアドレスに関する変数です．この変数は決して変更してはいけないと公式のユーザーガイドに記されています．

iparmとdparmはPARDISOに関する様々なパラメータを格納している配列です．pardisoinitではiparm(1)，iparm(2)そしてiparm(4)~iparm(64)，dparm(1)~dparm(64)にデフォルトの値を代入します．このデフォルトの値については公式ユーザーガイドのsection 2.3およびsection 2.4を参照してください．

iparm(3)はプロセッサ数を代入します．これは使用する計算機によって値が異なるので自分で設定する必要があります．

分散メモリ型計算機を使用する場合はiparm(52)の値も設定する必要があります．

msglvlはソルバーからメッセージを出力するかを設定することができます．

msglvl	対応する動作
0	なにも出力しない
1	ソルバーから統計的な情報を標準出力する．ただし，multi-recursive iterative solverを利用する場合は`pardiso-ml.out`というファイルに出力される

maxfctとmnumは，多くの場合1に設定しておけばよいようです．公式ユーザーガイドによると，それぞれ変数の説明は以下の通りになっています．

"Maximal number of factors with identical nonzero sparsity structure that the user would like to keep at the same time in memory. It is possible to store several different factorizations with the same nonzero structure at the same time in the internal data management of the solver."

"Actual matrix for the solution phase. With this scalar the user can define the matrix that he would like to factorize. The value must be: $1\leq \mathrm{MNUM}\leq \mathrm{MAXFCT}$"

連立一次方程式の設定

今回はこんな$\boldsymbol{A}$と$\boldsymbol{b}$を用意してみました．
その他の例題も記事の末尾に載せようと思います．

問題設定

   !------------------------------------------------問題設定
   !Ax=b
   !A:
   !    2  -1  0 
   !   -1   2 -1 
   !    0  -1  2 
   !b:
   !    1 
   !    0 
   !   -1 
   !x:
   !    0.5 
   !    0
   !   -0.5

   !行列の次元
   n   = 3
   !右辺の数
   nrhs= 1
   !行列の非零要素の数
   nnz = 7

   !!CSR形式で行列を格納
   allocate(ia(n+1))
   allocate(ja(nnz))
   allocate( A(nnz))
   
   !ia:各行最初の非零要素の位置
   ia = (/ 1, 3, 6, nnz+1/)
   !ja:非零要素の列インデックス
   ja = (/ 1, 2, 1, 2, 3, 2, 3/)
   !A:行列の非零要素
   A  = (/2.d0,-1.d0,-1.d0,2.d0,-1.d0,-1.d0,2.d0/)

   !右辺ベクトルを設定
   allocate(b(n,nrhs))
   b(:,1) = (/1.d0,0.d0,-1.d0/)

   !解ベクトルを設定
   allocate(x(n,nrhs))
   !-------------------------------------------------------

n，nrhs，nnzはそれぞれ行列の次元，右辺の数，行列の非零要素の数です．

PARDISOでは，行列$\boldsymbol{A}$の格納形式としてCompressed Sparse Row(CSR)形式(CRS形式とも)呼ばれる方式を使用します．CSR形式は，疎行列の非零要素のみを記憶するようなデータ形式であり，3つの一次元配列ia，jaおよびAで疎行列$\boldsymbol{A}$を表現します．CSR形式について詳細は例えば記事⁴⁵を参照してください．PARDISOではiaの最後の要素にはnnz+1が格納されることに注意が必要です．

`pardiso`を呼び出して連立一次方程式を解く

いよいよpardisoサブルーチンを呼び出して連立一次方程式を解きます．

solver

   phase= 13
   call pardiso(pt, maxfct, mnum, mtype, phase, n, A, ia, ja, perm, nrhs, iparm, msglvl, b, x, error, dparm)

ここで新たに登場する変数はphase，permですが，permはiparm(5)=1と設定しない限り気にする必要はありません．iparm(5)=1とする場合については公式のユーザーガイドを参照してください．

文献³でも概説されているように，pardisoも解析フェーズ，数値分解フェーズ，求解フェーズの3つのフェーズに分かれています．これらのフェーズを順にすべて処理することで連立一次方程式を解くことができます．一方で，phaseを適当な値に設定することで，一部のフェーズのみ処理させることも可能です．

phaseは2桁の整数で，十の位と一の位の数字によってpardisoサブルーチンの動作を指定します．詳細は公式のユーザーガイドを参照してください．連立一次方程式を解くためにはすべてのフェーズを順に処理するように，phase=13としておけばよいです．

呼び出し後に返ってくるerrorの値が0であれば連立一次方程式の求解に無事成功しているはずです．

解の確認

check

   if(error/=0)then
      print *, "Pardiso error. error:",error
      stop
   endif
   
   print *, "x"
   print *, x
endprogram main

以上のプログラムを実行すると，

 x
  0.500000000000000       0.000000000000000E+000 -0.500000000000000

と表示されてきちんと解を計算できていることがわかりました．

おわりに

本記事ではPARDISOを利用してとりあえず連立一次方程式を解けるようになることを目的として，最低限気を配らなければならないパラメータのみに注目しました．プログラムの試作段階ではソルバーの性能を最大限活用するよりも，とりあえず解けるようにして次の段階へ進みたいことがあると思います．この記事がそのような方の役に立てれば嬉しいです．

また，今回が初めてのQiitaへの投稿になります．もし記事の内容に何か誤りや不適切な点がありましたら教えて頂けると幸いです．

プログラム

今回作成したプログラムをまとめたものがこちらになります．

main.f90

program main
   implicit none
   !PARDISO関連
   integer::pt(64)    !PARDISO ハンドル
   integer::maxfct
   integer::mnum
   integer::mtype
   integer::solver
   integer::phase 
   integer,allocatable::perm(:)
   integer ::iparm(64)!制御パラメータ
   integer ::msglvl
   integer ::error
   double precision::dparm(64) 

   !行列とベクトル　行列はCSR形式で格納
   integer::n         !行列の次元
   integer::nnz       !行列の非零要素の数
   integer::nrhs
   integer,         allocatable::ia(:)!各行最初の非零要素の通し番号
   integer,         allocatable::ja(:)!非零要素の列インデックス
   double precision,allocatable::A(:) !行列の非零要素
   double precision,allocatable::b(:,:) !右辺ベクトル
   double precision,allocatable::x(:,:) !解ベクトル 
   integer::i
   
   !matrix type
   mtype  =  11! 
   solver =  1 !

   !pt, PARDISOの制御パラメータを初期化
   call pardisoinit(pt, mtype, solver, iparm, dparm, error)
   if(error/=0)then
      if(error==-10) print *, "No licence file pardiso.lic found"
      if(error==-11) print *, "License is expired"
      if(error==-12) print *, "Wrong username or hostname"
      stop
   endif   
   !number of processors (value of omp_num_threads)
   iparm(3)=1

   !msglvl: 0:no output, 1 output
   msglvl  =0
   
   
   !maxfct:maximal number of factors 
   !       with identical nonzero sparsity structure 
   !       that the user would like to keep 
   !       at the same time in memory
   maxfct  =1 !In many of applications mxfct=1

   !mnum  :Actual matrix for the solution phase.
   !       With this scalar the user can define the matrix
   !       that he would like to factorize.
   !       (1<=mnum<=maxfct)
   mnum    =1 !In many of application mnum=1   
   
   do i=1,64
      print *, "iparm(",i,"):",iparm(i)
   enddo


   !------------------------------------------------問題設定
   !Ax=b
   !A:
   !    2  -1  0 
   !   -1   2 -1 
   !    0  -1  2 
   !b:
   !    1 
   !    0 
   !   -1 
   !x:
   !    0.5 
   !    0
   !   -0.5

   !行列の次元
   n   = 3
   !右辺の数
   nrhs= 1
   !行列の非零要素の数
   nnz = 7

   !!CSR形式で行列を格納
   allocate(ia(n+1))
   allocate(ja(nnz))
   allocate( A(nnz))
   
   !ia:各行最初の非零要素の位置
   ia = (/ 1, 3, 6, nnz+1/)
   !ja:非零要素の列インデックス
   ja = (/ 1, 2, 1, 2, 3, 2, 3/)
   !A:行列の非零要素
   A  = (/2.d0,-1.d0,-1.d0,2.d0,-1.d0,-1.d0,2.d0/)

   !右辺ベクトルを設定
   allocate(b(n,nrhs))
   b(:,1) = (/1.d0,0.d0,-1.d0/)

   !解ベクトルを設定
   allocate(x(n,nrhs))
   !-------------------------------------------------------

   !Solve
   phase= 13
   call pardiso(pt, maxfct, mnum, mtype, phase, n, A, ia, ja&
               , perm, nrhs, iparm, msglvl, b, x, error, dparm)

   
   if(error/=0)then
      print *, "pardiso error. error:",error
      stop
   endif
   
   print *, "x"
   print *, x


endprogram main

その他の例題

疎行列とは言えませんが．記事⁶に載っていた問題です．

例題

   !------------------------------------------------問題設定
   !Ax=b
   !A:
   !   2   1 -3 -2
   !   2  -1 -1  3
   !   1  -1 -2  2
   !  -1   1  3 -2
   !b:
   !   -4
   !    1
   !   -3
   !    5
   !x:
   !    1
   !    2
   !    2
   !    1
   !行列の次元
   n   = 4
   !右辺の数
   nrhs= 1
   !行列の非零要素の数
   nnz = 16

   !CSR形式で行列を格納
   allocate(ia(n+1))
   allocate(ja(nnz))
   allocate( A(nnz))
   
   !ia:各行最初の非零要素の位置
   ia = (/ 1, 5, 9, 13, nnz+1/)
   !ja:非零要素の列インデックス
   ja = (/ 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4/)
   !A:行列の非零要素
   A  = (/2.d0,1.d0,-3.d0,-2.d0,2.d0,-1.d0,-1.d0,3.d0,1.d0,-1.d0,-2.d0,2.d0,-1.d0,1.d0,3.d0,-2.d0/)

   !右辺ベクトルを設定
   allocate(b(n,nrhs))
   b(:,1) = (/-4.d0,1.d0,-3.d0,5.d0/)

   !解ベクトルを設定
   allocate(x(n,nrhs))
   !-------------------------------------------------------

implicit_none, "Intel oneAPI Base Toolkit + HPC Toolkit (Intel Fortran)のインストール（Linux版）" (https://qiita.com/implicit_none/items/35bc4be8f2022903747a) ↩
Panua Technologyのページで提供されているFortranのサンプルコード内では，0の場合も1の場合も"use sparse direct solver"というコメントがついている．さらに，10を代入している例もある．しかしユーザーマニュアルでは0と1のみ記載されていたため，ここでもそのように紹介します． ↩
緒方隆盛, "疎行列直接法ソルバ入門", 計算工学, vol.25, No.2, 2020. ↩ ↩²
Hishimuma_t, "疎行列とベクトルを掛けたい貴方に", (https://zenn.dev/hishinuma_t/books/sparse-matrix-and-vector-product). ↩
Harmat, "「CSR形式を理解するのに時間がかかったので自分なりにまとめてみた」", (https://qiita.com/Harmat/items/793d65c3fe2fee2e7114). ↩
理数アラカルト, "連立一次方程式を掃き出し法で解く6つの例題" (https://risalc.info/src/simultaneous-linear-equations-examples.html#uni1) ↩