Programming Example 11
Program report11
implicit none
character(len=50):: matrixfile
integer:: status=0, n
real(8):: r
real(8), dimension(:,:), allocatable:: Matrix, Inv_Matrix, Identity
integer:: i, j, row
write(*,*) "Enter the file name of matrix"
read(*,*) matrixfile
!===== Read file =====
open(1, file=matrixfile)
n=0
do
read(1, *, iostat=status) r
if (status/=0) exit
n=n+1
end do
rewind(1)
allocate(Matrix(n,n))
do i=1, n
read(1, *) Matrix(i, :)
end do
close(1)
!===== Process data =====
allocate(Inv_Matrix(n,n), Identity(n,n))
Inv_Matrix = Minverse(Matrix)
Identity = matmul(Matrix, Inv_Matrix)
!Check
print*, "The matrix before inversion."
do i=1, n
write(*,'(a)',advance="no") "|"
do j=1, n
write(*,'(E12.3)',advance="no") Matrix(i,j)
end do
write(*,*) "|"
end do
!===== Write file =====
open(2, file="invmat.dat")
print*,"The matrix after inversion"
do i=1, n
write(*,'(a)',advance="no") "|"
do j=1, n
write(*,'(E12.3)',advance="no") Inv_Matrix(i,j)
end do
write(*,*) "|"
end do
print*, "Check the identity matrix"
do i=1, n
write(*,'(a)',advance="no") "|"
do j=1, n
write(2,'(E12.3)',advance="no") Inv_Matrix(i,j)
write(*,'(E12.3)',advance="no") Identity(i,j)
end do
write(2,*) " "
write(*,*) "|"
end do
close(2)
Contains
Function Minverse(A)
real(8), dimension(:,:):: A
real(8), dimension(:,:), allocatable:: Minverse, B, V
real(8), dimension(:), allocatable:: temp
real(8), parameter:: TOL=1D-6
real(8):: pivot
integer:: i, j, jmax, k, n,l
n=size(A)**0.5
allocate(Minverse(n,n),B(n,n),V(n,n),temp(n))
B=A
V(:,:)=0.0
do i=1, n
V(i,i)=1.0
end do
do i=1, n
jmax=i
do j=i+1, n
if (abs(B(j,i))>abs(B(jmax,i))) then
jmax=j
temp=B(i,:);B(i,:)=B(jmax,:);B(jmax,:)=temp
temp=V(i,:);V(i,:)=V(jmax,:);V(jmax,:)=temp
end if
end do
pivot=B(i,i)
if (abs(pivot)<TOL) then
print *, "The matrix is singular!"; stop
end if
B(i,:)=B(i,:)/pivot !All pivots on B(i,:) are divide by diagnole in order to set diagnole 1.
V(i,:)=V(i,:)/pivot !Apply the above process for V which leads to Inverse of the matrix.
!This is the wrong code I wrote.
! do j=1,n
! B(i,j)=B(i,j)/B(i,i)
! V(i,j)=V(i,j)/B(i,i)
! end do
!If we failed to set pivot=B(i,i), the value is changeable so we can't get correct value.
do k =1,n
if (i .ne. k) then
pivot=B(k,i) !Set pivot B(k,i) which means the number located on the column the above process works.
do l=1,n
B(k,l)=B(k,l)-B(i,l)*pivot !B(k,i) is changed to new row in order to change one column into 0 .
V(k,l)=V(k,l)-V(i,l)*pivot !Apply above process for V.
end do
end if
enddo
end do
!This is the wrong code I wrote.
! do l=1,n
! B(k,l)=B(k,l)-B(i,l)*B(k,i)
! V(k,l)=V(k,l)-V(i,l)*B(k,i)
! end do
!If we failed to set pivot=B(k,i),the same thing as above will happen.
Minverse=V
end Function Minverse
End Program report11
京大地球工学科国際コースについて【1】
国際コースに関する記事は移行しました(2020/02/13)
こちらへ
http://yuma2012.hateblo.jp/entry/2020/02/13/042247
眠い、前文省略します。
もともと文がうまくないんだから、ここで面白さを求めても駄目だわ。
今はインドネシア人の友人と原発のプレゼンを作ってます。あまりに僕のやる気がなく、挙句の果てにブログを書き始めたのでとうとうキレてます。
ここでは僕の所属する京都大学工学部地球工学科国際コースとは何たるかを書きます。なんでこんなの書くかといったら簡単。来年人が全然入らないと寂しいからです。では早速。
１．地球工学科とは？
まずは、大前提から。京大の工学部には６つの科があります。その一つが地球工学科で、土木・資源・環境の３分野で成り立っています。土木というとよくわからないかもしれませんが、端的に言ってインフラです。僕の勝手なイメージですが、土木＝インフラ整備（ハードの整備）＋計画系（ソフトの整備）な感じ。インフラはご想像の通り、ダムとか道路とか沿岸整備とか。計画は都市計画とか交通マネジメントとかですかね。ちなみに僕は交通か防災に行きたいなぁ。そういえば、内閣官房付になった藤井教授は交通系だったはず。
まあそんな感じ。ちなみに僕は合格してから内容を知った。
2.国際コースとは？
先ほども言った通り、地球工には３つ分野がありまして、その配属が3年に上がるときにあります。国際コースは土木科の中にあるので、土木しかできませーーん。しかし、国際コースは一回生の入学前に分属されてしまいます。そのあとは出ることしかできません！逆に言えば、研究室配属まで何のコース分けもないです。
ではまず基本的な情報をまとめていきます。
・すべての授業（一般教養などすべて。第二外国語を除く）が英語。
・留学生がクラスの半分～いる。国籍はアジアがほとんど。たまにエジプトとかモンゴルとか。
・課題が過大
・発表系の授業が多め。プレゼンたくさん。熱力学でプレゼン。専門でプレゼン。実験でプレゼン。一般教養でプレゼン。少人数セミナーでプレゼン。プレレレレレレ。
・出席点の比率高め。
・ほぼすべての授業が20人前後の少人数で行われる。
・3年次に海外インターン
・卒業論文が英語
国際コースのここが素晴らしい！！！
・英語は当然スーパー上達する。わしゃもうネイティブレベル。
・面白い人が多い。
・クラス仲が良い。
・いろいろな文化の違いを感じれる
・就職余裕
・世界に羽ばたく「グローーーバル人材www」になれる。
・他の学部学科の人からかなり羨ましがられる
・ドヤれる
・イキれる
・学畜になれる
・いつでも辞めて日本語コースに移動できる。
・GORILAとかいう呪いみたいな制度から解放された自由地区
（上記は75％の事実と25％の偏見を基に作成）
国際コースのここがクソ
・課題が大変。僕みたいに要領悪いとサークル行けへん
・したがって当然部活やるとくそ大変。（ただし、ギリギリ単位だけ取るならまあいける？）
・常に立ちはだかる言語の壁
・就職そんな楽じゃない
・もうなんか色々とFUCK
・その他悪口を順次受け付けています。国際コースの人は@yuma2011nまで
はあもうなんかそんな感じ。詳しい説明は以下の通り。
☆授業～～～
ほぼすべての授業が基本的に全部英語で行われます。
土木の専門授業はもちろん地球工学の主に外国人教授が授業をされます。が、日本人の教授も英語で授業をされます。外国人教授の国籍はあまりにバラバラなので、特に傾向を感じません。日本人の教授はかなり英語で苦労されている方も多いらしく、先輩曰く教授の英語力と業績に特に相関はないらしいです。基本的に常に語学で苦労しているのは教える方も教わる方も日本人です。専門科目はほとんどペーパーのテストで評価されます。
さらに、自然科目（例えば力学や熱力学）、英語科目、情報科目なども地球工の外国人教授が担当しています。唯一、一年次の数学は数学科の先生が担当します。二年次からの数学はこれまた地球工の先生が担当します。これは驚いた。何の科目とっても同じ人たち出てくんの。草でしょ。
一般教養の人文社会系や統合・健康群などは他の学科の学生と同じように受けますが、それらもすべてE２という英語の授業の中で選んで取ります。だいたいの授業がレポートかディスカッション、授業態度で評価されます。
ここで問題があって、人社の授業は前後期合わせて700ありますが、我々が取れる英語の授業(E2科目)は90個あまりしかありません。そのうち他の科目と被らないように取ると、選択肢はかなり狭まります。しかも、来年から京都大学はE2科目の人数制限を厳しくするそうです。ますます、人社E2事情は厳しくなるのです。クソ京都大学マジで死
ええと...気を取り直して。
唯一、第二外国語だけは他の学生と完全に同じです。ちなみに僕は韓国語。二コマのうち一つは日本人、もう一つは韓国人の教授なので、今季の僕が取っている授業で日本人の教授は一人だけです。
面倒になったのでこれから先のことは次回に回します。文体が荒れてすみません。なにぶん暇つぶしに書いたので。反省しません。
えええと国際コースは今年でできてから7年目ですので、まだまだ僕たちはマウスですね。これからどうなるかわかりませんが、超面白いので（笑）ぜひ入ってください。
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
