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
-----------------------------------------------------------------
Programming Report9
program sele
implicit none
integer,parameter::max=100000
integer::i,j,n=0,status=0
integer:: a(max),b(max),d
!Here, open two files and named each of them like 10 and 20.
open(10,file="unsort.txt");open(20,file="sort.txt")
!Do loop until the unsort.txt becomes empty.
do while (status ==0)
n = n+1
!Here, we get a(1),...a(n) from unsort.txt.
read(10,*,IOSTAT=status) a(n)
end do
n = n-1;close(10); b = a
do i = n-1,2,-1 !If N=5, i = 4.3.2.
do j = 1,i !If N=5(like above), j=1~4,1~3,1~2. This is just same as Bubble sort.
if (b(j)>b(j+1)) then !Exchange if the previous value is bigger than next one.
d = b(j+1);b(j+1) =b(j);b(j) =d !Exchange process.
end if
end do
end do
do i = 1, n
write(*,'(i3)') b(i) !Write on the screen.
write(20,*) b(i) !Write on sorted.txt here.
end do
close(20)
end program sele
programming report7 別ver
program report_7
implicit none
integer ::a,x,y,z,i,n=0
logical prime
write(*,*) 'Enter the limit of hypotenuse'
read(*,*) a
do z=1,a
do y=1,a-1
do x= 1, a-2
if (x**2+y**2==z**2) then
if (x<y .and. y<z)then
do i = 2,y
prime = .true.
if (mod(x,i)==0 .and. mod(y,i)==0) then
prime=.false.
exit
endif
end do
if (prime) then
n=n+1
write(*,*) n,')', x,y,z
endif
endif
endif
enddo
enddo
enddo
write(*,*) 'The total of primitive Pythagorean triples is', n
end program report_7
Programming Midterm
program midterm
implicit none
real:: Ax=0.0, Ay=0.0, Mxl=0.0, Mxr=0.0, Myl=0.0, Myr=0.0, SumA=0.0, SumMx=0.0, SumMy=0.0,xc=0.0,yc=0.0
integer:: i=0, N
print*, "Inout N here"
read(*,*) N
do i= 1,N-1
Ax = Ax + 1 - real(i)/real(N)
Ay = Ay + sqrt(1-(- 1+real(i)/real(N))**2)
Mxl = Mxl + 0.5*(1-(- 1+real(i)/real(N))**2)
Mxr = Mxr + 0.5*(1-real(i)/real(N))**2
Myl = Myl + (- 1+real(i)/real(N))*sqrt(1-(- 1+real(i)/real(N))**2)
Myr = Myr + real(i)/real(N)*(1-real(i)/real(N))
end do
SumA = (Ax + Ay + 1.0)/real(N)
SumMx = ( Mxl+ Mxr+ 0.5 )/real(N)
SumMy = ( Myl+ Myr)/real(N)
xc = SumMy/SumA
yc = SumMx/SumA
print *, "xc=",xc,"yc=",yc
end program midterm
Programming Example 7 (訂正版）
!This program was change on November 18.
program report7
implicit none
real:: b
integer:: a,c,i,j,k,N,t=0
logical judge
print*, "Enter the limit of hypotenuse"
read(*,*) N
print*, "The following numbers are all primitive Pythagorean triples with hypotenuse no longer than",N
do i= 2,N
c = i
do j = 1, N
a = j
b = sqrt(real(c**2-a**2))
if (b==real(int(b)) .and. b/=0 .and. real(a)<=b ) then
judge = .true.
do k =2,int(b)
if (mod(a,k)==0 .and. mod(int(b),k)==0) then
judge = .false.
end if
end do
if ( judge ) then
t = t + 1
print*, t,")",a,int(b),c
end if
end if
end do
end do
print*, "The total number of primitive Pythagorean triplese is",t
end program report7
Programming Example
program repo
implicit none
real(8)::up,xold=0.0,xnew !Up is border. Other two are terms of series.
character(16)::e !This e will be used to say conclusion correctly.
integer::i=0 !counter
print*,"Assign the numerical tolerance (e.g. 1D-6)"
read*,up !Read up here.
print*,"Input the guessed value"
read*,xnew !Read x_1 here.
write (e,'(f0.15)') up
!1d-6 can't be shown correctly if you "print up".
!So I changed it to character here.
do while (abs(xnew-xold)>=up) !This shows this loop will stop when xn converges enough.
xold = xnew !substitue xnew provided before this step for xold.
xnew = xnew - res(xnew)/dres(xnew,up) !The xnew provided before this step is changed into "new" xnew.
i = i + 1 !Add counter by 1.
print*,"Iteration =",i,"Convergent result =",xnew,"Residual =",res(xold)
end do
!This is conclusion.
print*,"The numerical solution with the precision up to",e,"is",xnew
contains
!These functions below are ones used above.
function res(x)
real(8) res, x
res = x**2.0 - 2.0**x
end function
function dres(x,up)
real(8) dres,x,up
dres = (res(x+up)-res(x))/up
end function
end program repo
Programming Example
program isbn
implicit none
logical judge
character(len=10)::a
integer::i,b,h,s=0,t,j,y,q,l
print*, 'Type in a 10-digit ISBN marked with "?" for the missing digit.'
read*, a
judge = .false. !This judge will be used to check whether including "?" or not.
b= ichar(a(10:10)) !This b will be used to check whether including "X" or not.
!This is ERROR message when user input wrong number.
if (len_trim(a)/=10 ) then
print*, "Program is terminated deu to incorrect input"
stop
end if
!This "do" checks whether input includes "?" or not.
do i =1,10
if (ichar(a(i:i))==63 ) then ! "?" is 63 in command
judge = .true.
h = i !This shows there is "?" at "h"th place.
end if
end do
if (judge) then
if ( b==120 .or. b==88 ) then !x=120 X=88
!This is when input include "?" and "X".
do i = 1,h-1
read(a(i:i),*) t
s = s + t*(11-i)
end do
do j = h+1,9
read(a(j:j),*) t
s = s + t*(11-j)
end do
s = s + 10
else
!This is when input include "?" but "X".
do i = 1,h-1
read(a(i:i),*) t
s = s + t*(11-i)
end do
do j = h+1,10
read(a(j:j),*) t
s = s + t*(11-j)
end do
end if
!This "do" finds what is "?" from 0 to 10.
do q =0,10
s = s + q*(11-h)
if ( mod(s,11)==0 ) then
l = q
exit
end if
s = s - q*(11-h)
end do
!This is final solution when input includes "?"
if (l==10) then
print*,"The missing digit is X"
else
print*,"The missing digit is",l
end if
stop
end if
!Situation above is When the setences include "?"
!Below one is When the input doesn't include "?"
if ( b==120 .or. b==88 ) then !Checking whether input includes x or not.
do i = 1,9
read(a(i:i),*) t
s = s + t*(11-i)
end do
s = s + 10
else
do i = 1,10
read(a(i:i),*) t
s = s + t*(11-i)
end do
end if
!Final answer
if ( mod(s,11)==0 ) then
print*, "This ISBN is valid"
else
print*, "This ISBN is invalid"
end if
end program isbn
