Keyword: スパースソルバー, sparse solver, RGMRES, 疎行列, SSOR
概要
本サンプルは以下に示す非対称実スパース連立方程式を解くFortranによるプログラムです。 連立方程式は係数行列と右辺ベクトルとして与えます。
以下のプログラム例ではRGMRES法 (restarted generalized minimal residual method) を用いていますが、その他にCGS法 (conjugate gradient squared method)、及びBICGSTAB法 (bi-conjugate gradient stabilised method) に対応しています。（入力データの'RGMRES'の部分を'CGS'、'BiCGSTAB'にそれぞれ置き換える事で手法を選ぶことが可能です）
またプリコンディショナーにはSSOR (symmetric successive-over-relaxation) を用いていますが、その他にJacobi の指定及びプリコンディショナーを使用しない設定が可能です。（入力データの'S'の部分をそれぞれ'J'、'N'に置き換えることで指定可能です）
※本サンプルはnAG Fortranライブラリに含まれるルーチン f11def()のExampleコードです。本サンプル及びルーチンの詳細情報はf11def のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はf11def のマニュアルページを参照）| このデータをダウンロード |
F11DEF Example Program Data 5 1 N, M 16 NNZ 'RGMRES' 'S' METHOD, PRECON 1.05 OMEGA 1.D-10 1000 TOL, MAXITN 2. 1 1 1. 1 2 -1. 1 4 -3. 2 2 -2. 2 3 1. 2 5 1. 3 1 5. 3 3 3. 3 4 1. 3 5 -2. 4 1 -3. 4 4 -1. 4 5 4. 5 2 -2. 5 3 -6. 5 5 A(I), IROW(I), ICOL(I), I=1,...,NNZ 0. -7. 33. -19. -28. B(I), I=1,...,N 0. 0. 0. 0. 0. X(I), I=1,...,N
- １行目はタイトル行で読み飛ばされます。
- ２行目は係数行列の大きさ（n=5,5x5なので5)、RGMRES法の再開サブスペース空間の次元 m (=1) を指定しています。数字より後ろは行末まで読み飛ばされます。
- ３行目では係数行列（スパース）に含まれる非ゼロ要素の数(nnz=16)を与えます。数字より後ろは行末まで読み飛ばされます。
- ４行目は用いる計算手法 (method='RGMRES') とプリコンディショナー (precon='S'（SSOR）) を指定しています。 プリコンディショナーの指定より後ろは行末まで読み飛ばされます。
※ここで他の手法（'CGS' もしくは 'BiCGSTAB'）やプリコンディショナー （'J'（ヤコビ）もしくは 'N'（プリコンディショナー無し）） の指定も可能です。 - ５行目ではSSORプリコンディショナーのリラックスパラメータであるω (=1.05) を与えています。SSOR以外のプリコンディショナーを用いる場合にはこの数値は無視されます。 数字より後ろは行末まで読み飛ばされます。
- ６行目では計算精度 tol (=1.D-10) 、最大反復数 maxitn (=1000) をそれぞれ与えています。 ２つの数値より後ろは行末まで読み飛ばされます。
- ７行目〜２２行目にはスパース行列の非ゼロ係数の値を指定しています。それぞれの行は値(a)、行番号(irow)、列番号(icol)で全部で16個（３行目で与えた値）の非ゼロ係数が指定されています。行番号及び列番号は1から最大5（２行目で指定した値）までの値を取ります。 例えば行列の一番左の一番上の値が3.5であった場合には
3.5 1 1
と与えます。 - ２３行目と２４行目は右辺ベクトルの値(b)を指定しています。全部で５個（２行目で指定した値）の数値が指定されています。
- ２５行目と２６行目は解の初期値(x)を与えています。 全部で５個（２行目で指定した値）の数値が指定されています。 ここでは初期値としてすべて 0 を与えています。
出力結果
（本ルーチンの詳細はf11def のマニュアルページを参照）| この出力例をダウンロード |
F11DEF Example Program Results
Converged in 13 iterations
Final residual norm = 5.087E-09
X
1.0000E+00
2.0000E+00
3.0000E+00
4.0000E+00
5.0000E+00
- １行目はタイトルです。
- ３行目に解を求める際の反復数が示されます。（この例では１３回）
- ４行目は求まった解の残差ノルムが出力されます。
- ６行目は解のタイトルである x が出力されます。
- ７行目〜１１行目に求まった解が表示されます。
ソースコード
（本ルーチンの詳細はf11def のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
PROGRAM f11defe
! F11DEF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
! .. Use Statements ..
USE nag_library, ONLY : f11def, nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6
! .. Local Scalars ..
REAL (KIND=nag_wp) :: omega, rnorm, tol
INTEGER :: i, ifail, itn, l, lwork, m, maxitn, &
n, nnz
CHARACTER (8) :: method
CHARACTER (1) :: precon
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: a(:), b(:), work(:), x(:)
INTEGER, ALLOCATABLE :: icol(:), irow(:), iwork(:)
! .. Intrinsic Functions ..
INTRINSIC max
! .. Executable Statements ..
WRITE (nout,*) 'F11DEF Example Program Results'
WRITE (nout,*)
! Skip heading in data file
READ (nin,*)
! Read algorithmic parameters
READ (nin,*) n, m
READ (nin,*) nnz
READ (nin,*) method, precon
l = n
IF (precon=='N' .OR. precon=='n') l = 0
lwork = max(4*n+m*(m+n+5)+l+101,8*n+l+100,2*n*(m+3)+m*(m+2)+l+100, &
11*n+l+100)
ALLOCATE (a(nnz),b(n),work(lwork),x(n),icol(nnz),irow(nnz), &
iwork(2*n+1))
READ (nin,*) omega
READ (nin,*) tol, maxitn
! Read the matrix A
DO i = 1, nnz
READ (nin,*) a(i), irow(i), icol(i)
END DO
! Read right-hand side vector b and initial approximate solution x
READ (nin,*) b(1:n)
READ (nin,*) x(1:n)
! Solve Ax = b using F11DEF
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
CALL f11def(method,precon,n,nnz,a,irow,icol,omega,b,m,tol,maxitn,x, &
rnorm,itn,work,lwork,iwork,ifail)
WRITE (nout,'(A,I10,A)') ' Converged in', itn, ' iterations'
WRITE (nout,'(A,1P,E16.3)') ' Final residual norm =', rnorm
WRITE (nout,*)
! Output x
WRITE (nout,*) ' X'
WRITE (nout,'(1X,1P,E16.4)') x(1:n)
END PROGRAM f11defe