Nagoya Institute of Technology

Nagoya Institute of Technology
連続体シミュレーション工学
Fortran 入門
名古屋工業大学 機械工学科
計測プログラム
後藤俊幸
Nagoya Institute of Technology
はじめに
Fortran （FORmula TRANslator）
IBM が数値計算を目的に開発した世界初の高級プログラミング言語 (1950）
・ 数学演算に強い (行列演算、複素数、各種組み込み関数、etc）
・ C、C++などよりも数値計算ではパフォーマンスに優れている
・ 多くの数値計算ライブラリの資産
・ 並列化も容易
・ スパコンでのシミュレーションでも高パフォーマンスを発揮
Nagoya Institute of Technology
Fortran の歴史
Fortran66
世界初のプログラミング言語の標準
Fortran77
固定形式 広く普及 多くの数値計算ライブラリが充実
(ベクトル型スパコン）
Fortran90
抜本的変更 自由形式 構造化 （Fortran77文法を包含）
配列演算の強化
(ベクトル並列型スパコン）
Fortran95
マイナーチェンジ
Fortran2003 大きく変更
(Fortran90/95
を包含）
Fortran2008 並列処理の標準化
（大規模スカラー並列型スパコン）
Nagoya Institute of Technology
簡単な例
（自由形式）
! Computation of three kinds of averages 注釈行
Program Example_1_1
プログラム文 （なくてもよい）
Real:: a,b,av1,av2,av3
型宣言文 すべての変数には型がある
print *, ‘enter two numbers’
画面（＊：標準出力）に表示
Read *, a,b
データを読み込む（＊：キーボード）
av1=(a+b)/2.d0
計算
av2=(a*b)**0.5
av3=2.0/(1.0/a+1.0/b)
print *, av1,av2,av3
データを画面に書き出す
End Program Example_1_1
End 文
Nagoya Institute of Technology
Fortran まずは使ってみよう
手続き
１．エディターでプログラムを書く
各自が好きなエディターで先ほどの例題（教科書 ｐ１）を入力
フォルダーに格納 p11.f90 （例えば「Fortran」という名称のフォルダー）
２．プログラムをFortranでコンパイル
コマンドプロンプトを立ち上げる （講義中はこのウィンドウ内でしか使えません）
スタート ＞ すべてのプログラム ＞ アクセサリ ＞ コマンドプロンプト
コマンドプロンプト： 右クリック ＞ 送る ＞ デスクトップ（ショートカット作成）
ここで、cd コマンドで Fortran のあるディレクトリに移動
g95 p11.f90 （PCでは g95 という名称のコマンドになっています）
エラーがなければ a.exe という実行ファイルができる
３．実行
a.exe と入力
Nagoya Institute of Technology
ディレクトリ変更
コンパイル
実行
Nagoya Institute of Technology
Fortran 文法
第1回
Nagoya Institute of Technology
使える文字
英字：A～Z（大文字・小文字どちらも使用可能。大文字と小文字は区別されない）
見やすいことが大事
数字：0～9
下線：_
特殊文字：（空白） ＝ ＋ － ＊ / ( ) ， ． ’ ： ！ ” ％ ＆ ； ＜ ＞ ？$
書き方
固定形式： 先頭から7コラム目から書き始め72コラムまで
ファイル名は sample.for (Windows) sample.f （Unix)
自由形式： 1行をどこから始めてもよい
ファイル名は sample.f90
この講義では自由形式で解説
Nagoya Institute of Technology
自由形式
行： 1行 132文字以内 ； で区切ると1行に複数の文が書ける
x=1.d0 ; y=2.d0 ; z=3.d0
注釈行： 先頭が ！ で始まる行、または空白行は注釈
文の途中の ！ があると、それより右が注釈となる
! This is an example
y=sin(x)
! Sine function
継続行： 長い行の場合には 行の最後に & をつけると次の行に継続される
y= 1.0+ 2.0*x*x+ cos(x) &
+ 0.5*sin (2.0*x)
空白：
実行文や命令語 以外では空白を入れてもよい
○ endif
○ end if
○ forall
× for all
Nagoya Institute of Technology
プログラム構成
main program
[Program プログラム名]
…….
…….
end [Programプログラム名]
sub program
プログラム単位
省略可
Example
Program main
…..
call time_advance (u1,u2,u3)
…..
End program main
プログラム単位
副プログラム種別名 副プログラム名 省略不可
…….
…….
end [副プログラム種別名 [副プログラム名] ]
Subroutine time_advance (u1,u2,u3)
……
return
End
Nagoya Institute of Technology
文の種類
宣言文： 実行に必要なメモリ領域の大きさなどを設定
real*8 :: a,b,c ; integer :: i,j,k
dimension (10,10):: array
非実行文
実行文： 実行時に計算や処理の流れを制御
do i= 1,10
a(i)=2*I
if (I > 3) exit
end do
文の順序
プログラム単位開始文 （program 文、subroutine 文, function文）
宣言文 (型宣言文、Common 文）
Data 文
実行文
Contains 文
End (1つのプログラム単位内に必ず１つ必要）
Nagoya Institute of Technology
プログラム単位名、変数名の名前の付け方
英字で始まる31文字以下の 英字、下線、数字で構成
vel1, vel2, stream_function, Potential_1
random_number_generator
変数と定数
変数： プログラムの実行によってデータの値が変わるもの
スカラー変数、 ベクトル変数、 配列変数
どの変数がどのような型 (整数、単精度実数、倍精度実数、…)かが大事
定数： プログラムの実行中でもデータの値が変わらないもの
整定数
通常 4バイト (8×4＝３２bit）
-231=-2147483648 ～ 231-1= 2147483647
実変数
単精度 ４バイト、、倍精度 ８バイト
小数形式 3.141592
指数形式 7.199847e-12
複素定数 ２つの４バイト実数の組
論理定数 ４バイト
.true.
2+3i=(2.e0, 3.e0)
.false.
文字定数 １文字が１バイト、 日本語は１文字２バイト
‘Sample’ “例題”
の範囲
Nagoya Institute of Technology
構文名
loop1 :
loop2 :
do j=1,n2
do i=1,n1
…………..
この部分がとても長いときや内部に
多くの do loop が混在するときに有効
…………..
end do: loop2
end do: loop1
Nagoya Institute of Technology
数値計算
代入文
変数＝変数、定数、数式
= は右辺の計算値を左辺の変数に代入するという意味
Program summation
integer :: sum
sum=0
do i=1,100
sum = sum+1
end do
print *, sum
end
和を計算する典型的な計算スキーム
Nagoya Institute of Technology
算術演算
足し算
引き算
掛算
割り算
べき乗
+
–
*
/
**
a+b
a-b
a*b
a/b
x**3 = x*x*x
演算の優先順位
足し算、引き算 < 掛算、割り算 < べき乗
（ ）がなければ左から右に向かって演算、 怪しいときには（ ）を使う
x+y*z= x+(y*z)
x/y/z= (x/y)/z
x*y**2= x*(y**2)
x**y**2=x**(y**2)
Nagoya Institute of Technology
実行の順序
文字列の出力
上から順位実行される
print *, ` 出力したい文字列‘
integer :: s
s=10
s=s+1; print *, s
s=s+1; print *, s
s=s+1; print *, s
end
real, parameter :: pi=3.14159265
real:: r,s
read *,r
s=pi*r*r
print *, `radius, area`, r, s
end
Nagoya Institute of Technology
練習問題 １
正の実数 x と正の整数 n を読みこんで、 x**n と x**(1.e0/n)
を求める
練習問題 2
摂氏の温度 C から華氏温度 F を求めること
F=1.8C +32.0
練習問題 ３
伸長 L (m) と体重 w (kg) から BMI （肥満度指数） を求めること
BMI = w / L2
Nagoya Institute of Technology
練習問題 プログラム例
! ex1
real:: x,xn,xn1
integer:: n
read *,x,n
x1=x**n
xn1=x**(1.e0/n)
print *,x,n,x1,xn1
end
! ex2
real::c,f
read *,c
f=1.8*c +32.0
print *,c,f
end
! ex3
real::L,w,bmi
print *,'Enter W and L'
read *,w,L
bmi=w/(L**2)
print *,'weight[kg]=',w,', Hight[m]=',L, ', BMI=',bmi
end
Nagoya Institute of Technology
第２回
Nagoya Institute of Technology
流れを変える
if, stop, goto, case
1 if 文
if
（ 論理式、関係式、論理演算式 ） 実行文
（
） の中が真の時だけ実行文を実行する
Integer n, s
n=32
s=1
if ( n>=10 ) s=2*s +1
2 stop 文
stop [ ‘メッセージ文’]
end 文にまで行かなくてもその場で停止
Integer n
if ( n<0 ) stop ‘n は正の整数でないといけません’
Nagoya Institute of Technology
3 goto 文
goto 文番号 （ go to でもよい）
指定された文番号の実行文までジャンプする
10
Integer:: n, s, k, sum
k=0; s=0; sum=0
print *,’enter students mark (quit, if negative integer)‘
read *, n
if ( n<0 ) go to 999
s=s+1
sum=sum+n
if ( n >=60)
k=k+1
go to 10
999 print *,‘Total and No . of students passed this exam. =‘,s,k
print *,’Average of score =‘, sum/s
end
Nagoya Institute of Technology
関係演算子
>
<
<=
==
/=
または
.gt.
.lt.
.le.
.eq.
.ne.
大きい
小さい
小さいか等しい
等しい
等しくない
if( n+m>0 ),
if( a<4.5 )
関係演算子の両辺の変数の型は一致していること
論理演算子
.not.
.and.
.or.
.eqv.
.neqv.
論理否定
論理積
「かつ」
論理和
「または」
論理等価
論理非等価
n=0 または 1
0 < k < 8
x=y=z
(n==0).or.(n==1)
(0<k).and.(k<8)
(x==y).and.(y==z)
Nagoya Institute of Technology
演算子の優先順位
.neqv. と .eqv. < .or. < .and. < .not. < 関係演算子 < 算術演算子
できるだけ （ ）を用いる
優先順位が同等なら 左から優先的に判断される
x+y> 1.0e0  (x+y) > 1.e0
(x<3.14) .or. (z<sqrt(2.e0))
Nagoya Institute of Technology
If 構文
もっと機能的な if 文
［構文名：］ if ( 論理式 ) then
………
実行文
………
endif ［構文名］
構文名を付けたときにはいつもペアで用いる
if ( 論理式 ) then
………
else
………
endif
飛出し可
if ( 論理式 ) then
………
else if ( 論理式 ) then
………
else
飛び込み禁止！
………
endif
• 分岐はいくらあってもよい
• if 構文の中にいくつもの if 構文があってもよい （多重if 構文）
• if 構文の中に外から飛び込みがあってはいけない
Nagoya Institute of Technology
入れ子構造
block1: if ( 論理式 ) then
………
block2:
if ( 論理式 ) then
………
end if blcok2
………
end if block1
block1: if ( 論理式 ) then
………
block2:
if ( 論理式 ) then
………
end if blcok1
………
end if block2
Nagoya Institute of Technology
例１
! Calculating various functions
real :: x
print *,'input a real number:'
read *,x
if(x>0.0) then
print *, sqrt(x), log(x), log10(x)
else
print *,' Not defined for x<0‘
end if
end
組み込み関数
Nagoya Institute of Technology
例２
1
0
! Triangle pulse
real:: x,f
print *,'input x:'
read *,x
if(x<0.0) then
f=0.0
else if (x<1.0) then
f=x
else if (x<2.0) then
f=-x+2.0
else
f=0.0
end if
print *, x,f
end
1
2
x<0 はすでに除かれている
x<1 まではすでに除かれている
x<2 まではすでに除かれている
Nagoya Institute of Technology
case 構文
多数の場合わけに便利
［構文名：］ select case ( 式 )
case ( 式の値、または式の値の範囲 ) ！A
実行文
case ( 式の値、または式の値の範囲 ) ！B
実行文
case default
！ （上のA,B以外すべて）
実行文
end select ［構文名］
場合分けに重なりが
ないこと
Nagoya Institute of Technology
例
integer :: i
print *, ‘ Enter mark ‘
read *, I
if ( i<0 .or. i>100 ) stop
select case (i)
case (90: )
print *, ‘秀’
case (80:89)
print *, ‘優’
case (70:79)
print *, ‘良’
case (60:69)
print *, ‘可’
case default
print *, ‘不可’
end select
end
Nagoya Institute of Technology
組み込み関数
数学演算でよく使われる関数
ｘ, y：実数、
ｚ：複素数、
exp (x)
総称名
sqrt(x)
abs(x)
sin(x)
exp(x)
log(x)
log10(x)
引数
i, j, k ： 整数
総称名：引数 xの型に合わせて値を返す
ｘ：単精度実数
exp（ｘ）：単精度実数
ｘ：倍精度実数
exp（ｘ）：倍精度実数
ｚ：単精度複素数
exp (z) ：単精度複素数
ｚ：倍精度複素数
exp (z) ：倍精度複素数
z=cmplx(x,y)
real(z)
aimag(z)
conjg(z)
aint(x)
int(x)
anint(x)
Nint(x)
xの整数部 を実数で返す
xの整数部を整数で返す
xを四捨五入して実数で返す
xを四捨五入して整数で返す
Nagoya Institute of Technology
例
! Calculating various functions
real :: x,y
complex::z
print *,'input a real number:'
read *,x
if(x>0.0) then
print *, sqrt(x), log(x), log10(x)
else
print *,' Not defined for x<0‘
end if
print *,'input a complex number x,y (=x+iy):‘
read *, x,y
z=cmplx(x,y)
print *, sqrt(z), exp(z), log(z)
end
ノート
複素関数の定義
例えば z= i にとると
Nagoya Institute of Technology
X＝2
と (x,y)=(0,1）を代入した例
Nagoya Institute of Technology
練習問題 １
正の実数 x を読みこんで、 整数部 と小数部 をもとめること
練習問題 2
2次方程式 ax2+bx+c=0 の a,b,c を入力し、判別式の値をもとに
解を実数の範囲で求めよ。
Nagoya Institute of Technology
練習問題 プログラム例
! ex4
real:: x,y
integer:: n
read *,x
n=int(x)
y=x-aint(x)
print *,x,n,y
end
! ex5
real::a,b,c,d,x1,x2
print *,'Enter a,b,c'
read *, a,b,c
d=b*b-4.0*a*c
if(d>=0) then
x1=0.5*(-b+sqrt(d))/a
x2=0.5*(-b-sqrt(d))/a
print *, x1,x2
else
print *, 'No real solution'
endif
end
Nagoya Institute of Technology
第３回
Nagoya Institute of Technology
Data の型
すべての変数には型がある
integer :: 変数名のリスト
（単精度整数 4バイト）
real :: 変数名のリスト
（単精度実数 4バイト）
real*8 :: 変数名のリスト
（倍精度実数 8バイト）
complex :: 変数名のリスト
（単精度複素数 8バイト）
complex*16 :: 変数名のリスト
（倍精度複素数 16バイト）
real :: x,y
real*8 ::d1,d2
complex::z
integer :: n,m,l
Nagoya Institute of Technology
初期値設定
型宣言文の中で変数の初期値を与える
（プログラム単位が呼ばれたときに1回だけ行う）
Real :: pi=3.14159265,s,r
Integer :: n=30
Print * ‘Enter r’
Read *, r
S=pi*r*r
n=n+1
Print *, n,r,s
End
Real :: pi=3.14159265,s,r
Integer :: n
Print * ‘Enter r’
Read *, r
S=pi*r*r
n=30
n=n+1
Print *, n,r,s
End
左は n に1回だけ３０という値が設定される
右は n に３０という値が代入される
Nagoya Institute of Technology
暗黙の型宣言
i,j,k,l,m,n
で始まる変数 整数型
それ以外
Kmin, Kmax, jnum, num
型の変換
組み込み関数を用いる
x=real(k)
n=integer(s)
Z=cmplx(a,b)
実数型
Nagoya Institute of Technology
!----Solar year
----real ::year,day,jikan,hun,byo
integer :: d,h,m
year=365.2422
day=aint(year)
d=int(day)
jikan=24.0*(year-day)
h=int(jikan)
hun=60.0*(jikan-aint(jikan));
m=int(hun)
byo=60.0*(hun-aint(hun))
print*,'1yearis',d,'days',h,'hours',m,'minutes', &
byo,'seconds‘
end
Nagoya Institute of Technology
暗黙の型宣言の変更
implicit 型 （先頭文字、またはその範囲），・・・・
implicit real*8 (a-h, o-z)
implicit complex (z)
先頭文字がa からh, o-zで始まる変数は
すべて倍精度実数型
先頭文字がzで始まる変数はすべて複素数型
暗黙の型宣言の無効化
implicit none
すべての変数の型宣言をしないとコンパイルエラーになる
プログラムミスが減るので、積極的に利用するとよい
Nagoya Institute of Technology
文字型と論理型
文字型宣言
文字列データ宣言には型と長さ（文字数）必要
character ([len=]長さ) :: 変数名リスト
character :: 変数名１*長さ，変数名２*長さ，・・・・
character (len=8) :: c=‘doraemon’
character (2) :: fname, file, a
論理型宣言
logical :: 変数名リスト
logical :: flag, index, range
flag=.true.
index=.false.
range =(8.0<x).and.(x<10.0)
cが文字変数 doraemonがその内容
Nagoya Institute of Technology
例題
Program Logical_example
Implicit none
logical :: hantei
real:: a=0.e0, b=1.e0, x
read *, x
hantei=(a<x).and.(x<=b)
print *, hantei
End
判定の値は真(ture) なら T
偽(false) なら Fを出力
Nagoya Institute of Technology
名前付き定数
parameter 宣言
プログラム中で用いる定数に名前を付ける
型，parameter :: 定数の名前＝定数の値または式
integer, parameter :: mask=2147483467
real , parameter :: pi=3.14159265, pi2=2.e0*pi
character(1), parameter :: sp=‘ ‘, ap=“’”
character(len=*), parameter :: f=‘This is Fortran95’
* にすると、実際に設定された長さになるので便利
Nagoya Institute of Technology
複素数型, 高精度型
complex :: 変数名リスト
complex*16 :: 変数名リスト
real*8 :: 変数名リスト
real([kind=]8) :: 変数名リスト
real(8) :: deg,rad,s,c
real(8),parameter :: pi=3.141592653589793
print*,'Input angle variable in degree:‘
read*, deg
rad=pi*deg/180.0
s=sin(rad)
c=cos(rad)
print*,'SIN=',s,', COS=',c
print*,'COS^2+SIN^2=',s**2+c**2
end
Nagoya Institute of Technology
!----Complex number
----real :: x
complex :: z
integer :: n
complex, parameter :: i=(0.0,1.0)
!
print*,'Input a real number x and an integer n:'
read *,x,n; z=exp(i*x)
print*,z**n
print*,cos(n*x),sin(n*x)
end
Nagoya Institute of Technology
異なる型の混合演算と代入
強い方の型に変換して型をそろえたのち、演算、代入を行う
整数型 < 実数型 < 複素数型,
単精度 < 倍精度 < 4倍精度
整数型 / 整数型 は整数型
・データが短い変数は下位に０が代入される
・実数を複素数にするときには虚数部に０が入る
7/3
 2
1/3
 0
2**(-2)
 1/(2**2) 0
N-m*(n/m) 正の整数nをmで割ったときの余り
X**n 
x*x*....*x
X**y 
exp(y*log(x))
xをn回かける
Nagoya Institute of Technology
異なる型の代入
右辺を計算した後左辺の型に変換して代入
real :: a,b
real*8 :: x
a=0.0
b=2.14
x= a+b
a+b は4バイト実数同士の足し算,
xには下位に0を付けて8バイトに変換して格納される
print *, a,b,x
end
Real ::x
X=1.0e0/3
 1.e0/3.e0  xに代入
X=(3.14e0, 2.56e0)
X=5
 xに3.14e0がはいる
 xに5.e0 が入る
Nagoya Institute of Technology
Real::x
Real(8) ::y
y=1.0e0/3.e0  1.e0/3.e0 単精度計算してさらに下位に0をつけてyに代入
y=2.0d0/3.e0  2.0d0/3.0d0 倍精度計算をしてyに代入
X=(3.14e0, 2.56e0)  xに3.14e0が入る
X=5
 xに5.e0 が入る
integer ::i
i=-3.623  -3が入る
i=(3.14e0, 2.56e0)  iに3が入る
i=5.0e0/3
 iに1が入る
complex ::z
z=1.0e0  z=(1.0e0, 0.0e0)
z=5.e0
 z=(5.0e0, 0.0e0)
Nagoya Institute of Technology
関数と引数の型
関数が総称名の場合
y=f(x) : 引数xに応じた型が f に当てはまる。その後 yの型に合わせて代入
x: 倍精度
f: 倍精度
x: 倍精度複素数
f: 倍精度複素数
Nagoya Institute of Technology
レポート課題１
テキスト p.42, 問題 4-2、4-3，4-4、4-5
自分でプログラムを作成し、これをレポートに書きだす。
値を入力し、得られた計算結果をレポートに書き出す。
表紙と考察を付けること
締切
： 次回の講義時 10:30 まで
提出先 ： ２号館Ｂ棟４階４２４Ｂ室前の
レポート回収箱
レポート課題番号
学籍番号
氏名
提出日
Nagoya Institute of Technology
第4回
Nagoya Institute of Technology
繰り返し
膨大な計算を高速に行う 計算機の最も強力な能力
単純 do ループ
［構文名：］ do do変数=開始値, 限界値［, 増分値］
・・・・
end do ［構文名］
省略すると １
負の値も可能
ループを抜け出た時点では n=101
開始値が限界値を超えて
いるときは実行しない
Nagoya Institute of Technology
多重do ループ と
入れ子構造
Nagoya Institute of Technology
ループ から抜けでる
End do の直後の実行文に移る
End do に移る
Nagoya Institute of Technology
integer :: m,n,mn
var1:do m =2,9
var2: do n =1,9
print *,m,'*',n,'?'
read *,mn
if(mn==m*n)then
print *,‘Correct'
else
print *,‘wrong'
end if
end do var2
end do var1
end
Nagoya Institute of Technology
出力の書式指定
膨大な計算結果を見やすく配置する
Print 書式識別子, 出力データリスト
Write(*,書式識別子) 出力データリスト
書式識別子
* ： 標準の書式
‘(書式項目のリスト)’
（入力：キーボード、 出力：ディスプレイ）
(書式項目のリスト) は文字定数になっている
空白
文番号
書式項目のリスト
を含めて
文字
文はプログラム単位の
どこにあってもよい
Nagoya Institute of Technology
書式
編集記述子
Iw
Fw.d
w, d, n, r を整数として
整数データ 文字幅w 右詰め
I5 : 56789
実数データ 文字幅w 右詰め 少数点以下d桁 w≥ d+2
d
F5.2 ： □2.35 , -4.12
5
w
Ew.d
実数データ 10の指数形式
w≥ d+7
E12.5: -0.23452E-01
12
Dw.d
倍精度実数データ 10の指数形式
Aw
文字データ 文字幅w 右詰め
A5 : abcde
wを省略：文字データの幅に合わせて書く
r …..
…. を r 回繰り返して書く
nX
n 桁スキップ
w≥ d+7
n=1の時でも省略できない
‘……’ 文字定数を書式の中に直接書く
“……” 同上
1x : □
3x : □ □ □
‘This is an example.’
‘(1x, I3,1x, “Hx=“, F5.3)’
Nagoya Institute of Technology
編集記述子 つづき
/
(
改行
)
書式項目のリストの一部分を囲む
反復などに利用
‘(/, 3(1x, e12.5))’
print *
write(*,*)
出力行の間に1行の空白行を入れる
read *
read ‘(A)’
文字列の時
Goｌgo13
‘AKB48’ と ‘ ‘ で囲んで入力する
直接文字列を入力すればよい
Nagoya Institute of Technology
!ex5-3
!----- Table of SIN and COS ----real, parameter :: rad=3.1415926/180.0
real :: r,s,c,p
integer ::k
print *,"kakudo
SIN
COS
CHECK“
do k=0,90,10
r=k*rad
s=sin(r)
c=cos(r)
p=s*s+c*c
print '(1x,i3,2x,3f10.5)',k,s,c,p
end do
end
Nagoya Institute of Technology
!ex5-4
!----- Mean and Standard Deviation ---integer :: n=0,mark
real :: w=0.0,v=0.0,mean,sd
do
print *,'score ( 0<x end)'
read *,mark
if (mark<0)exit
n=n+1
w=w+mark
v=v+mark**2
end do
if(n>0)then
mean =w/n; sd=sqrt(v/n-mean**2)
print '("average=",f5.1)',mean
print '("sd=",f5.1)',sd
end if
print '("total=",i5)',n
end
平均
標準偏差
Nagoya Institute of Technology
無限 do ループ
!ex5-5
!----
saidai kouyaku suu
---
integer :: m,n,k
print *,'Imput two integer m,n(>0):'
read *,m,n
if(m>0 .and. n>0)then
do
k=mod(m,n)
do ループから抜け出ることを
if (k==0) exit
確かに保証する
m=n
n=k
end do
if (n/=1)then
print *,'saidai kouyaku suu=',n
else
print *,'nothing'
end if
end if
end
Nagoya Institute of Technology
! --- Prime Number upto 100 --Integer :: n, m, q, max=1000
print *, 2
loop: do n = 3, max, 2
q=Nint(sqrt(real(n)))
do m = 3, q, 2
if( mod(n,m)==0 ) cycle loop
end do
print *, n
end do loop
end
! ２以外の偶数は調べる必要がない
! q以下のすべての奇数を調べて
! 割り切れたときには素数でないので
! ラベル loop のループへ出る
! ループを通れば素数である
! cycle loop で，ここへ来る
m
1 2 3
1
2
3
3からの奇数をqまで虱潰しに調べる． このqは以下のようにして
求めることができる．いま， qまで調べたとすると、このqの整数倍
は素数ではない． すなわち q*1,q*2,q*3,…, q*m は除外さ
れる．こうして調べられるmの最大はqまで．mがq以上ならmとqを
取り換えればよい．このq*qがn以上であればすべて網羅したこと
になるので、 q*q ≥ nと見積もれる．
3
n
q
n
Nagoya Institute of Technology
do while
[構文名] : do while ( 論理式 )
…………
end do [構文名]
do i=1,100
………
end do
同等
i=1
do while (i<=100)
………
i=i+1
end do
i=i+1 の位置に注意
i=1
do while (i<=100)
i=i+1
………
end do
i=101 でも実行される
Nagoya Institute of Technology
!ex 5.7
! Computation of sqrt(x) by recursive equation
implicit none
Real:: r=0.e0, x, d=1.e-5
Do
print *,'Input a positive real number'
read *, r
if(r>0.e0) exit
End do
x=r
do while (abs(x**2-r)/r >d)
x=0.5e0*(x+r/x)
end do
print '(1x,e12.5)’, x
end
収束条件：
Nagoya Institute of Technology
レポート課題２
テキスト p.56, 問題 5-2, 5-3, 5-6
自分でプログラムを作成し、これをレポートに書きだす
得られた計算結果をレポートに書き出す
表紙と考察を付けること
締切
： 次回の講義時 10:30 まで
提出先 ： ２号館Ｂ棟４階４２４Ｂ室前の
レポート回収箱
レポート課題番号
学籍番号
氏名
提出日
Nagoya Institute of Technology
数値積分
f(x)
fj
階段公式
h
x0 x1
台形公式
Simpson公式
xj
xn
x
Nagoya Institute of Technology
レポート課題３
テキスト p.57, 問題 5-5 の数値積分を以下に示す手順に沿って
行うこと
自分でプログラムを作成し、これをレポートに書きだす
得られた計算結果をレポートに書き出す
表紙と考察を付けること
締切
： 次回の講義時 10:30 まで
提出先 ： ２号館Ｂ棟４階４２４Ｂ室前の
レポート回収箱
レポート課題番号
学籍番号
氏名
提出日
Nagoya Institute of Technology
以下の積分を
（１） 解析的に実行せよ （計算式をきちんと書くこと）
（２） 数値積分を指定された2つの方法で求め、（１）で求めた結果と比較せよ
（３） 分点数を変えたとき数値計算結果がどのように変化するかを考察せよ
階段公式と台形公式
台形公式とシンプソン公式
台形公式とシンプソン公式
積分の上端が∞の時には、関数にもよるがある有限な値 Ｌ を仮定し、これを少しずつ
大きくしていくことで求める精度内収束するかどうかを確認して積分値を求める。
exp(-x*x/2)の場合 L=5 と取ると exp(-12.5)=3.72e-06だからこの程度のＬで十分であ
る。Ｌを変えて積分値がどのように収束していくかをグラフで見てみるのも面白いだろう。
Nagoya Institute of Technology
第5回
Nagoya Institute of Technology
配列
ベクトル、行列などの大量の変数を扱う
数値計算の要
配列宣言文 （1次元）
型, dimension(［下限: ]上限) :: 配列名のリスト
型:: 配列名(［下限: ]上限) ,配列名(［下限: ]上限) , ….
省略すると １
x,r は配列ではない（スカラー）
寸法＝上限ー下限＋１
配列要素 添え字で指定する
の
の 番目の要素
番目の要素
Nagoya Institute of Technology
配列の初期値設定
配列定数 (/定数, 定数, 定数, …/)
Do 型反復
配列構成子
(/(データのリスト, do 変数=初期値, 限界値, ［増分値］)/)
データのリスト
が配列構成子
実行文中でも使える
Nagoya Institute of Technology
配列の演算
配列に対する、代入、足し算、引き算、掛算を同時に行う
（１） 部分配列
配列名 ([添え字の開始値］： ［限界値］［：刻み幅］)
は要素数が１の1次元配列とみなす
ＯＫ
ＮＯ
はスカラー
は1次元配列
Nagoya Institute of Technology
（２） 配列代入文
配列＝同じ形状の他の配列
（形状：寸法と次元数）
すべての要素が０
奇数番が１
偶数番が０
Nagoya Institute of Technology
並び替え
!----Sorting Array Elements
----integer::m,mark(100),n=0,max,maxi,i,j
do while(n<=99)
print*,"Input integer,or negative one to stop:"
read*,m; if(m<=0) exit
n=n+1; mark(n)=m
end do
!----Sorting
do i=1,n-1
1 2
max=mark(i)
1
maxi=i
2
do j=i+1,n
if(mark(j)>max) then
max=mark(j)
i
maxi=j
end if
end do
if(maxi/=i) then
mark(maxi)=mark(i)
mark(i)=max
n
end if
end do
!----Output
do i=1,n
print"(' Rank',i3,' =',i10)",i,mark(i)
end do
end
n
i+1
最大値を探索
i
j
入れ替え
Nagoya Institute of Technology
（３） 配列演算式
配列名のままで要素間の演算を一括して書く
real::x(0:4)=(/(2.0**n,n=0,4)/)
x=x(0:4)*x(4:0:-1)
print*,x
end
Nagoya Institute of Technology
（４） ベクトル演算
配列に対する演算を同時に行う
と
は演算の内容が異なる！
A
time
a(0)
a(1)
a(2)
a(3)
a(4)
a(5)
a(6)
a(7)
a(8)
a(0)
a(1)
a(2)
a(3)
a(4)
a(5)
a(6)
a(7)
a(8)
a(9)
0
0
1
2
3
4
5
6
7
8
B
a(9)
time
a(0)
a(1)
a(2)
a(3)
a(4)
a(5)
a(6)
a(7)
a(8)
a(9)
0
0
0
0
0
0
0
0
0
0
Nagoya Institute of Technology
（５） 配列の入出力
配列名でまとめて出力するやり方と do 構文で書くやり方のちがいに注意
１行に10個
１行に1個を10行にわたって書く
Nagoya Institute of Technology
!----Pascal's Triangle
----integer, parameter::max=12
i
integer::i, c(0:max)=0
c(1)=1
0
do i=1,max
１
c(1:i)=c(1:i)+c(0:i-1)
print'(i2,":",50i5)',i,c(1:i) max １
end do
end program
１
１
１
１
２
３
４
１
３
６
0
１
４
１
Nagoya Institute of Technology
多次元配列
型, dimension(［下限: ]上限, ［下限: ]上限, …. ) :: 配列名のリスト
型:: 配列名(［下限: ]上限, ［下限: ]上限, …. ) , ….
または
または
または
多次元配列要素の順序と構成
a(1,1)
a(1,2)
a(2,1)
a(2,2)
a(3,1)
a(3,2)
メモリー上に記録される仕方
結果は同じだが右のほうが計算が速い
Nagoya Institute of Technology
配列の動的割り付け
実行時に必要な分だけ配列を割り当て、いらなくなったら開放する
型, allocatable :: 配列変数名(:, :, ….) のリスト または
型, dimension (:, :, …) , allocatable :: 配列変数名のリスト
……
allocate (配列変数名（添え字範囲の指定）, …. )
……
deallocate(配列変数名, …. )
Nagoya Institute of Technology
!----Matrix Product
----integer,dimension(:,:),allocatable::a,b,c
integer::m=0,n=0,k=0,i,j,h
!
do while(m<=0.or.n<=0.or.k<=0)
print*,"Input Array Size m, n, k ( > 0) :"
read*,m,n,k
end do
allocate(a(m,n),b(n,k),c(m,k))
!
print*,"Input Elements of A in a Matrix Form:"
read*,((a(i,j),j=1,n),i=1,m)
print*,"Input Elements of B in a Matrix Form:"
read*,((b(i,j),j=1,k),i=1,n)
!
do i=1,m
do j=1,k
c(i,j)=0
do h=1,n
c(i,j)=c(i,j)+a(i,h)*b(h,j)
end do
end do
end do
do i=1,m
print"(10i5)",(c(i,j),j=1,k)
end do
end
Nagoya Institute of Technology
where 文
配列のある条件を満たした要素の位置だけ、その配列に配列代入文を実行する
where （選別式）配列代入文
同時的代入
where 構文
［構文名：］ where （選別式）
配列代入文群
elsewhere （選別式）
配列代入文群
end where ［構文名］
Nagoya Institute of Technology
forall 文
forall （添え字範囲式，［選別式］）指標の書かれた配列代入文
指標変数＝［開始値］：［限界値］［：刻み値］
同時的代入
Nagoya Institute of Technology
forall 構文
［構文名：］ forall （添え字範囲式，［選別式］）
指標のついた配列代入文群
end forall ［構文名］
Nagoya Institute of Technology
Nagoya Institute of Technology
第６回
Nagoya Institute of Technology
文字列データ
人とコンピュータとのインターフェース、 データの効率的、機能的管理
文字型宣言
character ([len=]長さ) :: 変数名リスト
character :: 変数名１*長さ，変数名２*長さ，・・・・
character :: 配列名(［下限: ]上限) *長さ，・・・・
character (len=8) :: a,b,c=‘doraemon’
character :: fname(1:10)*8
文字部分列
文字列の1部分
文字変数名（［文字位置の下限］：［上限］）
character (len=8) ::
C(3:3)
C(4:8)
r
emon
c=‘doraemon’
1文字
Nagoya Institute of Technology
文字代入
文字は左詰めで代入
はみ出た分は切り捨て
character (len=8) ::
c
左辺の文字変数の長さ ＝ 右辺の文字列の長さ D o r
a
e
m o
n
＜
K o l
m o g o r
＞
D i
a
character (len=8) ::
a=‘abcdefgh’
a(1:5)=a(2:6)
‘bcdeffgh’
a(2:6)=a(1:5)
‘aabcdegh’
r
c
o v
Nagoya Institute of Technology
整数の１６進表示
!---- Hexadecimal Representation ---Implicit none
Integer :: dec, i
character(len=1) :: h(0:15)=(/'0','1','2','3','4','5','6','7' &
,'8','9','A','B','C','D','E','F'/)
character(len=8) :: hex
232=168=(24)8
ext: do
hex=" "
print*, "input a positive integer,or 0 to stop:"
read*, dec; if(dec<=0) exit ext
hxd:
do i=0,7
hex(8-i:8-i)=h(mod(dec,16))
dec=dec/16
if(dec==0) exit hxd
end do hxd
print ‘(1x,A8)’, hex
end do ext
end
右から埋めていく
次の桁へ行く
４ ０ Ｅ
Nagoya Institute of Technology
文字列の連結
文字列１ // 文字列２
character:: a*9,f*7='Fortran', g*2='95'
a=f//g
Print '(1x,a9)', a
Fortran95
f='boy'
g='s'
a=f//g
boy____s_
print'(1x,a9)', a
end
Character (1) :: f(9)=(/'1','2','3','4','5','6','7','8','9'/)
Character (22):: fmt
Integer :: n,k
n=2; k=3
Read *, int
fmt="( 2x, 'Int = ', "//f(n)//"I"//f(k)//" )"
Print fmt, int
end
Print ( 2x, ‘Int = ’, nIk ), int
と書いたことになる
Nagoya Institute of Technology
単語の並び替え
!---- rearrange characters of a word---character(len=4) :: w0,w1,w2,w(1:24)
integer :: i,j,k,n=1
print*,“input a word in 4 letters: "
read '(a)',w0
1番目に来る文字を先頭から順に選ぶ
do i=1,4
w1=w0(i:i)//w0(1:i-1)//w0(i+1:4)
b l u
２番目に来る文字を2番目以降から順に選ぶ
do j=2,4
w2=w1(1:1)//w1(j:j)//w1(2:j-1)//w1(j+1:4)
do k=3,4
w(n)=w2(1:2)//w2(k:k)//w2(3:k-1)//w2(k+1:4)
n=n+1
end do
end do
end do
print '(/12a5))', w
end
e
u l b e
Nagoya Institute of Technology
ファイル入出力
１．少量の計算データの格納と取り出し
コンソール入出力
（キーボードとディスプレイ）
print 書式識別子, 出力リスト
read 書式識別子, 出力リスト
f=‘( A10, F10.5 )’
write(6,f) x
write(*,‘( A10, F10.5 )’ ) x
print‘( A10, F10.5 )’, x
f=‘( 3I10 )’
read(5,f) i,j,k
read(*,‘( 3i10 )’) i,j,k
read(*,100) i,j,k
100 format(3i10)
同じ
同じ
Nagoya Institute of Technology
ファイル入出力
２．ハードディスク上での大量の計算データの格納と取り出し
入出力ファイルの割り当て
open 文, close 文
ハードディスクとの結合と分離
open（［unit=］ファイル装置番号, file=ファイル名 ［，制御項目］）
close（［unit=］ファイル装置番号, file=ファイル名 ［，制御項目］）
ファイル指定入出力
write（［unit=］装置識別子，［fmt=］書式識別子 ［，制御項目］）出力リスト
read（［unit=］装置識別子，［fmt=］書式識別子 ［，制御項目］）出力リスト
CPU
open(10,file=‘ data1’)
open(20,file=‘ data2’)
read(10) u1
write(20) u2
close(10)
close(20)
DISK
open(10,file=‘ data1’)
data1
u1
data2
u2
open(20,file=‘ data2’)
Nagoya Institute of Technology
open 文, close 文
open（［unit=］ファイル装置番号, file=ファイル名 ［，制御項目］）
close（［unit=］ファイル装置番号, file=ファイル名 ［，制御項目］）
装置識別子
•
•
•
unit を省略するときはファイル装置番号を先頭に書く
ファイル装置番号 は正の整定数または整数型変数（式）
ファイル名は 入力ファイル名を文字定数または文字変数で与える
制御項目
iostat= 整数変数 i
正常
i=0
エラー i > 0
status=‘文字列 ‘
old : 既存のファイルでなければならない
new : 同名のファイルが存在してなならない
replace : 同名のファイルあれば上書き
unknown : 入力ではnew, 出力ではreplace
delete : 終了後ファイルを消去（close 文のみ）
status= を省略すると unknown
position=‘文字列 ‘
rewind : 先頭位置から開始
append : 最終行のあとに追加
Nagoya Institute of Technology
ファイル指定入出力
write（［unit=］装置識別子，［fmt=］書式識別子 ［，制御項目］）出力リスト
read（［unit=］装置識別子，［fmt=］書式識別子 ［，制御項目］）出力リスト
装置識別子
*
番号
6
文字変数
コンソール入出力
open文で割り当てられるファイル装置番号
コンソール
内部ファイル
制御項目
iostat= 整数変数 i
正常
i=0
エラー
i>0
ファイルが終了 i < 0
advance=‘no’
入出力後に改行しない
Nagoya Institute of Technology
書式なし入出力文
書式識別子のない入出力文
バイナリ―データを入出力するので高速
大量のデータの入出力に使う
real*8, dimension(-10000:10000) :: a,b,c,d
read(10) a,b
write(20) c,d
内部ファイルの方法
character::d*10
Integer ::p=123456789
write(d,’(i10)’ ),p
d=‘_123456789’
dに数字列が右詰めで入る
Nagoya Institute of Technology
例１
Open (8,file=‘result’)
Write(8,’(8f10.7)’) (x(i),i=1,8)
close(8)
character(80) :: line(1:1000)
open(11,file=‘../mydata/report.txt’)
do j=1,1000
read(11,’(a)’,iostat=io) line(j)
if(io<0) exit
データの終わりに来たら do を出る
end do
close(11)
character:: filename*10
print*, ‘Input file name?’
read ‘(a)’, filename
open(10,file=filename)
Nagoya Institute of Technology
例２
iw.m : w桁の整数を右詰で書く
m桁より短い時には0を数字の前に書き加える
irun=32 のとき 0032
character(4) :: fn
. . . .
irun=32
write(fn,'(i4.4)') irun
内部ファイルの方法
open(10,file=‘int_'//fn//‘.dat')
Int_0032.dat というファイルが新たに作られる
open(20,file=‘spect.dat’)
. . . .
write(10,’(1x, 2(f7.4,1x))’) time,esum
. . . .
write(20) (energy(i), i=1,kmax)
close(10)
close(20)
spect.dat というファイルが新たに作られる
例３
read(8,’(i10)’,iostat=i) k
If(i>0) stop ‘Input file error’
Int_0032.dat に time, esum のデータを書き出す
spect.dat にenergy のデータを書き出す
Nagoya Institute of Technology
その他 ファイル入出力制御文
１つ前の記録の先頭に戻る
backspace（［unit=］ファイル装置番号, isostat=整数変数）
ファイルの先頭に戻る
rewind（［unit=］ファイル装置番号, isostat=整数変数）
ファイル終了記録を書き込む
endfile（［unit=］ファイル装置番号, isostat=整数変数）
Nagoya Institute of Technology
!------ Character Code Table -----integer ::i,j
open(8,file='Code_table')
do i=0,31
write(8,"(7(2x,i3,2x,a))")(i+32*j,char(i+32*j),j=1,7)
end do
close(8)
end
Nagoya Institute of Technology
第７回
Nagoya Institute of Technology
副プログラム
構成的プログラムをつくる
主プログラム
副プログラム
Subroutine
あるまとまった計算・作業のみを行う
…..
call sub_A
call sub_B
…..
call sub_C
….
call sub_A
….
end
sub_A
よくできたサブルーチンは他の
プログラムでも使える
sub_B
sub_C
call sub_G
sub_G
routine
決まりきった仕事, 日常の仕事
Nagoya Institute of Technology
subroutine example
３つの数字を読んで、大きい順に並び替える
元プログラム
サブルーチン
implicit none
integer :: i,j,k
read *,i,j,k
if(i<j) call swap(i,j)
if(j<k) call swap(j,k)
if(i<j) call swap(i,j)
print *,i,j,k
end
!
!
!
subroutine swap(p,q)
integer :: p,q,r
r=p; p=q; q=r
end
実引数
（実際に出入りする変数）
仮引数
（仮に与えられた変数）
Nagoya Institute of Technology
subroutine文
subroutine サブルーチン名［（仮引数のリスト）］
………..
end subroutine ［サブルーチン名］
配列も許される
仮引数のないものもある
call 文
call サブルーチン名［（仮引数のリスト）］
return 文
return
元プログラムに戻る位置を設定する
………
call sub(…,code)
select case(code)
case(1)
……
case(2)
……
default
……
end select
……
subroutine sub(…,return_code)
……
return_code=1
return
……
return_code=2
return
……
return_code=999
end subroutine sub
Nagoya Institute of Technology
副プログラム
・ 主プログラムはただ一つ
・ 副プログラムはいくつあってもよい
・ 副プログラムから他の副プログラムをよんでもよい
・ 副プログラムが自分自身を呼ぶことはできない（再起呼び出し指定があれば別）
・ 実引数と仮引数の型、個数、形状は一致していなければいけない
Nagoya Institute of Technology
実引数と仮引数の型、個数、形状は一致していなければいけない
integer, parameter::n=100
real, dimension (0:n):: x,y
………
call set_initial(x,y,n)
………
End
!
subroutine set_initial(a,b,m)
integer :: m,i
real, dimension(0:m)::a,b
do i=0,m
a(i)=0.0
b(i)=0.0
end do
end set_initial
実引数の個数、型、形状
仮引数の個数、型、形状
変数の順序に注意！
並んだ順に渡される
call set_initial(y,x,n)
subroutine set_initial(a,b,m)
・変数の個数、型、形状が一致していること
・変数の順序に注意
・変数間の対応が取れていればよく、同じ変数名である必要はない
Nagoya Institute of Technology
変数の有効範囲
変数名は１つのプログラム単位内でのみ有効
文番号、構文名も１つのプログラム単位内でのみ有効
主プログラム
…..
call sub_A
….
y=sin(2.0*x)
….
print *, x,y
end
副プログラム
変数 x,y は
共有しない
Subroutine sub_A
….
x=2.0*k
y=exp(-x)
…..
end subroutine sub_A
どうやって変数を共有するのか？
文番号、構文名は共有する方法がない （なくても困らない）
Nagoya Institute of Technology
変数の共有方法
１ 引数として共有する
２ common文を用いる
３ module文を用いる
１．引数として共有
call set_initial(x,y,z,n)
(x,y,z,n)
(a,b,c,m)
2． common文
subroutine set_initial(a,b,c,m)
・変数の個数、型、形状が一致していること
・変数の順序に注意
・変数間の対応が取れていればよく、同じ変数名である必要はない
common［/共通ブロック名/］変数のリスト
プログラム１
common/param1/x,y,z
Common/param2/a,pi
….
プログラム２
common/param1/p,q,r
….
(x,y,z)
(p,q,r)
プログラム３
common/param2/s,pi
….
(a,pi)
(s,pi)
・変数の個数、型、形状が一致していること
・変数の順序に注意
・変数間の対応が取れていればよく、同じ変数名である必要はない
Nagoya Institute of Technology
save 宣言
副プログラムが1度呼び出されたとき定義された変数の値は、次回にその副プログラム
が呼び出されたとき、その値を保持しているかどうかはわからない
確実に保証する場合に save 文を用いる
主プログラム
integer :: n
real :: z
do n=1,100
call sub_A(z,n)
……
end do
end
副プログラム
subroutine sub_A(w,k)
real, save :: pi
……
if(k==1) then
……
pi=3.14159265
save宣言がないと２回目以降
end if
piの値が確定しているかどう
w=cos(0.1*pi*k)
か保障されない
……
end subroutine sub_A
Nagoya Institute of Technology
仮引数の授受属性
引き渡される変数が、入力（in）、出力（out）、入出力（inout）のいずれかを特定
間違い防止に効果がある
intent 宣言
型, intent（in, out, inoutのいずれか）：：変数名のリスト
implicit none
integer :: i,j,k
read *,i,j,k
if(i<j) call swap(i,j)
if(j<k) call swap(j,k)
if(i<j) call swap(i,j)
print *,i,j,k
end
subroutine swap(p,q)
implicit none
integer, intent(inout) ::p,q
integer :: r
r=p; p=q; q=r
end
Nagoya Institute of Technology
例題１
!------Drawing a Histogram -------integer :: nn(0:9),j
character(len=80) :: cn
character(len=30) :: fmt
nn=(/1,5,12,22,28,35,48,28,10,3/)
fmt='(1x,i2,"-",i2,1x,a50)'
do j=0,9
call str(nn(j),cn)
50-59 I############
という具合に書き出す
print fmt,10*j,10*j+9,cn
end do
end
!
subroutine str(n,c)
integer, intent(in) :: n
character(len=80), intent(out) :: c
integer ::k
コメント文にしてみる ! c=‘ ‘
c=‘ ‘
c(1:1)='I'
do k=2,n+1
I######### を書く
c(k:k)='#'
end do
end
Nagoya Institute of Technology
例題２
subroutine count(n)
implicit none
integer, intent(out) :: n
integer :: i, sum=0
! sum=0
do i=1,10
sum=sum+i
end do
n=sum
sum に初期値を設定 （一番最初だけ）
毎回 sum=0 にする
end subroutine count
!
Program main
implicit none
integer :: n
call count(n)
print '(1x,i6)',n
n=55
call count(n)
print '(1x,i6)',n
n=110
end
前のsum=55がそのまま残っているから
Nagoya Institute of Technology
外部関数
プログラマが副プログラムを用いて自分で定義する関数
外部手続 ： サブルーチン、外部関数
組み込み手続き ： Fortran システムであらかじめ用意されている数学関数
関数副プログラム
［型］ function 関数名（仮引数のリスト）［result（変数名）］
………
end function 関数名
function 関数名（仮引数のリスト）［result（変数名）］
型 :: 関数名 または result 変数名
………
end function 関数名
関数の値の代入の仕方
関数名 または result 変数名= 関数の値を与える式
Nagoya Institute of Technology
関数の値の代入
関数名 または result 変数名= 関数の値を与える式
外部関数の引用
関数名（実引数のリスト）
integer :: n=100
real :: x,dx,sum=0.0,a=1.0
dx=1.0/n
do i=1,n-1
x=i*dx
sum=sum+fr(x,a)
end do
台形公式
sum=(sum+0.5*(fr(0.0,a)+fr(1.0,a)))*dx
print ‘(1x,e10.3)’,sum
end
!
real function fr(x,p)
・変数の個数、型、形状が一致していること
・変数の順序に注意
real :: x,z,p
・変数間の対応が取れていればよく、同じ変数名である必要はない
z=(p-x*x)/(p+x*x)
fr=log(z)
end function fr
Nagoya Institute of Technology
例題
１６進表現
!-- Generating a Hexadecimal Representation ------function hex(n)
character(len=8) :: hex
character(len=1) :: h(0:15)=(/ '0','1','2','3','4','5','6','7',&
'8','9','A','B','C','D','E','F'/)
integer, intent(in) ::n
n は入力
integer :: nin,j,nn
hex=' '
nin=n
do j=8,1,-1
nn=nin/16
hex(j:j)=h(nin-16*nn)
hex は文字型変数で出力
右から埋めていく
if(nn==0) exit
nin=nn
end do
end function hex
!---- Main ---program main
implicit none
character(len=8) :: hex
integer :: i
print *, 'Input a positive Integer or negative one to stop:'
do
read *, i; if(i<0) exit
print *,hex(i)
end do
end
Nagoya Institute of Technology
関数属性宣言
使用する関数の名前が外部関数名（組み込み関数名）であることを
コンパイラに知らせる
external 宣言
型, external :: 関数名のリスト
external :: 関数名のリスト
intrinsic 宣言
外部関数名がたまたま組み込み関数名と同じになったとき
外部関数の定義が優先される
型, intrinsic :: 関数名のリスト
intrinsic :: 関数名のリスト
call integral(gamma,0.0,1.0,s)
…………
real function gamma(x)
…………
関数が引数として入るときには属性宣言が必要
Nagoya Institute of Technology
例題
implicit none
real:: h, phix,err
real, external ::Gaussian
print *, 'What is your deviation value?'
read *, h
real function Gaussian(s)
implicit none
real, parameter :: rn=2.506628e0
real, intent (in) :: s
Gaussian=exp(-s*s/2.0e0)/rn
end function Gaussian
h=(h-50.0e0)/10.e0
call simpson(Gaussian, 0.0e0,h, 16,phix)
err=(0.5e0-phix)*100.0e0
print '(1x,a,f6.2,a)', &
'From the top you are ', err, ' %.'
end
subroutine simpson(func,a,b,kubun,integ)
implicit none
integer, intent(in)::kubun
integer:: i
real, external :: func
real, intent(in) :: a,b
real, intent(out) :: integ
real :: x,dx
x
dx=(b-a)/(2.0e0*kubun)
integ=func(a)+4.0e0*func(a+dx)+func(b)
do i=1,kubun-1
x=a+(2.0e0*i)*dx
integ=integ+2.0e0*func(x)+4.0e0*func(x+dx)
end do
integ=integ*dx/3.0e0
end subroutine simpson
Simpsonの公式
! \sqrt(2*pi)
Nagoya Institute of Technology
内部手続き
サブルーチンや関数を1つのプログラム単位の中に含めてしまうこと
integer :: n
real :: z,g,pi
pi=3.14159265
do n=1,100
call sub_A(z,n)
……
g=func(z)
……
end do
………
Contains
subroutine sub_A(x,m)
………
z=log(func(x))
………
end subroutine sub_A
!
real function func(x)
………
func=0.5*exp(-pi***0.25d0)
end function func
end
contains文
サブルーチン、関数の規則が当てはまる
１．親プログラム中で、内部手続き関数名に対する型宣言および
関数属性宣言は不要
real, external :: func
エラーとなる
２．親プログラム中で使われている変数名、配列名は内部手続き中でも
有効．親プログラムと同じ値を共有し内部手続き中で宣言する必要なし
pi=3.14159265
３．親プログラム中で使われている変数名と同じ名前を内部手続き中で
型宣言すると、その変数は内部手続き中でのみ有効な局所変数と
みなされ、親プログラムと共通の値を持たなくなる
４．内部手続き中に contains 文を書くことはできない
５．内部手続きから同じ親プログラム単位中の別の内部手続きを引用する
ことは可能
func
６．内部手続きの end 文では、後ろの手続き種別（subroutine, function）
などを省略できない（省略すると親プログラムの end 文と区別がつかな
くなる）
Nagoya Institute of Technology
例題
!----- All Round Calendar -----program main
character(len=9) :: name(0:6)=(/'Sunday
', &
'Monday
', 'Tuesday ','Wednesday', &
'Thursday ','Friday
','Saturday '/)
integer :: md(1:12)=(/31,28,31,30,31,30,31,31,&
30,31,30,31/)
integer :: y,m,d,ld
do
print *,'year?'
read *,y; if(y<=0) exit
md(2)=28+leap(y)
print *,'month?‘
do
read *,m; if(0<m .and. m<=12) exit
print *, 'm<=12'
end do
print *, 'Day?'
read *, d
if(d<=0 .or. d>md(m)) then
print *, 'nothing'
else
call count(ld)
print "(1x,i5,',',i2,',',i2,',',a9,a/)", &
y,m,d,name(ld)
end if
end do
!
contains
!
西暦１年１月１日から何日目かを数えて曜日を算出
subroutine count(n)
integer, intent(out) :: n
n=yday(y)+sum(md(1:m-1))+d
n=mod(n,7)
end subroutine
!
!
前年の１２月３１日が西暦１年１月１日から何日目か
function yday(x) result(days)
integer :: x,days
days=365*(x-1)+(x-1)/4-(x-1)/100+(x-1)/400
end function yday
!
integer function leap(x) 閏年なら１、平年なら０
integer :: x
leap=yday(x+1)-yday(x)-365
end function leap
!
end program
Nagoya Institute of Technology
配列の共有
数値計算では非常によく出てくる
親プログラム中の配列をその形状や寸法に応じて副プログラムに渡す
１．形状明示配列
添え字の上下限が明示されている配列 （静的）
２．形状引き継ぎ仮配列
副プログラムの宣言文で寸法は指定しないでおき、引用した段階で送られてく
る実引数配列の寸法に連動して割り付ける配列 （動的）
Nagoya Institute of Technology
形状明示配列
添え字の上下限が明示されている配列 （静的）
添え字の上下限（寸法）の渡し方
１．仮引数になっていて元プログラムから値が送られてくる整数変数
２．仮引数となっている配列から size 関数によって求められる寸法の値
３．Common文により元プログラムと値を共有する整数変数
４．内部手続きにおいては、親プログラム中の整数変数
Subroutine sub(x,n)
integer, intent(in):: n
integer :: a(1:n),j
real, intent(out) :: x
real, allocatable ::temp(:)
allocate(temp(0:n))
a(1:n),temp(:)
仮引数ではなく、サブルーチン内だけで使用される配列
Nagoya Institute of Technology
例題
応募者 ｎ 人の中から抽選で ｍ 人を選ぶ
program main
integer, allocatable :: number(:)
integer ::n,m,i
print *,‘応募者数は？’
read*,n
応募者数
print *, ‘当選者数は?’
read*,m
当選者数
allocate(number(n))
call chusen(number,n,m)
!
print *,‘順位 当選番号’
do i=1,min(n,m)
print '(2x,i4,i8)',i,number(i)
end do
end
subroutine ran(i,r)
(0,1)の一様疑似乱数生成
integer, intent(inout):: i
real, intent(out) :: r
integer, parameter ::mask=2147483647,a=48828125
i=iand(a*i,mask)
r=real(i)/real(mask)
end subroutine
subroutine chusen(kk,n1,n2)
integer, intent(in) :: n1,n2
integer, intent(out) :: kk(n1)
integer :: i,ir,j,temp
real ::x
kk=(/(i,i=1,n1)/)
print *,'ir'
read *, ir
do i=1,min(n1,n2)-1
call ran(ir,x)
j=i+int(x*(n1-i+1))
temp=kk(i)
kk(i)=kk(j)
kk(j)=temp
end do
end subroutine
乱数の種（任意の正の整数）を入れる
（0,1）の乱数
残りから選ぶ
i,j を入れ替える
Nagoya Institute of Technology
形状引き継ぎ仮配列
副プログラムの宣言文で寸法は指定しないでおき、引用した段階で送られてくる
実引数配列の寸法に連動して割り付ける配列 （動的）
親プログラム
副プログラム
integer, parameter:: n=100
real :: u(0:n,0:n),x(0:n),y(0:n)
interface
subroutine fmake(w,a,b)
real, intent(in)::w(0:,0:),a(0:)
real, intent(out)::b(0:)
end subroutine fmake
end interface
…………
call fmake(u,x,y)
…………
end
実引数の実際の大きさに
合わせて割り付けられる
Interface 構文でコンパイラに予め副プログラムの
構造を知らせておく
subroutine fmake(w,a,b)
real, intent(in)::w(0:,0:),a(0:)
real, intent(out)::b(0:)
integer ::m
Real :: sum
m=ubound(w,1)
W の第1次元方向の上限値
do i=0,m
sum=0.0
do j=0,m
sum=sum+w(i,j)*a(j)
end do
b(i)=sum
end do
end subroutine fmake
Nagoya Institute of Technology
interface 構文
Interface
引用する手続き名
その中の仮引数の形状、属性など
end interface
例1
Subroutine test(x)
implicit none
real :: x(:)
! real :: x(0:)
print ‘(2i3)’, lbound(x), ubound(x)
end subroutine test
Program main
implicit none
integer :: I
real :: s(-10:10)=(/(i,i=1,21)/)
interface
subroutine test(x)
real :: x(:)
end subroutine test
end interface
call test(s)
call test(s(0:8))
End
1 21 を返す
1
9 を返す
Nagoya Institute of Technology
例２
正方行列のトレース
program main
integer :: n,i
real, allocatable :: a(:,:)
interface
function trace(x) result(tr)
real, intent(in) :: x(1:,1:)
real :: tr
end function
end interface
print *,'n?'
read *,n
allocate(a(1:n,1:n))
do i=1,n
! print '(1x,a,i2)','line',i
do j=1,n
print '(1x,a,i2,a,i2,a)', &
'Enter element a(',i,' ,',j,')'
read *,a(i,j)
end do
end do
print '(1x,a,F7.4)','trace=',trace(a)
end
!----Trace of a Square Matrix
function trace(x) result(tr)
Implicit none
integer :: n1,n2,i
real, intent(in) :: x(1:,1:)
real :: tr
n1=size(x,1)
n2=size(x,2)
! テキスト p.82
if(n1/=n2) then
print *, 'not square';return
else
tr=0.0
do i=1,n1
tr=tr+x(i,i)
end do
end if
end function
Nagoya Institute of Technology
内部手続き中の配列
親プログラム中で定義されている整数変数を用いて配列の寸法を決めることができる
例 スカラー３重積
ベクトル３重積
implicit none
real, dimension(1:3) :: a,b,c,bc
real, external ::
vector_prod, scalar_prod
はいらない
print *, "Input three vector's components"
read *, a
read *, b
read *, c
bc = vector_prod(b,c)
print '( "Sp3 = ",F10.3)', scalar_prod(a,bc)
print '( "Vp3 = (",3F10.3,")")', vector_prod(a,bc)
b x c
a . (b x c)
a x (b x c)
contains
内部手続き 復習
function vector_prod(x,y) result(z)
implicit none
real, dimension(1:3), intent(in) :: x,y
real, dimension(1:3) :: z
intent はいらない
z(1)=x(2)*y(3)-x(3)*y(2)
z(2)=x(3)*y(1)-x(1)*y(3)
z(3)=x(1)*y(2)-x(2)*y(1)
end function vector_prod
function scalar_prod(x,y) result(sp)
implicit none
real, dimension(1:3), intent(in) :: x,y
real :: sp
sp=x(1)*y(1)+x(2)*y(2)+x(3)*y(3)
end function scalar_prod
end
１．親プログラム中で、内部手続き関数名に対する型宣言および 関数属性
宣言は不要
２．親プログラム中で使われている変数名、配列名は内部手続き中でも
有効．親プログラムと同じ値を共有し内部手続き中で宣言する必要なし
３．親プログラム中で使われている変数名と同じ名前を内部手続き中で
型宣言すると、その変数は内部手続き中でのみ有効な局所変数と
みなされ、親プログラムと共通の値を持たなくなる
６．内部手続きの end 文では、後ろの手続き種別（subroutine, function）
などを省略できない
Nagoya Institute of Technology
再帰的呼び出し
副プログラムが自分自身を呼び出す
無限ループの罠にはまらない！
再帰的手続きの指定
recursive subroutine サブルーチン名（仮引数のリスト）
recursive function 関数名（仮引数のリスト） result（変数名）
必須
recursive : 再帰的 繰り返して用いられる
Nagoya Institute of Technology
ルジャンドル多項式
n=0
n=1
n=2
n=3
n=4
チェビシェフ多項式
n=0
n=1
n=2
n=3
n=4
エルミート多項式
n=4
n=2
n=1
n=0
n=3
Nagoya Institute of Technology
例 ルジャンドル多項式
偏微分方程式の解などによく出てくる代表的な特殊関数
漸化式
program main
implicit none
integer :: n,i
real :: x,y,p
character(len=1) :: c(-25:25)
print *,'What is the order of Polynomial, n?'
read *,n
print '( "Legendre Polynomial P(",i2,",x)")',n
do i=-20,20
x=0.05*i
! dx=1/n=0.05
y=p(n,x)
call presentation
print '(f5.2,2x,f10.6,1x,51a1)',x,y,c
end do
contains
subroutine presentation
integer ::j
if(i/=0) then
c=' '
else
c=(/('-',j=-25,25)/)
end if
c(-25)='I';c(0)='I';c(25)='I'
c(nint(25*y))='*'
end subroutine
end
recursive function p(n,x) result(pnx)
real :: pnx
integer, intent(in) ::n
real,intent(in) ::x
if(n==0) then
pnx=1.0
else if(n==1) then
pnx=x
else
pnx=(2.0-1.0/n)*x*p(n-1,x) &
-(1.0-1.0/n)*p(n-2,x)
end if
end function
Nagoya Institute of Technology
Legendre, Chebyshev, Hermite polynomials 参考例
! Computation of special functions
integer, parameter:: fun_kind=3 ! Choose type of special function
! fun_kind=1: Legendre, =2: Chebyshev, =3: Hermite
real*8, parameter:: x_max=1.d0
Integer, parameter:: n_grid=50, n_order=4, k_write=2
character(len=1):: cn(0:9)=(/'0','1','2','3','4','5','6','7','8','9'/)
character(len=2):: fn(1:3)=(/'p_','t_','h_'/)
real*8 :: dx,x,y,fun
integer :: i,n,n_n,n_p,n_extend
select case (fun_kind)
case(1) ; print *,'
case(2) ; print *,'
case(3) ; print *,'
end select
Legendre Polynomial P_n(x) '
Chebyshev Polynomial T_n(x) '
Hermite Polynomial H_n(x) '
dx=x_max/n_grid
n_extend=1
if(fun_kind==3) n_extend=4
n_n=-n_grid*n_extend
n_p= n_grid*n_extend
! extend the x-domain for Hermite polynomial
do n=0,n_order
write(6,'(/," n=",i2,/)') n
open(10+n,file=fn(fun_kind)//cn(n)//'.dat')
do i=n_n,n_p
x=i*dx
y=fun(fun_kind,n,x)
write(10+n,'(1pe10.3, 2x, 1pe10.3)') x, y
if(mod(i,k_write)==0) write(6,'(1pe10.3, 2x, 1pe10.3)') x,y
end do
end do
end
1/2
Nagoya Institute of Technology
recursive real*8 function fun(kind,n,x) result (sf)
integer, intent(in) :: kind, n
real*8, intent(in) :: x
if(n==0) then
sf=1.d0
else if(n==1) then
select case(kind)
case(1); sf=x
case(2); sf=x
case(3); sf=2.d0*x
end select
else
select case(kind)
case(1); sf=(2.d0-1.d0/n)*x*fun(kind,n-1,x)-(1-1.d0/n)*fun(kind,n-2,x)
case(2); sf=2.d0*x*fun(kind,n-1,x)-fun(kind,n-2,x)
case(3); sf=2.d0*x*fun(kind,n-1,x)-2.d0*(n-1)*fun(kind,n-2,x)
end select
end if
end
2/2
Nagoya Institute of Technology
レポート課題４
サブルーチン、外部関数、再帰的呼び出しの方法を用いて次の積分をそれぞれ指定された
積分法で求めよ
自分でプログラムを作成し、これをレポートに書きだす
得られた数値計算結果をレポートに書き出す また解析計算と比較する
表紙と考察を付けること
締切
： 6月6日（月）10:30 まで
提出先 ： ２号館Ｂ棟４階４２４Ｂ室前の
レポート回収箱
（１）
Simpson 公式
（２）
台形公式
レポート課題番号
学籍番号
氏名
提出日
Nagoya Institute of Technology
課題 参考例
!
Integration of special functions
implicit none
real:: xmin,xmax,int
real, external :: func
! kind=1: Legendre, =2: Tshebysheff, =3: Hermite
integer :: kind,n,scheme,nmax,i
real, parameter :: pi=3.1415926535897, s=sqrt(pi/2.0e0)
real, dimension(3) :: P_exact(3)=(/1.0,0.0,0.0/), H_exact(3)=(/s,2.e0*s,12.e0*s/)
print *
print *, '
*** Integral of Legendre Polynomials on [0,1] *** '
kind=2; scheme=2; xmin=0.e0; xmax=1.0e0; nmax=200
do i=0,2
n=2*i
call integral(func,kind,n,xmin,xmax,scheme,nmax,Int)
print '(/," n=",i1,1x,", I=",1pe12.5,", I_exact=",1pe12.5 )', n,Int,P_exact(i+1)
end do
print *
print *, '
*** Integral of Hermite Polynomials times exp(-x*x/2) on [0,infinity] *** '
kind=3; scheme=1; xmin=0.e0; xmax=5.e0; nmax=200
do i=0,2
n=2*i
call integral(func,kind,n,xmin,xmax,scheme,nmax,Int)
print '(/," n=",i1,1x,", I=",1pe12.5,", I_exact=",1pe12.5)', n,Int, H_exact(i+1)
end do
end
Nagoya Institute of Technology
subroutine integral(fun,kind,n,a,b,scheme,kubun,integ)
implicit none
integer, intent(in):: kind,n,scheme,kubun
integer:: i
real, external :: fun
real, intent(in) :: a,b
real, intent(out) :: integ
real :: x,dx
if(scheme==1) then
! Trapezoidal rule
dx=(b-a)/kubun
integ=0.5*(fun(kind,n,a)+fun(kind,n,b))
do i=1,kubun-1
x=a+i*dx
integ=integ+fun(kind,n,x)
end do
integ=integ*dx
else
! Simpson rule
dx=(b-a)/(2.0e0*kubun)
integ=fun(kind,n,a)+4.0e0**fun(kind,n,x+dx)+fun(kind,n,b)
do i=1,kubun-1
x=a+(2.0e0*i)*dx
integ=integ+2.0e0*fun(kind,n,x)+4.0e0*fun(kind,n,x+dx)
end do
integ=integ*dx/3.0e0
endif
end subroutine integral
Nagoya Institute of Technology
real function func(kind,n,x)
integer, intent(in) :: kind, n
real, intent(in) :: x
real, external :: sp_fun
func=sp_fun(kind,n,x)
if(kind==3) func=sp_fun(kind,n,x)*exp(-0.5e0*x*x)
end
recursive real function sp_fun(kind,n,x) result (sf)
integer, intent(in) :: kind, n
real, intent(in) :: x
if(n==0) then
sf=1.e0
else if(n==1) then
select case(kind)
case(1); sf=x
case(2); sf=x
case(3); sf=2.e0*x
end select
else
select case(kind)
case(1); sf=(2.e0-1.e0/n)*x*sp_fun(kind,n-1,x)-(1-1.e0/n)*sp_fun(kind,n2,x)
case(2); sf=2.e0*x*sp_fun(kind,n-1,x)-sp_fun(kind,n-2,x)
case(3); sf=2.e0*x*sp_fun(kind,n-1,x)-2.e0*(n-1)*sp_fun(kind,n-2,x)
end select
end if
end
Nagoya Institute of Technology
Nagoya Institute of Technology
Nagoya Institute of Technology
Nagoya Institute of Technology
例題
implicit none
real:: h, phix,err
real, external ::Gaussian
print *, 'What is your deviation value?'
read *, h
real function Gaussian(s)
implicit none
real, parameter :: rn=2.506628e0
real, intent (in) :: s
Gaussian=exp(-s*s/2.0e0)/rn
end function Gaussian
h=(h-50.0e0)/10.e0
call simpson(Gaussian, 0.0e0,h, 16,phix)
err=(0.5e0-phix)*100.0e0
print '(1x,a,f6.2,a)', &
'From the top you are ', err, ' %.'
end
subroutine simpson(func,a,b,kubun,integ)
implicit none
integer, intent(in)::kubun
integer:: i
real, external :: func
real, intent(in) :: a,b
real, intent(out) :: integ
real :: x,dx
dx=(b-a)/(2.0e0*kubun)
integ=func(a)+4.0e0**func(x+dx)+func(b)
do i=1,kubun-1
x=a+(2.0e0*i)*dx
integ=integ+2.0e0*func(x)+4.0e0*func(x+dx)
end do
integ=integ*dx/3.0e0
end subroutine simpson
Simpsonの公式
! \sqrt(2*pi)
Nagoya Institute of Technology