ナビゲーション：前へ 上へ 次へ
3 区分求積法を用いて定積分を行うプログラム
ここでは関数 1/x**2 の区間 [1, 10] 上の定積分を行うプログラムを取り上げます。
計算を行い、その結果と、計算にかかった時間を画面に出力します。
3.1 逐次計算版のプログラムコード
最初に共配列を用いない逐次計算プログラムの例を示します。本プログラムは以下の２つのソースファイルで構成されます。
- serial.f90 - 主プログラム、被積分関数
- quad1.f90 - 区分求積法を行うサブルーチンを含むモジュール
1: Program example
2: Use quad1
3: Implicit None
4: Integer,Parameter :: wp = Selected_Real_Kind(15)
5: Real(wp) res,true
6: Double Precision secs
7: Real(wp),Parameter :: a = 1.0_wp, b = 10.0_wp
8: Integer,Parameter :: steps = 1000*1000*1000
9: Integer(Selected_Int_Kind(18)) t1,t2,irate
10: Call System_Clock(Count_Rate=irate)
11: Call System_Clock(t1)
12: Call rectangle_rule(f,a,b,steps,res)
13: Call System_Clock(t2)
14: Print *,'Calculated value: ',res
15: true = 1 - 0.1d0
16: Print *,'True value (approx):',true
17: Print *,'Relative error',Abs((res-true)/true)
18: secs = (t2-t1)/Dble(irate)
19: If (secs>=0.1d0) Then
20: Print 1,secs,'seconds'
21: Else
22: Print 1,secs*1000,'milliseconds'
23: End If
24: 1 Format(1X,'Time taken ',F0.3,1X,A)
25: Contains
26: Real(wp) Function f(x)
27: Real(wp),Intent(In) :: x
28: f = 1/x**2
29: End Function
30: End Program
------------------
10行目：システムクロックが１秒間に刻むカウント数を得ます。
11行目：計算開始時のクロックカウントを得ます。
12行目：区分求積分を行う rectangle_rule という関数を呼び出します。
（rectangle_rule 関数は quad1.f90 に定義されています。）
13行目：計算終了時のクロックカウントを得ます。
14-23行目：計算結果及び経過時間を出力します。
26-29行目：被積分関数 1/x**2
[
quad1.f90
] - 区分求積を行うrectangle_ruleを含むモジュール
1: Module quad1
2: Implicit None
3: Private
4: Public rectangle_rule
5: Interface rectangle_rule
6: Module Procedure rectangle_dp
7: End Interface
8: Contains
9: Subroutine rectangle_dp(f,from,to,nstep,integral)
10: Interface
11: Double Precision Function f(x)
12: Double Precision,Intent(In) :: x
13: End Function
14: End Interface
15: Double Precision,Intent(In) :: from,to
16: Integer,Intent(In) :: nstep
17: Double Precision,Intent(Out) :: integral
18: Double Precision x1,x2,y,stepval
19: Integer i
20: If (to<=from) Error Stop 'Invalid from-to'
21: stepval = (to-from)/nstep
22: x1 = from
23: integral = 0
24: Do i=1,nstep
25: x2 = from + i*stepval
26: y = f((x1+x2)/2)
27: integral = integral + y
28: x1 = x2
29: End Do
30: integral = integral*stepval
31: End Subroutine
32: End Module
------------------
4-7行目： rectangle_ruleという名前の総称引用仕様名を定義します。
今回のコードでは、以下に示すようにこの総称引用仕様名は一つの個別手続きを持ちます。
9-31行目：上記総称引用仕様名rectangle_ruleの個別手続きrectangle_dpを定義します。
この手続きは倍精度の引数を取りますが、必要に応じて例えば単精度の引数を取る
別の個別手続きrectangle_spというようなサブルーチンを後々追加するような事も可能です。
このサブルーチンで実際に区分求積法を用いて、被積分関数fについて区間[from,to]の定積分を計算します。
（計算において区間[from,to]はnstep等分されます）
3.1.1 出力結果例
本プログラムを実行すると例えば以下のような出力が得られます。（実行環境により得られる数値は異なります）
Calculated value: 0.9000000000003117 True value (approx): 0.9000000000000000 Relative error 3.4626622556920158E-13 Time taken 4.950 seconds
各数値は以下の通りです。
- Calculated value - 計算によって得られた結果
- True value (approx) - 解析的に得られる正しい結果
- Relative error - 相対誤差
- Time taken - 計算に要した時間（秒）
3.2 逐次計算プログラムを共配列プログラムに改変する
このセクションでは、上述の区分求積分プログラムを変更して、共配列プログラムに仕立て上げます。
3.2.1 主プログラムの変更
以下に変更後の主プログラムを示します。
[
coarrays.f90
] - 共配列プログラム（変更適用済み）の主プログラム
1: Program coexample 2: Use coquad1 3: Implicit None 4: Integer,Parameter :: wp = Selected_Real_Kind(15) 5: Real(wp) res[*],true 6: Double Precision secs 7: Real(wp),Parameter :: a = 1.0_wp, b = 10.0_wp 8: Integer,Parameter :: steps = 1000*1000*1000 9: Integer(Selected_Int_Kind(18)) t1,t2,irate 10: Call System_Clock(Count_Rate=irate) 11: Call System_Clock(t1) 12: Call co_rectangle_rule(f,a,b,steps,res) 13: Call System_Clock(t2) 14: If (This_Image()==1) Then 15: Print *,'Calculated value: ',res 16: true = 1 - 0.1d0 17: Print *,'True value (approx):',true 18: Print *,'Relative error',Abs((res-true)/true) 19: secs = (t2-t1)/Dble(irate) 20: If (secs>=0.1d0) Then 21: Print 1,secs,'seconds',Num_Images(),Num_Images()*secs 22: Else 23: Print 1,secs*1000,'milliseconds',Num_Images(),Num_Images()*secs 24: End If 25: 1 Format(1X,'Time taken ',F0.3,1X,A,' by ',I0,' images, = ',F0.2,' computing power') 26: End If 27: Contains 28: Real(wp) Function f(x) 29: Real(wp),Intent(In) :: x 30: f = 1/x**2 31: End Function 32: End Program
主な変更点を以下に示します。
5行目 - 結果を格納するスカラー変数 res を共配列変数にする
変更前: Real(wp) res,true 変更後: Real(wp) res[*],true
共配列、他の像からも直接アクセスできる変数です。 （配列の各要素を参照する場合と同じように）角括弧と像番号を用いて、どの像の変数であるかを示します。
例１）res[3] = 1.23_wp ! 像3 の res に 1.23 を代入する 例２) Print *, res[1], res[2] ! 像1と像2 の res を出力する
共配列という名前が少し誤解を招いてしまうのですが、共配列はスカラーであっても配列であっても構いません。
（補足：表記方法が配列と似ていた事から「共配列」という言葉が使われるようになりましたが、
意味をとらえる上では「並列」と置き換えた方がイメージがつかみやすいかもしれません。例）並列スカラー変数、並列配列変数）
共配列の宣言は変数名に続けて [*] を記述することで行います。
例） REAL x[*], y(10)[*]ここでxはスカラーの共配列、yは1次元配列の共配列を示します。
今回は各像でそれぞれ得られる部分的な計算結果（スカラー値）を格納する目的で利用するので、 スカラー共配列を用います。
12行目 - 呼び出している求積計算ルーチンを共配列版のものに変更する
変更前: Call rectangle_rule(f,a,b,steps,res) 変更後: Call co_rectangle_rule(f,a,b,steps,res)
区分求積法の計算を行うサブルーチン rectangle_rule を co_rectangle_rule というサブルーチンに置き換えます。 co_rectangle_rule は全体の計算を各像で分担するような役割を持ち、それぞれの像に置いての計算は、 今まで通り quad1.f90 の rectangle_rule を呼び出す事により行います。（詳細後述）
14-26行目 - 結果の出力を像１のみに限定し、像数及び簡略な並列化指標も出力する
共配列プログラムは多くの場合複数の像で実行されます。 そのため結果の出力は一つの像からのみ行うようにします。 この変更を行わないと、Print文の出力がすべての像で行われてしまう事になります。
各像は一意の像番号（整数）を持っていますが、その範囲は１-像の総数と決まっています。 以下では（必ず存在する）像1の場合のみPrint文が実行されるように変更しています。 尚、組込み関数Num_images()は像の総数を返し、This_Image()は現在実行されている自身の像番号を返します。
変更前:
Print *,'Calculated value: ',res
...
Print 1,secs,'seconds'
...
1 Format(1X,'Time taken ',F0.3,1X,A)
変更後:
If (This_Image()==1) Then
Print *,'Calculated value: ',res
...
Print 1,secs,'seconds',Num_Images(),Num_Images()*secs
...
1 Format(1X,'Time taken ',F0.3,1X,A,' by ',I0,' images, = ',F0.2,' computing power')
End If
補足： 上記変更では、並列効率を表す簡便的な指標となる数値（像の総数×経過時間）も併せて出力するようにします。
3.2.2 計算処理を各像に分担するサブルーチン co_rectangle_rule を追加
上記で求積計算ルーチンrectangle_ruleへの呼び出しをco_rectangle_ruleに変更しましたが、 ここでは、そのco_rectangle_ruleを追加します。
co_rectangle_ruleは積分の範囲と像の数、及び自身の像番号を用いて、自身が担当する計算範囲を得て、
その部分のみの求積計算を quad1.f90 の rectangle_rule を呼び出すことにより行います。
尚、以下のコードが示す通り、今回はcoquad1というモジュールにco_rectangle_ruleという総称引用仕様名を定義し、そ
の個別仕様手続きの一つとして、倍精度の引数を取るco_rectangle_dpを定義します。
[
coquad1.f90
] - co_rectangle_ruleを含むモジュール
1: Module coquad1
2: Use quad1
3: Implicit None
4: Private co_rectangle_dp
5: Public co_rectangle_rule
6: Interface co_rectangle_rule
7: Module Procedure co_rectangle_dp
8: End Interface
9: Contains
10: Subroutine co_rectangle_dp(f,dfrom,dto,dnstep,integral)
11: Interface
12: Double Precision Function f(x)
13: Double Precision,Intent(In) :: x
14: End Function
15: End Interface
16: Double Precision,Intent(In) :: dfrom,dto
17: Integer,Intent(In) :: dnstep
18: Double Precision,Intent(Out) :: integral[*]
19: Double Precision from,to
20: Integer i,nstep
21: ! Calculate steps and bounds for this image.
22: nstep = dnstep/Num_Images()
23: from = (dto*(This_Image()-1)+dfrom*(Num_Images()-This_Image()+1))/Num_Images()
24: to = (dto*(This_Image())+dfrom*(Num_Images()-This_Image()))/Num_Images()
25: ! Use the single-image code to calculate the integral on our part of the interval.
26: Call rectangle_rule(f,from,to,nstep,integral)
27: ! We have integrated our interval, now synchronise.
28: Sync All
29: If (This_Image()==1) Then
30: ! Everyone has finished, so sum all the values.
31: Do i=2,Num_Images()
32: integral = integral + integral[i]
33: End Do
34: ! Broadcast the result back to everyone.
35: Do i=2,Num_Images()
36: integral[i] = integral
37: End Do
38: End If
39: Sync All
40: End Subroutine
41: End Module
------------------
5-8行目： co_rectangle_ruleという名前の総称引用仕様名を定義します。
今回のコードでは、一つの個別手続きを定義します。
10-40行目：co_rectangle_ruleの個別手続きco_rectangle_dpを定義します。
co_rectangle_ruleルーチンは逐次版のrectangle_ruleとほぼ同じインターフェースを持ちます。
唯一の違いは計算結果を格納する出力変数 integral が共配列変数である点です。
rectangle_rule: Double Precision,Intent(Out) :: integral
co_rectangle_rule: Double Precision,Intent(Out) :: integral[*]
22-24行目：像の総数と自身の像番号を用いて、自身が担当する計算範囲を得ます。
26行目：自身の計算範囲のみの求積計算をrectangle_ruleを呼び出すことにより行います。
28行目：すべての像を同期します。（すべての像で担当範囲の計算が終わるまで待ちます）
31-33行目：自身が像1の場合は、すべての像の計算結果を足し合せて最終結果を求めます。
35-37行目：像１で求まった最終結果を他の像にも伝搬します。
（今回のサンプルでは最終結果を像１でのみ出力するので、この部分の処理は特に必要ありませんが、
このモジュールに汎用性を持たせるために、このようなコードを含めています）
39行目：すべての伝搬が終わるまで待ちます。
3.2.3 出力結果例
以下に像の数を２、４、８に指定して得られた出力の例を示します。また参考までに、逐次計算版の出力結果も併せて示します。
【逐次計算版の出力結果例（再掲）】 Calculated value: 0.9000000000003117 True value (approx): 0.9000000000000000 Relative error 3.4626622556920158E-13 Time taken 4.950 seconds
【像の数＝２の場合の出力結果例】 Calculated value: 0.9000000000003480 True value (approx): 0.9000000000000000 Relative error 3.8660432879724895E-13 Time taken 2.542 seconds by 2 images, = 5.08 computing power
【像の数＝４の場合の出力結果例】 Calculated value: 0.9000000000002202 True value (approx): 0.9000000000000000 Relative error 2.4461913975907617E-13 Time taken 1.292 seconds by 4 images, = 5.17 computing power
【像の数＝８の場合の出力結果例】 Calculated value: 0.9000000000000347 True value (approx): 0.9000000000000000 Relative error 3.8487731520338762E-14 Time taken .680 seconds by 8 images, = 5.44 computing power