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

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

新規登録して質問してみよう
ただいま回答率
86.12%
Windows 10

Windows 10は、マイクロソフト社がリリースしたOSです。Modern UIを標準画面にした8.1から、10では再びデスクトップ主体に戻され、UIも変更されています。PCやスマホ、タブレットなど様々なデバイスに幅広く対応していることが特徴です。

FORTRAN

FORTRAN(フォートラン)は科学時術計算に向いた手続き型プログラミング言語です。 並列計算の最適化が行いやすい特性上、数値予報および気候モデルなどの大規模な計算を行う分野のスーパーコンピュータで使われています。

受付中

FORTRANによる逆行列の計算

danbo2000314
danbo2000314

総合スコア0

Windows 10

Windows 10は、マイクロソフト社がリリースしたOSです。Modern UIを標準画面にした8.1から、10では再びデスクトップ主体に戻され、UIも変更されています。PCやスマホ、タブレットなど様々なデバイスに幅広く対応していることが特徴です。

FORTRAN

FORTRAN(フォートラン)は科学時術計算に向いた手続き型プログラミング言語です。 並列計算の最適化が行いやすい特性上、数値予報および気候モデルなどの大規模な計算を行う分野のスーパーコンピュータで使われています。

1回答

0グッド

0クリップ

1405閲覧

投稿2021/11/13 12:26

前提・実現したいこと

多項式曲線フィッティングを行うために、wベクトルの計算が最終目標ですが、その過程でA_matの逆行列を求める必要がある。

発生している問題・エラーメッセージ

ガウス消去法の過程で、0割りが発生してしまっているのか、A_matの逆行列の値がすべて、NaNとなってしまっている。また、A_mat行列に関しても、一か所だけNaNとなっている個所がある。

エラーメッセージ

該当のソースコード

FORTRAN

1program main 2 3 implicit none 4 integer :: i, j 5 real(8) :: k, xn(10,1), tn(10,1), X_mat(10,10), A_mat(10,10), mid_mat(10,10), w(10,1), b 6 7 open(20, file='tn_data.txt') 8 9 do i = 1,10 10 b = real(i) 11 xn(i,1) = (b-1) / 9.0 12 end do 13 14 do i = 1,10 15 read(20,*) tn(i,1) 16 end do 17 18 close(20) 19 20 write(*,*) 'xn=' 21 do i = 1,10 22 write(*,*) xn(i,1) 23 end do 24 25 write(*,*) 'tn=' 26 do i = 1,10 27 write(*,*) tn(i,1) 28 end do 29 30 call X_calculation(xn, X_mat) 31 32 write(*,*) 'X_mat=' 33 do i = 1,10 34 write(*,*) (X_mat(i,j), j=1,10) 35 end do 36 37 call A_calculation(xn, A_mat) 38 39 write(*,*) 'A_mat=' 40 do i = 1,10 41 write(*,*) (A_mat(i,j), j=1,10) 42 end do 43 44 call matinv(A_mat, 10) 45 46 write(*,*) 'A_mat-1=' 47 do i = 1,10 48 write(*,*) (A_mat(i,j), j=1,10) 49 end do 50 51 mid_mat = matmul(A_mat,X_mat) 52 w = matmul(mid_mat,tn) 53 54 open(30,file='w_data.csv') 55 56 write(*,*) 'w=' 57 do i = 1,10 58 write(*,*) w(i,1) 59 write(30,*) w(i,1) 60 end do 61 62 close(30) 63 64stop 65end program main 66 67 68subroutine X_calculation(xn, X_mat) 69 70 implicit none 71 real(8) :: xn(10,1), X_mat(10,10), b 72 integer :: i, j 73 74 do i = 1,10 75 X_mat(1, i) = 1.0 76 end do 77 78 do i = 1,10 79 do j = 1,9 80 b = dble(j) 81 X_mat(j+1,i) = xn(i,1)**j 82 end do 83 end do 84 85end subroutine X_calculation 86 87 88subroutine A_calculation(xn,A_mat) 89 90 implicit none 91 real(8) :: xn(10,1), A_mat(10,10), sum, b, c 92 integer :: i, j, k 93 94 do i = 1,10 95 do j = 1,10 96 do k = 1,10 97 b = dble(i) 98 c = dble(j) 99 sum = xn(k,1)**(b+c-2.0) 100 A_mat(i,j) = A_mat(i,j) + sum 101 end do 102 end do 103 end do 104 105end subroutine A_calculation 106 107 108subroutine matinv(A_mat,n) 109 110real(8) :: A_mat(10,10),c,dum,anmax 111integer :: i,j,k,km,row(10) 112 113do i = 1, 10 114 115 anmax = 0.0 116 km = i 117 118 do k = i,10 119 if(abs(A_mat(k,i))>anmax)then 120 anmax = abs(A_mat(k,i)) 121 km = k 122 end if 123 end do 124 125 do j = 1,10 126 call swapr(A_mat(i,j),A_mat(km,j)) 127 end do 128 129 row(i) = km 130 131 c = A_mat(i,i) 132 133 if(abs(c)<1.0d-6)then 134 write(*,*) "Pivot is cloze to Zero at ",i 135 stop 136 end if 137 138end do 139 140do k = 1,10 141 142 c = A_mat(k,k) 143 A_mat(k,k) = 1.0 144 145 do j = 1,10 146 A_mat(k,j) = A_mat(k,j)/c 147 end do 148 149 do i = 1,10 150 if(i/=k)then 151 dum = A_mat(i,k) 152 A_mat(i,k) = 0.0 153 do j = 1,10 154 A_mat(i,j) = A_mat(i,j) - dum*A_mat(k,j) 155 end do 156 end if 157 end do 158 159end do 160 161do j = 10,1,-1 162 163 if(j/=row(j))then 164 do i = 1,10 165 call swapr(A_mat(i,j),A_mat(i,row(j))) 166 end do 167 end if 168 169end do 170 171end subroutine matinv 172 173 174subroutine swapr(a,b) 175 176implicit none 177real(8) :: a,b,c 178 179c = a 180a = b 181b = c 182 183end subroutine swapr

試したこと

変数の型を変えた。結果として、明らかに大きい数が結果として得られていた。

補足情報(FW/ツールのバージョンなど)

PCはWindowsで、パワーシェルを使ってFortranを利用している。

以下のような質問にはグッドを送りましょう

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

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

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

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

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

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

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

適切な質問に修正を依頼しましょう。

回答1

0

A_mat行列に関しても、一か所だけNaNとなっている個所がある。

ということであれば、どういうソフトやライブラリを使っても逆行列を求めることは不可能です。
なぜならば、NaNを含む行列の逆行列は存在しないからです。

A_matを求める過程で誤りがあるのでしょうから、それを修正しましょう。

逆行列についてはLapackの計算ルーチン: 実行列の逆行列 : (DGETRF により既に分解された行列)を使うのが簡単です。

投稿2021/11/14 04:29

ppaul

総合スコア24620

良いと思った回答にはグッドを送りましょう。
グッドが多くついた回答ほどページの上位に表示されるので、他の人が素晴らしい回答を見つけやすくなります。

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

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

このような回答には修正を依頼しましょう。

まだベストアンサーが選ばれていません

会員登録して回答してみよう

アカウントをお持ちの方は

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

ただいまの回答率
86.12%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問

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

Windows 10

Windows 10は、マイクロソフト社がリリースしたOSです。Modern UIを標準画面にした8.1から、10では再びデスクトップ主体に戻され、UIも変更されています。PCやスマホ、タブレットなど様々なデバイスに幅広く対応していることが特徴です。

FORTRAN

FORTRAN(フォートラン)は科学時術計算に向いた手続き型プログラミング言語です。 並列計算の最適化が行いやすい特性上、数値予報および気候モデルなどの大規模な計算を行う分野のスーパーコンピュータで使われています。