4 Harwell Boeing形式の読み込みサンプルコード２（複素エルミート行列、反復法）
ここではHarwell Boeing形式の複素エルミート行列ファイル（MXTYPE=CHA）を読み込み、 列方向のインデックスを展開し並べ替えを行った後にnAGルーチンに渡して解くサンプルを示します。 ここで利用するnAGルーチンは F11JSF（複素エルミート反復法ソルバー） です。 F11JSF には列方向のインデックスが圧縮されていない形式で渡す必要がありますのでHarwell Boeing形式のデータを読み込んだ後に、 列方向のインデックスの展開（及び並べ替え）を行い、そして変換後のデータを用いて解を求めます。
4.1 読み込むデータ
読み込むデータは係数行列（Harwell Boeing形式の複素エルミート行列）と右辺ベクトル（テキスト形式）の２つです。
以下にHarwell Boeing形式で上記係数行列を記述した例を示します。 （青い記号は空白を表しています）
係数行列ファイルの例： f11jsf_matrix.hb
右辺ベクトルはN行（N=行列のオーダー。上記例では5）のテキストデータで、 それぞれの行に一つの数値を持つ形式のものです。 上記例の場合以下のようになります。
(8,54) (-10,-92) (25,27) (26,-28) (54,12) (26,-22) (47,65) (71,-57) (60,70)
右辺ベクトルファイルの例: f11jsf_uhen.txt
4.2 サンプルコード
以下にサンプルコードを示します。 以下のサンプルはHarwell Boeing形式の複素エルミート行列を読み込む必要最小限の処理を行っています。 （※実行列、複素非エルミート行列の読み込み等Harwell Boeing形式のすべてに対応してはいませんのでご注意ください）
※以下のサンプルはnAG数値計算ライブラリ（トライアルのご案内 ）を利用します。
PROGRAM solve_hermitian
IMPLICIT NONE
CHARACTER filename*256, mxtype*3, ptrfmt*16, indfmt*16, valfmt*20
CHARACTER method*6, precon
INTEGER i, j, k, ifail, nnz, n, nrows, ncols, lwork, maxitn, itn
DOUBLE PRECISION rnorm, omega
COMPLEX (kind(0D0)), ALLOCATABLE :: values(:), rhs(:), x(:), work(:)
DOUBLE PRECISION, ALLOCATABLE :: rdiag(:)
INTEGER, ALLOCATABLE :: iwork_f11zpf(:), iwork(:), istr(:)
INTEGER, ALLOCATABLE :: row_idx(:), compressed_col_idx(:), col_idx(:)
! ---------------------------------------------------------
! Read Harwell Boeing Sparse Matrix File
! ---------------------------------------------------------
PRINT *, 'Enter Harwell Boeing Matrix File Name'
READ *, filename
OPEN (10, FILE=filename, STATUS='old')
READ (10, '(A72,A8)') ! title, key
READ (10, '()') ! Skip 5 I14 format: totcrd, ptrcrd, indcrd, valcrd, rhscrd
READ (10, '(A3,11X,4I14)') mxtype, nrows, ncols, nnz ! ,neltvl
PRINT '(I0, " x ", I0, " (",I0," non-zeros) mxtype is ",A)', nrows, ncols, nnz, mxtype
IF (mxtype(1:1)/='C' .AND. mxtype(1:1)/='c') STOP 'Only complex matrix is accepted...'
IF (mxtype(2:2)/='H' .AND. mxtype(2:2)/='h') STOP 'Only hermitian matrix is accepted...'
n = ncols ! n is the order of the matrix
! Allocate and read coordinates and data (nnz = number of non-zeros)
ALLOCATE (compressed_col_idx(n+1), row_idx(nnz), values(nnz))
READ (10, '(2A16,2A20)') ptrfmt, indfmt, valfmt ! ,rhsfmt
READ (10, ptrfmt)(compressed_col_idx(i), i=1, n+1) ! read compressed column indecies
READ (10, indfmt)(row_idx(i), i=1, nnz) ! read row indecies
READ (10, valfmt)(values(i), i=1, nnz) ! read values
CLOSE (10)
! ---------------------------------------------------------
! Allocate and read RHS
! ---------------------------------------------------------
! open file
PRINT *, 'Enter Vector (txt) File Name'
READ *, filename
ALLOCATE (rhs(n))
OPEN (10, FILE=filename, STATUS='old')
DO i = 1, n
READ (10, *) rhs(i)
END DO
CLOSE (10)
! ---------------------------------------------------------
! Convert CCS to CS(SCS)
! ---------------------------------------------------------
PRINT *, 'Converting to CS(SCS)...'
ALLOCATE (col_idx(nnz))
j = 1
DO i = 1, ncols
k = compressed_col_idx(i+1) - compressed_col_idx(i)
IF (k>0) THEN
col_idx(j:j+k-1) = i
j = j + k
END IF
END DO
! ---------------------------------------------------------
! Sort CS(SCS)
! ---------------------------------------------------------
PRINT *, 'Sorting CS(SCS)...'
ALLOCATE (iwork_f11zpf(n), istr(n+1))
ifail = 0
CALL f11zpf(n, nnz, values, row_idx, col_idx, 'R', 'K', istr, iwork_f11zpf, ifail)
! ---------------------------------------------------------
! Solve Ax = b using f11jsf
! ---------------------------------------------------------
PRINT *, 'Solving...'
lwork = 7*n + 120
ALLOCATE (x(n), rdiag(n), work(lwork), iwork(n+1))
ifail = 0
omega = 1.0D0 ! Omega parameter: only for precon='S', 0.0<=omega<=2.0
maxitn = 30000 ! Maximum number of iteration
precon = 'N' ! Preconditioner: 'N'=none, 'J'=Jacobi, 'S'=SSOR
method = 'CG' ! Method: 'CG'=Conjugate Gradient, 'SYMMLQ'=Lanczos
x = 0D0 ! Set initial approximation to 0
CALL f11jsf(method, precon, n, nnz, values, row_idx, col_idx, omega, rhs, 0D0, maxitn, x, rnorm, &
itn, rdiag, work, lwork, iwork, ifail)
! ---------------------------------------------------------
! Write results
! ---------------------------------------------------------
PRINT *, 'Results: number of iterations =', itn
PRINT '("(",f8.5,", ",f8.5," )")', x
END PROGRAM solve_hermitian
サンプルコード: solve_complex_hermitian.f90
注意点等：
※本サンプルは複素エルミート行列（Harwell Boeing MXTYPE="CHA"）のみに対応しています。
※本サンプルで用いるnAGルーチンはSCS形式（圧縮されていないカラム形式）に対応したルーチンですので、
Harwell Boeing形式から読み込まれたcompressed_col_idx（圧縮された列インデックス）はまず圧縮されていない形に変換しています。
（"Convert CCS to CS(SCS)"の部分）
そしてその後にnAGルーチンが必要とする順番に並べ替えています。（"Sort CS(SCS)"の部分）
※引数詳細等につきましては各ルーチン（ f11jsf , f11zpf ）のマニュアルページをご参照下さい。
4.3 実行結果例
以下に本サンプルプログラムの実行結果例を示します。
Enter Harwell Boeing Matrix File Name f11jsf_matrix.hb 9 x 9 (23 non-zeros) mxtype is CHA Enter Vector (txt) File Name f11jsf_uhen.txt Converting to CS(SCS)... Sorting CS(SCS)... Solving... Results: number of iterations = 9 ( 1.00000, 9.00000 ) ( 2.00000, -8.00000 ) ( 3.00000, 7.00000 ) ( 4.00000, -6.00000 ) ( 5.00000, 5.00000 ) ( 6.00000, -4.00000 ) ( 7.00000, 3.00000 ) ( 8.00000, -2.00000 ) ( 9.00000, 1.00000 )
ナビゲーション：前へ 上へ 次へ