3.8　演習問題 (3-3) 区分求積プログラム

【問題】

　出題プログラム (3-3) は、図 3.8の斜線で示された領域の面積を近似的に計算するプログラムです。このプログラムを複数のSPEを利用して計算するプログラムに書き換えなさい。

　求める領域の面積は、グラフの区間で表すことができ、関数は以下の数式で与えられています。

　このプログラムは、図 3.9に示すように求める面積の区間をn等分に分割して矩形面積()の総和を求める区分求積法を利用しています。func()関数は関数の値を計算します。calc_integral()関数は、細かく分割された矩形面積の総和を計算します。

　なお、関数は、以下のように書き直すことで、vec_madd()関数を用いて簡単にSIMD演算に変換することができます。

【出題プログラム (3-3)】

 1 #include <stdio.h>
 2 #include <altivec.h>
 3
 4 vector float func(vector float vx)
 5 {
 6     vector float vy;
 7     vector float va = (vector float) {   5.0f,   5.0f,   5.0f,   5.0f };
 8     vector float vb = (vector float) { -16.0f, -16.0f, -16.0f, -16.0f };
 9     vector float vc = (vector float) { -36.0f, -36.0f, -36.0f, -36.0f };
10     vector float vd = (vector float) {  64.0f,  64.0f,  64.0f,  64.0f };
11     vector float ve = (vector float) { 192.0f, 192.0f, 192.0f, 192.0f };
12
13     vy = vec_madd(va, vx, vb);
14     vy = vec_madd(vy, vx, vc);
15     vy = vec_madd(vy, vx, vd);
16     vy = vec_madd(vy, vx, ve);
17
18     return vy;
19 }
20
21 float calc_integral(float start, float end, float delta)
22 {
23     int i;
24     float *sum;
25
26     vector float vx     = (vector float) { start+delta*0, start+delta*1,
27                                            start+delta*2, start+delta*3 };
28     vector float vsum   = (vector float) { 0.0f, 0.0f, 0.0f, 0.0f };
29     vector float vdelta = (vector float) { delta, delta, delta, delta };
30     vector float vstep  = (vector float) { 4.0f, 4.0f, 4.0f, 4.0f };
31
32     for (i = 0; i <= (end-start)/delta; i += 4) {
33         vsum = vec_madd(func(vx), vdelta, vsum);
34         vx   = vec_madd(vdelta, vstep, vx);
35     }
36
37     sum = (float *) &vsum;
38
39     return (sum[0] + sum[1] + sum[2] + sum[3]);
40 }
41
42 int main(int argc, char **argv)
43 {
44     float start = 0.0f;
45     float end   = 4.0f;
46     float delta = 0.0001f;
47     float result;
48
49     printf("start = %f, end = %f\n", start, end);
50
51     result = calc_integral(start, end, delta);
52
53     printf("result = %f\n", result);
54
55     return 0;
56 }

【解答方針】

　この問題は、複数のSPEにアプリケーションの処理を分割し、割り振る方法を理解しているかがポイントになります。複数のSPEを利用する方法については第3.5節を参考にしてください。

　なお、PPEプログラムのmain()関数は書き換えず、calc_integral()関数からSPEプログラムを実行し計算結果を返すように書き換えてください。