質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.61%

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

受付中

回答 3

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 1,194

shuim

score 8

 前提・実現したいこと

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

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • seastar3

    2018/10/20 21:58

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

    キャンセル

  • shuim

    2018/10/21 00:11

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

    キャンセル

  • seastar3

    2018/10/21 04:04

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

    キャンセル

回答 3

0

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

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


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

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

0

また、24行目の

600 FORMAT(1H ,4F10.5)


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

600 FORMAT(1H ,10F10.5)


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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

0

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

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

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.61%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る