FORTRANによる逆行列の計算

前提・実現したいこと

多項式曲線フィッティングを行うために、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