台形公式では,関数を閾値[x_min,x_max]をn個の領域に分割し,各領域を台形に近似することで積分値を算出する. 台形公式による積分算出プログラムを以下に示す.まずヘッダ. そして本体. これを呼び出すドライバを以下に示す. 実行結果.

x = [0.000000,100.000000] n = 120
S = 4999.999512

上の台形公式プログラムを並列プログラム化することを考える.この問題では,積分領域をP*N個に分割する場合,P個のプロセッサに領域の1/Pの積分値を計算させて合計値を算出することで並列化できる.各プロセッサでは領域をさらに適当な個数に分割する.以下に台形公式による積分算出プログラムの並列プログラム版ドライバを示す. 以下実行結果.

> mpirun -np 2 a.out 
P0 :: x = [0.000000,50.000000] : n = 60 : s = 1250.000000
P1 :: x = [50.000000,100.000000] : n = 60 : s = 3750.001221
TOTAL :: x = [0.000000,100.000000] : n = 120 : S = 5000.000977
> mpirun -np 3 a.out 
P0 :: x = [0.000000,33.333332] : n = 40 : s = 555.555542
P1 :: x = [33.333332,66.666664] : n = 40 : s = 1666.665771
P2 :: x = [66.666664,100.000000] : n = 40 : s = 2777.779541
TOTAL :: x = [0.000000,100.000000] : n = 120 : S = 5000.000977
> mpirun -np 4 a.out 
P0 :: x = [0.000000,25.000000] : n = 30 : s = 312.500000
P1 :: x = [25.000000,50.000000] : n = 30 : s = 937.499878
P2 :: x = [50.000000,75.000000] : n = 30 : s = 1562.499756
P3 :: x = [75.000000,100.000000] : n = 30 : s = 2187.500977
TOTAL :: x = [0.000000,100.000000] : n = 120 : S = 5000.000488

各プロセスが計算した閾値,分割数,積分値および全体の閾値,分割数,積分値を表示される.

台形公式による面積算出には,閾値が分割数で割り切れないと誤差が生じる.そもそも図形を台形に近似する方法なのだから,閾値が分割数で割り切れたとしても誤差は生じる.後者の問題は分割数を大きくすれば誤差を小さくできるが,その場合は,アンダーフローの問題を回避しなくてはならない.並列計算では,これら数値計算におけるよくある問題とは別の問題がさらに追加される.負荷分散の問題である.上のdaikei-para.cで

#define DIV_NUM 120		/* 閾値の分割数  */

を

#define DIV_NUM 100		/* 閾値の分割数  */

として実行してみる.

> mpirun -np 2 a.out 
P0 :: x = [0.000000,50.000000] : n = 50 : s = 1250.000000
P1 :: x = [50.000000,100.000000] : n = 50 : s = 3750.000000
TOTAL :: x = [0.000000,100.000000] : n = 100 : S = 5000.000000
> mpirun -np 3 a.out 
P0 :: x = [0.000000,33.000000] : n = 33 : s = 544.500000
P1 :: x = [33.000000,66.000000] : n = 33 : s = 1633.500000
P2 :: x = [66.000000,99.000000] : n = 33 : s = 2722.500000
TOTAL :: x = [0.000000,100.000000] : n = 100 : S = 4900.500000
> mpirun -np 4 a.out 
P0 :: x = [0.000000,25.000000] : n = 25 : s = 312.500000
P1 :: x = [25.000000,50.000000] : n = 25 : s = 937.500000
P2 :: x = [50.000000,75.000000] : n = 25 : s = 1562.500000
P3 :: x = [75.000000,100.000000] : n = 25 : s = 2187.500000
TOTAL :: x = [0.000000,100.000000] : n = 100 : S = 5000.000000

プロセス数3の時に誤差が大きくなっている.よく見ると,計算閾値が[0,100]のはずが[0,99]になっている.これは分割数がプロセス数で割り切れないことによる誤差である.ここで

#define DIV_NUM 1000		/* 閾値の分割数  */

としてみる.

> mpirun -np 3 a.out 
P0 :: x = [0.000000,33.299999] : n = 333 : s = 554.445984
P1 :: x = [33.299999,66.599998] : n = 333 : s = 1663.326660
P2 :: x = [66.599998,99.900002] : n = 333 : s = 2772.217285
TOTAL :: x = [0.000000,100.000000] : n = 1000 : S = 4989.990234

分割数を大きくすると計算されない閾値が小さくなる.しかし冒頭の誤差がある上に計算されない閾値による誤差があるというのはどうにも気持ちが悪い.では,余った[99,100]は誰が計算すべきだろうか.この半端な領域をあるプロセス,例えばプロセスランクが最大のプロセスに実行させようとすると,daikei-para.cは以下のようになる(assign関数を変更). これを

#define DIV_NUM 1000		/* 閾値の分割数  */

で実行すると,

> mpirun -np 3 a.out 
P0 :: x = [0.000000,33.299999] : n = 333 : s = 554.445984
P1 :: x = [33.299999,66.599998] : n = 333 : s = 1663.326660
P2 :: x = [66.599998,100.000000] : n = 334 : s = 2782.211914
TOTAL :: x = [0.000000,100.000000] : n = 1000 : S = 4999.984375

計算されない閾値はなくなり,先程より誤差が小さくなった.今回は,規模が小さいので,assign関数はこの程度のロジックでそれなりに効果をあげることができた.ただし,プロセス数が大きくなると,プログラムによっては,一部のプロセスにかかる負荷が大きくなり,並列化効率が下がることがあるので,問題に応じて注意する必要がある.