前提・実現したいこと
ここに質問の内容を詳しく書いてください。
FORTRANで行列のLU分解のコードを数値計算の理論と実際(河村哲也)という本を見ながら写経しました。次の3×3行列で出力結果が次のようになりました。(整数で書きましたが実数で入力してます)
[1 2 3] [1 0 0] [1 2 3]
[2 3 4] = [2 1 0] [2-1-2]
[3 4 5] [3 2 1] [0 0 0]
上三角行列の形になってません。
おかしな点についてわかる方がおられましたら誠にお手数ですがご回答していただけるとありがたいです。
発生している問題・エラーメッセージ
3×3行列などを設定した場合、LU分解の結果で、上三角行列でない行列が計算されてしまう。
該当のソースコード
C********************
C LU分解
C********************
IMPLICIT REAL8(A-H,L,O-Z)
PARAMETER(M=10)
DIMENSION A(M,M),L(M,M),U(M,M)
WRITE(,)'行列の大きさM((M,M)行列)を入力してください'
READ(,) N
DO 15 I=1,N
DO 10 J=1,N
WRITE(,) '行列の要素を入力してください.A(',I,J,')?'
READ(,) A(I,J)
10 CONTINUE
15 CONTINUE
PVT = 1.D-5
WRITE(,) '行列の要素'
WRITE(,600) ((A(I,J),J=1,N),I=1,N)
CALL LUDECM(A,M,N,PVT,L,U)
WRITE(,) '---------------計算結果---------------'
WRITE(,) '下三角行列の要素'
WRITE(,600) ((L(I,J),J=1,N),I=1,N)
WRITE(,) '上三角行列の要素'
WRITE(,600) ((U(I,J),J=1,N),I=1,N)
600 FORMAT(1H ,4F10.5)
STOP
END
C
C*****LU分解のサブルーチン
C
SUBROUTINE LUDECM(A,M,N,PVT,L,U)
IMPLICIT REAL8(A-H,L,O-Z)
DIMENSION A(M,M),L(M,M),U(M,M)
DO 10 J=1,N
DO 10 I=1,N
L(I,J) = 0.0
U(I,J) = 0.0
10 CONTINUE
C
DO 20 I = 1,N-1
IF(DABS(A(I,I)).LT.PVT) THEN
WRITE(,*)I,'ステップでピボットが小さくなりすぎました'
END IF
DIAG = 1.D0/A(I,I)
DO 30 K=I+1,N
AM = A(K,I)DIAG
L(K,I) = AM
DO 40 J=I+1,N
A(K,J) = A(K,J) - AMA(I,J)
40 CONTINUE
30 CONTINUE
DO 50 J=1,N
U(I,J) = A(I,J)
50 CONTINUE
20 CONTINUE
C
U(N,N) = A(N,N)
DO 60 I=1,N
L(I,I) = 1.0
60 CONTINUE
RETURN
END
FORTRAN77
試したこと
テキストとコードの比較
補足情報(FW/ツールのバージョンなど)
ここに貼り付けた際、先頭の6カラムはなくなってしまいました。
コンパイルはできています。