fortranで行列のLU分解のコードを写経したのですがLU分解ができていません

Question

前提・実現したいこと

ここに質問の内容を詳しく書いてください。
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カラムはなくなってしまいました。
コンパイルはできています。

Comments

[seastar3]

2X2,3X3以降の読み込むデータ例を示してもらえませんか。

[shuim]

例えば2×2→1,2,2,3 3×3→1,2,3,2,3,4,3,4,5などです。(1行の成分、2行の成分…の順)

[shuim]

著者名　河村哲也の誤りでした。

[seastar3]

いや、正確に１行目コマ数、以下個々のデータの形式で質問文に示してほしいのです。,はREAD文に対応していません。

[seastar3]

著者名も編集アイコンを押して書き直しましょう。編集理由を書かないと反映されませんので気をつけて。

[shuim]

著者名を修正いたしました。申し訳ありませんが質問文に示してほしいといわれたことの意味が分かりません。コンパイルはできると思いますし、exeを実行したら数字入力、Enterの繰り返しで再現できると思います。

[seastar3]

キー入力の前提なのですね。昔のパンチカード方式を思い浮かべていたのでテキストファイルのフォーマットがずれている可能性を考えていました。それはそれでFortranコードを```Fortran77 ・・・ ```で挟みMarkdown記法に合わせると表記が見やすくなりますよ。

Answer 1

普通、FORTRAN77 時代のルーチンでは結果として返される以外の配列部分は、ワーク領域に使ったりしているので、上三角以外の部分にゴミが残っているのだろう。

実際、下三角と上三角を掛ければ正しい結果が出ている。

もう一度よく本を読み直すと、そういう但し書きがあるのではないかと思う。

Answer 2

また、24行目の

600 FORMAT(1H ,4F10.5)

の4F10.5の実数4項のみ表示という書式が4項を越えると溢れるのではないかと考えます。

600 FORMAT(1H ,10F10.5)

なり増やす必要があるのではないでしょうか。

Answer 3

行列の変換を詳しくトレースしていませんが、簡単なデバッグ作業として、LUDECMサブルーチン内のループ内に、

WRITE(*,700) I,J
WRITE(*,710) ((U(a,b),b=1,N),a=1,N)
710 FORMAT(2I5)
710 FORMAT(1H ,10F10.5)

のようなダンプ用のコードを入れて確認していってはどうでしょうか。