Keyword: スパース, 不完全コレスキー分解
概要
本サンプルは実スパース対称行列の不完全コレスキー分解を行うFortranによるサンプルプログラムです。 本サンプルは実スパース対称行列Aを読み込み、不完全コレスキー分解を計算し、行列Aと以下に示される下三角行列Cの非ゼロ要素を出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン f11jaf() のExampleコードです。本サンプル及びルーチンの詳細情報は f11jaf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はf11jaf のマニュアルページを参照）| このデータをダウンロード |
F11JAF Example Program Data 7 N 16 NNZ 0 0.0 LFILL, DTOL 'N' 0.0 MIC, DSCALE 'M' PSTRAT 4. 1 1 1. 2 1 5. 2 2 2. 3 3 2. 4 2 3. 4 4 -1. 5 1 1. 5 4 4. 5 5 1. 6 2 -2. 6 5 3. 6 6 2. 7 1 -1. 7 2 -2. 7 3 5. 7 7 A(I), IROW(I), ICOL(I), I=1,...,NNZ
- １行目はタイトル行で読み飛ばされます。
- ２行目に行列Aの次数(n)を指定しています。
- ３行目に行列Aの下三角部分の非ゼロ要素の数(nnz)を指定しています。
- ４行目の１番目のパラメータ(lfill)は値が0以上の場合、コレスキー分解で許されたフィル（fill）の最大レベルを指定しています。２番目のパラメータ(dtol)は、１番目のパラメータの値が負の値の場合、フィルイン（fill-in）を制御するためのドロップ許容値を指定しています。１番目のパラメータの値が0以上の場合は２番目のパラメータの値は参照されません。
- ５行目には行の合計を保持するためにコレスキー分解が修正されるかどうか（mic='N':コレスキー分解が修正されない）と対角スケーリング引数（dscale=0.0）を指定しています。
- ６行目にピボット法（pstrat='M':Markowitz法を用いて対角ピボットが実行される）を指定しています。
- ７〜２２行目に行列Aの非ゼロ要素、行列Aの非ゼロ要素の行インデックス(irow)、列インデックス(icol)を指定しています。
出力結果
（本ルーチンの詳細はf11jaf のマニュアルページを参照）| この出力例をダウンロード |
F11JAF Example Program Results
Original Matrix
N = 7
NNZ = 16
1 0.4000E+01 1 1
2 0.1000E+01 2 1
3 0.5000E+01 2 2
4 0.2000E+01 3 3
5 0.2000E+01 4 2
6 0.3000E+01 4 4
7 -0.1000E+01 5 1
8 0.1000E+01 5 4
9 0.4000E+01 5 5
10 0.1000E+01 6 2
11 -0.2000E+01 6 5
12 0.3000E+01 6 6
13 0.2000E+01 7 1
14 -0.1000E+01 7 2
15 -0.2000E+01 7 3
16 0.5000E+01 7 7
Factorization
N = 7
NNZ = 16
NPIVM = 0
17 0.5000E+00 1 1
18 0.3333E+00 2 2
19 0.3333E+00 3 2
20 0.2727E+00 3 3
21 -0.5455E+00 4 3
22 0.5238E+00 4 4
23 -0.2727E+00 5 3
24 0.2683E+00 5 5
25 0.6667E+00 6 2
26 0.5238E+00 6 4
27 0.2683E+00 6 5
28 0.3479E+00 6 6
29 -0.1000E+01 7 1
30 0.5366E+00 7 5
31 -0.5345E+00 7 6
32 0.9046E+00 7 7
I IPIV(I)
1 3
2 4
3 5
4 6
5 1
6 2
7 7
- ２〜２０行目には元の行列の情報が出力されています。
- ４行目に行列Aの次数が出力されています。
- ５行目に行列Aの非ゼロ要素の数が出力されています。
- ７〜２２行目に行列Aの要素の値、行インデックスと列インデックスが出力されています。
- ２２〜４１行目にはコレスキー分解して得られた行列の情報が出力されています。
- ２３行目にコレスキー分解して得られた行列の次数が出力されています。
- ２４行目にコレスキー分解して得られた行列の非ゼロ要素の数が出力されています。
- ２５行目にコレスキー分解で修正されたピボットの数が出力されています。
- ４３〜５０行目にピボットインデックスが出力されています。
ソースコード
（本ルーチンの詳細はf11jaf のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
PROGRAM f11jafe
! F11JAF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
! .. Use Statements ..
USE nag_library, ONLY : f11jaf, nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6
! .. Local Scalars ..
REAL (KIND=nag_wp) :: dscale, dtol
INTEGER :: i, ifail, la, lfill, liwork, n, nnz, &
nnzc, npivm
CHARACTER (1) :: mic, pstrat
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: a(:)
INTEGER, ALLOCATABLE :: icol(:), ipiv(:), irow(:), istr(:), &
iwork(:)
! .. Executable Statements ..
WRITE (nout,*) 'F11JAF Example Program Results'
! Skip heading in data file
READ (nin,*)
! Read algorithmic parameters
READ (nin,*) n
READ (nin,*) nnz
la = 2*nnz
liwork = 2*la + 7*n + 1
ALLOCATE (a(la),icol(la),ipiv(n),irow(la),istr(n+1),iwork(liwork))
READ (nin,*) lfill, dtol
READ (nin,*) mic, dscale
READ (nin,*) pstrat
! Read the matrix A
DO i = 1, nnz
READ (nin,*) a(i), irow(i), icol(i)
END DO
! Calculate incomplete Cholesky factorization
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
CALL f11jaf(n,nnz,a,la,irow,icol,lfill,dtol,mic,dscale,pstrat,ipiv, &
istr,nnzc,npivm,iwork,liwork,ifail)
! Output original matrix
WRITE (nout,*) ' Original Matrix'
WRITE (nout,99997) 'N =', n
WRITE (nout,99997) 'NNZ =', nnz
DO i = 1, nnz
WRITE (nout,99999) i, a(i), irow(i), icol(i)
END DO
WRITE (nout,*)
! Output details of the factorization
WRITE (nout,*) ' Factorization'
WRITE (nout,99997) 'N =', n
WRITE (nout,99997) 'NNZ =', nnzc
WRITE (nout,99997) 'NPIVM =', npivm
DO i = nnz + 1, nnz + nnzc
WRITE (nout,99999) i, a(i), irow(i), icol(i)
END DO
WRITE (nout,*)
WRITE (nout,*) ' I IPIV(I)'
DO i = 1, n
WRITE (nout,99998) i, ipiv(i)
END DO
99999 FORMAT (1X,I8,E16.4,2I8)
99998 FORMAT (1X,2I8)
99997 FORMAT (1X,A,I16)
END PROGRAM f11jafe