PETSc - 超音波流体屋のプログラム備忘録

PETScとはC/Fortranで使える、反復解法スパースソルバで有力なオープンソースコード。
PETScの入門的な使い方として、連立方程式の求解の計算例を示す。
https://www.mcs.anl.gov/petsc/index.html

インストール

https://www.mcs.anl.gov/petsc/documentation/installation.html
make installがないので、使用予定の場所に解凍する。　(2016/10/20編集)
OpenMPI環境で$MPI_ROOTにインストールされている場合下記が多分よい。
(追記)blas,lapackがyumで手に入っているなら--downloadの部分も不要。
その際libblas.soおよびliblapack.soにアクセスできるよう、適切なシンボリックリンクを貼ること。

./configure --with-mpi-dir=$MPI_ROOT --download-fblaslapack --with-scalar-type=complex --prefix=/usr/local/petsc
make all test
sudo make install

解凍後に表示される環境変数などを.bashrcなどに登録しておくとよい

export PETSC_DIR=/usr/local/petsc
export PETSC_ADIR=$PETSC_DIR
export PETSC_INC=$PETSC_DIR/include
export PETSC_AINC=$PETSC_ADIR/include
export PETSC_LIB=$PETSC_ADIR/lib
export LD_LIBRARY_PATH=$PETSC_LIB:$LD_LIBRARY_PATH

コンパイルなど

ifortがFPATHを見てくれみたいなので、-Iオプションでインクルードパスを指定

ifort hoge.f90 -o hoge.out -fpp -lpetsc -I${PETSC_INC} -I${PETSC_AINC}

ベクトルの入力

ベクトル・行列ともにアクセス方法は三種類で、
VecSetValueで個別に指定するか、VecSetValuesでインデックスと値とともに指定するか、
ポインタを取得・再格納するかに限られると思う。
PetscScalarは上記設定でインストールした場合complex(8)として扱える。
行列の格納はほとんどベクトルと似た関数が用意されている。

注意点としてVecSetValue(s)は引数に直接1d0などのreal(8)型を指定すると、
周囲の要素を壊してしまう傾向にあるため、一旦PetScScalarに格納するのが安全である。

program petsc_test
implicit none
#include <petsc/finclude/petscsys.h>
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscvec.h90>
 
    integer, parameter :: n = 4
    integer :: ierr, idx(n)
    Vec :: v
    PetscScalar, pointer :: p_v(:)
    PetscScalar :: val(n), x
 
    val(:) = 1d0; val(3) = 4d0;
    idx = (/0,1,2,3/) !0 start indexに注意
 
    call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
    call VecCreateSeq(PETSC_COMM_SELF, n, v,ierr)
 
    call VecSetValues(v, n, idx, val, INSERT_VALUES, ierr)
    x = 3d0 !手間だが一旦格納する
    call VecSetValue(v, 1, x, INSERT_VALUES,ierr) !idx:1 (0 startに注意)
 
    call VecGetArrayF90(v,p_v,ierr)
    p_v(1) = (2d0, -4d0)
    call VecRestoreArrayF90(v,p_v,ierr)
 
    call VecView(v, PETSC_VIEWER_STDOUT_SELF, ierr)
 
    call VecDestroy(v,ierr)
    call PetscFinalize(ierr)
 
end program petsc_test
 

結果:

Vec Object: 1 MPI processes
  type: seq
2. - 4. i
3.
4.
1.

行列の入力・連立方程式の求解

スパース行列を扱う上ではポインタ取得をするよりかは、MatSetValuesを使用するとよいと思う。
データ管理はPETScサイドでやってくれるので、COOやらCSRやら悩まないで済む。
行列は入力を終えたらMatAssemblyBegin/Endをしないとエラーになる。
反復/直接解法はKSP(Krylov SubsPace method)の構造体から行う。
直接解法は前処理からLUやCHOLESKYを選択すれば実施できる。

program petsc_mat
implicit none
 
#include <petsc/finclude/petscsys.h>
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscvec.h90>
#include <petsc/finclude/petscmat.h>
#include <petsc/finclude/petscksp.h>
#include <petsc/finclude/petscpc.h>
 
    integer, parameter :: n = 3
    integer :: ierr, idx(n), iter
    PetscScalar :: v(n), zero
    Mat :: a
    Vec :: b, x
 
    KSP :: ksp !連立方程式ソルバ構造体
    PC  :: pc  !前処理構造体
    KSPConvergedReason :: res
    integer(8) :: clock1, clock2, crate, cmax !時間計測用
 
    call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Mat prep.
    call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,n,PETSC_NULL_INTEGER,a,ierr)
    v(:) = (/8d0, 1d0, 6d0/)
    call MatSetValues(a, 1, 0, n, (/0,1,2/), v, INSERT_VALUES,ierr)
    v(:) = (/3d0, 5d0, 7d0/)
    call MatSetValues(a, 1, 1, n, (/0,1,2/), v, INSERT_VALUES,ierr)
    v(:) = (/4d0, 9d0, 2d0/)
    call MatSetValues(a, 1, 2, n, (/0,1,2/), v, INSERT_VALUES,ierr)
    call MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY,ierr)
    call MatAssemblyEnd(a,MAT_FINAL_ASSEMBLY,ierr)
    call MatView(a, PETSC_VIEWER_STDOUT_SELF, ierr)
 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Vec prep.
    call VecCreateSeq(PETSC_COMM_SELF, n, b,ierr)
    v(:) = (/1d0, 1d0, 1d0/)
    call VecSetValues(b, n, (/0,1,2/), v, INSERT_VALUES,ierr)
    call VecView(b, PETSC_VIEWER_STDOUT_SELF, ierr)
 
    call VecDuplicate(b,x,ierr)
    zero = 0d0 !1d0/15d0 !で初期値ありにするとイタレーションが0になる
    call VecSet(x,zero,ierr) !あらかじめ初期値を指定できる
 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !KSP prep. (Krylov SubsPace method)
    call KSPCreate(PETSC_COMM_SELF,ksp,ierr)
    call KSPSetType(ksp, KSPBICG, ierr) !何も指定しないとGMRESが入る
    !KSPTYPEはKSPのあとに, CG, BICG, MINRES, GMRES, BCGS(=biCgStab)
 
    call KSPSetOperators(ksp,a,a,ierr) !前のaは解く行列,後ろは前処理用(同じでよい)
    call KSPSetInitialGuessNonzero(ksp,PETSC_TRUE,ierr) !指定しないとゼロ初期値で解く
 
    call KSPGetPC(ksp,pc,ierr) 
    call PCSetType(pc,PCJACOBI,ierr) 
    !PCTYPEはPCのあとにJACOBI, ILU, ICC(不完全cholesky), LU, CHOLESKY(直接法)
 
    call KSPSetTolerances(ksp,1d-9,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
    call KSPSetFromOptions(ksp,ierr)
    call KSPView(ksp,PETSC_VIEWER_STDOUT_SELF,ierr) !オプションのレビュー
 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Solve
    call system_clock(clock1,crate,cmax)
    call KSPSolve(ksp,b,x,ierr)
    call system_clock(clock2,crate,cmax)
 
    call KSPGetIterationNumber(ksp,iter,ierr)
    call KSPGetConvergedReason(ksp, res, ierr) !正なら収束, 負なら発散
    ! 1,2:rtol, 3,9:atol, -3:its, -4:dtol, -5,6:breakdown, -9:nan/inf
    ! -7:non-sym, -8,10:indefinite, -11:pcfail
    print *, 'Status:', res, 'Iterations:', iter, ', time [s]', real(clock2 -clock1)/real(crate)
 
    call VecView(x, PETSC_VIEWER_STDOUT_SELF, ierr)
 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Finalize
    call KSPDestroy(ksp,ierr)
    call MatDestroy(a,ierr)
    call VecDestroy(b,ierr)
    call PetscFinalize(ierr)
 
end program petsc_mat
 

結果

Mat Object: 1 MPI processes
  type: seqaij
row 0: (0, 8.)  (1, 1.)  (2, 6.) 
row 1: (0, 3.)  (1, 5.)  (2, 7.) 
row 2: (0, 4.)  (1, 9.)  (2, 2.) 
Vec Object: 1 MPI processes
  type: seq
1.
1.
1.
KSP Object: 1 MPI processes
  type: bicg
  maximum iterations=10000
  tolerances:  relative=1e-09, absolute=1e-50, divergence=10000.
  left preconditioning
  using nonzero initial guess
  using DEFAULT norm type for convergence test
PC Object: 1 MPI processes
  type: jacobi
  PC has not been set up so information may be incomplete
  linear system matrix = precond matrix:
  Mat Object:   1 MPI processes
    type: seqaij
    rows=3, cols=3
    total: nonzeros=9, allocated nonzeros=9
    total number of mallocs used during MatSetValues calls =0
      using I-node routines: found 1 nodes, limit used is 5
 Status:           2 Iterations:           3 , time [s]  4.4999999E-04
Vec Object: 1 MPI processes
  type: seq
0.0666667
0.0666667
0.0666667