4 Harwell Boeing形式の読み込みサンプルコード2(複素エルミート行列、反復法)

ナビゲーション:前へ   上へ   次へ

4 Harwell Boeing形式の読み込みサンプルコード2(複素エルミート行列、反復法)

ここではHarwell Boeing形式の複素エルミート行列ファイル(MXTYPE=CHA)を読み込み、 列方向のインデックスを展開し並べ替えを行った後にNAGルーチンに渡して解くサンプルを示します。 ここで利用するNAGルーチンは F11JSF(複素エルミート反復法ソルバー) です。 F11JSF には列方向のインデックスが圧縮されていない形式で渡す必要がありますのでHarwell Boeing形式のデータを読み込んだ後に、 列方向のインデックスの展開(及び並べ替え)を行い、そして変換後のデータを用いて解を求めます。

MatrixJSF

4.1 読み込むデータ

読み込むデータは係数行列(Harwell Boeing形式の複素エルミート行列)と右辺ベクトル(テキスト形式)の2つです。

以下にHarwell Boeing形式で上記係数行列を記述した例を示します。 (青い記号は空白を表しています)

HarwellBoeingExampleImg

係数行列ファイルの例: 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 )



ナビゲーション:前へ   上へ   次へ

関連情報
Privacy Policy  /  Trademarks