Programming Example 10
加筆:
Jacobi メソッドとGauss-Seidelの中でわざわざTempを使ってるのはA(i:j)の形で一々使うと、メモリ効率が悪い(どっかに書いてあった)ので、置き換えただけです。
I used Temp in two sessions Jacobi and Gauss-Seidel method in order to avoid the situation A(i:j) uses the high rate of RAM performances.
Someone, I don’t recall who well, said that using A(i:j) many times causes ineffective performance.
!Yuma
!This proram has many "dirty" expressions. Don't care.
program report10
implicit none
integer,parameter:: t = 10000
real(8)::A(t,t),B(t),X0(t)=0.0,X(t)=0.0,Jres,Temp
integer(8)::N,i,j,jmax,r,k,c=0.0
logical judge
print*, "Please set the number of columns and rows (n-n matrix).(N<=10000) "
read*, N
!Set format20
20 format(I2,(a),I2,3x,(a))
!Read Matrix A
do i = 1, N
do j = 1, N
write(*,20) i,"*",j, "On matrix A"
read(*,*) A(i,j)
end do
end do
!Read vector B
do i = 1, N
write(*,20) i,"*",n, "On vector B "
read*, B(i)
end do
!Set Format10 which will be used repeatedly
10 format(f10.3)
!Display the original A and B
print*, "Matrix and vector before pivoting"
do i =1, N
write(*,fmt='(a)',advance='no') "|"
do j = 1,N
write(*,10,advance='no') A(i,j)
end do
write(*,fmt='(a)',advance='no') "|"
write(*,10) B(i)
end do
print*, "-----------------------------------------------------------------"
!Judge whether diagnole has 0 or not on matrix A
judge = .false.
do i = 1, N
if (A(i,i)==0) then
judge = .true.
end if
end do
!Reorder the matrix A
if (judge) then
do r =1,N-1
jmax =A(r,r)
do j = r + 1,N
if (abs(jmax)<abs(A(j,r))) then
jmax = j
end if
end do
if (jmax /= A(r,r)) then
do k = 1,N
A(r,k)= A(r,k)+A(jmax,k);A(jmax,k)=A(r,k)-A(jmax,k);A(r,k)=A(r,k)-A(jmax,k)
end do
b(r)=b(jmax)+b(r);b(jmax)=b(r)-b(jmax);b(r)=b(r)-b(jmax)
end if
end do
end if
!Display the reordered A and B
print*, " "
print*, "Matrix and vector after pivoting"
do i =1, N
write(*,fmt='(a)',advance='no') "|"
do j = 1,N
write(*,10,advance='no') A(i,j)
end do
write(*,fmt='(a)',advance='no') "|"
write(*,10) B(i)
end do
print*, "-----------------------------------------------------------------"
!Jacobi Method
print*, " "
print*, "Jacobi method with the numerical tolerance 0.100E-02"
print*, "-----------------------------------------------------------------"
write(*,fmt='(a)',advance='no') "Iteration"
do i = 1, N
write(*,'(a)',advance='no') " X"
write(*,'(I1)',advance='no') i
end do
write(*,'(a)') " Residual"
print*, "-----------------------------------------------------------------"
Jres = 1.0
do while (sqrt(Jres) >1.0E-3)
write(*,'(I3)',advance='no') c
do k=1,N
write(*,10,advance='no') x(k)
end do
if (c==0) then
write(*,*) " -"
else
write(*,'(E10.3)') Jres
end if
Jres = 0
do i= 1, N
X0(i)= X(i)
end do
do i= 1, N
Temp = B(i)
do j= 1, N
if (j /= i) then
Temp = Temp - A(i,j)*X0(j)
end if
end do
X(i)= Temp/A(i,i)
end do
do i = 1,N
Jres = Jres + (x0(i)-x(i))**2
end do
c= c + 1
end do
print*, "-----------------------------------------------------------------"
!The end of Jacobi
print*, " "
c = 0
!Gauss-Seidel
print*, "Gauss-Seidel method with the numerical tolerance 0.100E-02"
print*, "-----------------------------------------------------------------"
write(*,fmt='(a)',advance='no') "Iteration"
do i = 1, N
write(*,'(a)',advance='no') " X"
write(*,'(I1)',advance='no') i
end do
write(*,'(a)') " Residual"
print*, "-----------------------------------------------------------------"
do i =1, N
x(i)=0
end do
Jres = 1.0
do while (sqrt(Jres)>1.0E-3)
write(*,'(I3)',advance='no') c
do k=1,N
write(*,10,advance='no') x(k)
end do
if (c==0) then
write(*,*) " -"
else
write(*,'(E10.3)') Jres
end if
Jres = 0
do i= 1, N
X0(i)= X(i)
end do
do i= 1, N
Temp= B(i)
do j= 1, N
if (j /= i) then
Temp = Temp - A(i,j)*X(j)
end if
end do
X(i)= Temp/A(i,i)
enddo
do i = 1,N
Jres = Jres + (x0(i)-x(i))**2
end do
c = c + 1
end do
print*, "-----------------------------------------------------------------"
end program report10
Result
Please set the number of columns and rows (n-n matrix).(N<=10000)
4
1* 1 On matrix A
2
1* 2 On matrix A
0
1* 3 On matrix A
2
1* 4 On matrix A
-3
2* 1 On matrix A
0.25
2* 2 On matrix A
0
2* 3 On matrix A
0.25
2* 4 On matrix A
6
3* 1 On matrix A
-5
3* 2 On matrix A
1
3* 3 On matrix A
-1.5
3* 4 On matrix A
2
4* 1 On matrix A
1
4* 2 On matrix A
1
4* 3 On matrix A
1
4* 4 On matrix A
1
1* 4 On vector B
9
2* 4 On vector B
7.5
3* 4 On vector B
-18
4* 4 On vector B
10
Matrix and vector before pivoting
| 2.000 0.000 2.000 -3.000| 9.000
| 0.250 0.000 0.250 6.000| 7.500
| -5.000 1.000 -1.500 2.000| -18.000
| 1.000 1.000 1.000 1.000| 10.000
-----------------------------------------------------------------
Matrix and vector after pivoting
| -5.000 1.000 -1.500 2.000| -18.000
| 1.000 1.000 1.000 1.000| 10.000
| 2.000 0.000 2.000 -3.000| 9.000
| 0.250 0.000 0.250 6.000| 7.500
-----------------------------------------------------------------
Jacobi method with the numerical tolerance 0.100E-02
-----------------------------------------------------------------
Iteration X1 X2 X3 X4 Residual
-----------------------------------------------------------------
0 0.000 0.000 0.000 0.000 -
1 3.600 10.000 4.500 1.250 0.135E+03
2 4.750 0.650 2.775 0.912 0.918E+02
3 3.263 1.562 1.119 0.936 0.579E+01
4 3.951 4.682 2.642 1.067 0.125E+02
5 4.171 2.339 2.150 0.975 0.579E+01
6 3.813 2.704 1.792 0.987 0.390E+00
7 3.998 3.408 2.167 1.016 0.671E+00
8 4.038 2.819 2.027 0.993 0.369E+00
9 3.953 2.942 1.952 0.997 0.281E-01
10 4.002 3.098 2.043 1.004 0.352E-01
11 4.008 2.951 2.004 0.998 0.232E-01
12 3.988 2.989 1.989 0.999 0.210E-02
13 4.001 3.023 2.011 1.001 0.181E-02
14 4.002 2.987 2.000 1.000 0.143E-02
15 3.997 2.998 1.997 1.000 0.157E-03
16 4.000 3.006 2.003 1.000 0.910E-04
17 4.000 2.997 2.000 1.000 0.874E-04
18 3.999 3.000 1.999 1.000 0.117E-04
19 4.000 3.001 2.001 1.000 0.448E-05
20 4.000 2.999 2.000 1.000 0.525E-05
-----------------------------------------------------------------
Gauss-Seidel method with the numerical tolerance 0.100E-02
-----------------------------------------------------------------
Iteration X1 X2 X3 X4 Residual
-----------------------------------------------------------------
0 0.000 0.000 0.000 0.000 -
1 3.600 10.000 4.500 1.250 0.135E+03
2 4.750 0.650 2.775 0.912 0.918E+02
3 3.263 1.562 1.119 0.936 0.579E+01
4 3.951 4.682 2.642 1.067 0.125E+02
5 4.171 2.339 2.150 0.975 0.579E+01
6 3.813 2.704 1.792 0.987 0.390E+00
7 3.998 3.408 2.167 1.016 0.671E+00
8 4.038 2.819 2.027 0.993 0.369E+00
9 3.953 2.942 1.952 0.997 0.281E-01
10 4.002 3.098 2.043 1.004 0.352E-01
11 4.008 2.951 2.004 0.998 0.232E-01
12 3.988 2.989 1.989 0.999 0.210E-02
13 4.001 3.023 2.011 1.001 0.181E-02
14 4.002 2.987 2.000 1.000 0.143E-02
15 3.997 2.998 1.997 1.000 0.157E-03
16 4.000 3.006 2.003 1.000 0.910E-04
17 4.000 2.997 2.000 1.000 0.874E-04
18 3.999 3.000 1.999 1.000 0.117E-04
19 4.000 3.001 2.001 1.000 0.448E-05
20 4.000 2.999 2.000 1.000 0.525E-05
-----------------------------------------------------------------
