3 Harwell Boeing形式の読み込みサンプルコード1(実非対称行列、直接法)

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

3 Harwell Boeing形式の読み込みサンプルコード1(実非対称行列、直接法)

ここではHarwell Boeing形式の実非対称行列ファイル(MXTYPE=RUA)を読み込み、 列方向のインデックスが圧縮されたままの形式(CCS形式)でNAGルーチンに渡して解くサンプルを示します。 前述のソルバールーチンである F11MFF(実非対称直接法ソルバールーチン)は Harwel-Boeing形式と同様の CCS形式に直接対応しています。ここで紹介するサンプルは以下に示される問題を、 標準的なHarwell Boeing形式の実非対称行列とテキスト形式で記述された右辺ベクトルを読み込んで解きます。

MatrixImage

3.1 読み込むデータ

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

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

MatrixImageInHarwellBoeing

係数行列ファイルの例: f11mff_matrix.hb

右辺ベクトルはN行(N=行列のオーダー。上記例では5)のテキストデータで、 それぞれの行に一つの数値を持つ形式のものです。 上記例の場合以下のようになります。

1.56
-0.25
3.60
1.33
0.52

右辺ベクトルファイルの例: f11mff_uhen.txt

3.1.1 サンプルコード

以下にサンプルコードを示します。 以下のサンプルはHarwell Boeing形式の実非対称行列を読み込む必要最小限の処理を行っています。 (※複素行列の読み込み等Harwell Boeing形式のすべてに対応してはいませんのでご注意ください)

※以下のサンプルはNAG数値計算ライブラリ(トライアルのご案内 )を利用します。

PROGRAM solve_real_nonsymm
    IMPLICIT NONE
    CHARACTER filename*256, mxtype*3, ptrfmt*16, indfmt*16, valfmt*20
    INTEGER i, ifail, nnz, n, nrows, ncols, nzlmx, nzlumx, nzumx, nnzu, nnzl, factor
    DOUBLE PRECISION flop
    DOUBLE PRECISION, ALLOCATABLE :: values(:), x(:), uval(:), lval(:)
    INTEGER, ALLOCATABLE :: iu(:), il(:), row_idx(:), compressed_col_idx(:), iprm(:)

    ! ---------------------------------------------------------
    ! 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)/='R' .AND. mxtype(1:1)/='r') STOP 'Only real matrix is accepted...'
    IF (mxtype(2:2)/='U' .AND. mxtype(2:2)/='u') STOP 'Only unsymmetric matrix is accepted...'
    IF (nrows/=ncols) STOP 'Only square matrix is accepted...'
    n = ncols  ! n is the order of the matrix
    ! Allocate and read coordinates and data (nnz = number of non-zero's)
    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
    ! ---------------------------------------------------------
    PRINT *, 'Enter Vector (txt) File Name'
    READ *, filename
    ALLOCATE (x(n))
    OPEN (10, FILE=filename, STATUS='old')
    DO i = 1, n
        READ (10, *) x(i)
    END DO
    CLOSE (10)

    ! ---------------------------------------------------------
    ! Calculate column permutation
    ! ---------------------------------------------------------
    PRINT *, 'Calculate column permutation...'
    ifail = 0
    ALLOCATE (iprm(n*7))
    CALL f11mdf('M', n, compressed_col_idx, row_idx, iprm, ifail)

    ! ---------------------------------------------------------
    ! Factorize
    ! ---------------------------------------------------------
    PRINT *, 'Factorize...'
    ifail = -1
    factor = 32             ! problem dependent 'factor', change this if ifail==4
    nzlmx = factor*nnz
    nzlumx = factor*nnz
    nzumx = factor*nnz
    ALLOCATE (uval(nzumx), lval(nzlumx), iu(2*n+nzumx+1), il(7*n+nzlmx+4))
    CALL f11mef(n, row_idx, values, iprm, 1D0, nzlmx, nzlumx, nzumx, il, lval, iu, uval, nnzl, nnzu, flop, ifail)
    if ( ifail==4 ) print *, "Set 'factor' >=", nzlumx/nnz+1
    if ( ifail/=0 ) stop

    ! ---------------------------------------------------------
    ! Compute solution using f11mff
    ! ---------------------------------------------------------
    PRINT *, 'Computing solution...'
    ifail = 0
    CALL f11mff('N', n, iprm, il, lval, iu, uval, 1, x, n, ifail)

    ! ---------------------------------------------------------
    ! Write results
    ! ---------------------------------------------------------
    PRINT *, "Result:"
    PRINT '(g16.10)', x

END PROGRAM solve_real_nonsymm

サンプルコード: solve_real_nonsymm.f90

注意点等:
※今回利用したNAGルーチンはCCS形式(圧縮カラム形式)に対応したルーチンですので、 Harwell Boeing形式から読み込まれたrow_idx(行インデックス)、compressed_col_idx(圧縮された列インデックス)、values(値)は 特に変換等行われずにそのまま引数として渡されています。

※ソルバールーチンであるf11mffはその呼び出し前に2つの前処理ルーチン(f11mdf, f11mef)を必要としますので、 データの読み込みが完了後に f11mdf, f11mef, f11mff の順に呼び出しが行われています。

※nzlmx, nzlumx, nzumxの設定値(現在factor*nnz,factor=32)は行列の分解の際に必要となる領域のサイズを決定するものですが、 解こうとする問題によってはより大きな値の設定が必要となる場合もあります。詳細はf11mefのマニュアルページをご参照下さい。

※引数詳細等につきましては各ルーチン( f11mdf , f11mef , f11mff )のマニュアルページをご参照下さい。

3.2 実行結果例

以下に本サンプルプログラムの実行結果例を示します。

 Enter Harwell Boeing Matrix File Name
f11mff_matrix.hb
5 x 5 (11 non-zeros)  mxtype is RUA
 Enter Vector (txt) File Name
f11mff_uhen.txt
 Calculate column permutation...
 Factorize...
 Computing solution...
 Result:
0.7000000000
0.1600000000
0.5200000000
0.7700000000
0.2800000000



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

関連情報
Privacy Policy  /  Trademarks