2.4　計算しやすいベクタデータの作成

　SIMD演算では、スカラ配列を入力値として計算する場合がほとんどです。しかし、配列上に格納されたデータが必ずしもSIMD演算できる形式になっているとは限りません。ここでは、SIMD演算に都合のよいベクタデータを作成する方法について解説します。

　画像処理で扱う色情報や3Dグラフィックスで扱う3次元座標などのデータは、必ずしもSIMD演算しやすいデータの並びになっているとは限りません。例えば、画像を構成する各画素のRGB値からその画素の輝度Yを算出する場合を考えてみます。RGB値から輝度Yを求める式は、以下のとおりです。

　この計算式に合わせてSIMD演算により輝度を算出するためには、図 2.7に示すように「R要素」、「G要素」、「B要素」ごとにまとめられたベクタデータが必要になります。

図 2.7 SIMD演算による輝度の計算

　しかし、画像ファイル等から読み込まれた色情報は、図2.7 (a)に示すように「RGB」の単位で保持されている場合があり、必ずしもSIMD演算できるデータの並びになっていません。このような場合に、図 2.8 (b)のようにベクタデータを並べ替える必要があります。

　それでは、実際にベクタデータを並べ替える方法について解説します。ベクタデータの並べ替えは、vec_perm()関数という組み込み関数を用いて実行します。vec_perm()関数は、並べ替えの対象となる2つのベクタデータと並べ替えのパターンを引数として、新しいベクタデータを生成して関数の値として返します。

用例 (2-7) vec_perm()関数の記述方法

vc = vec_perm(va, vb, vpat);

　第1引数vaと第2引数vbは並べ替え対象となるベクタデータを指定します。第3引数vpatには、並べ替えパターンを指定します。並べ替えパターンの各バイトは、並べ替え対象となるva、vbを連結して32バイトの配列とみなした場合の配列の要素番号を指定します。図 2.9にvec_perm()関数の処理イメージを示します。

　それでは、並べ替えの例題として転置行列を求めるプログラムを考えてみます。

　転置行列とは、m行n列の行列Aに対してAの要素A_ijと要素A_jiを入れ替えたn行m列の行列で、A^Tと表記します。主に3Dグラフィックスの座標変換などに使用されます。

図 2.10 転置行列

　例題プログラムは、図 2.10のような4×4行列の転置行列を求めるために、vec_perm()関数を用いて行列の各要素の並べ替えを行います。

リスト (2-3) 転置行列を求めるプログラム

 1 int a[16]  __attribute__((aligned(16))) = {  1,  2,  3,  4,
 2                                              5,  6,  7,  8,
 3                                              9, 10, 11, 12,
 4                                             13, 14, 15, 16 };
 5 int aT[16] __attribute__((aligned(16)));
 6
 7 vector signed int *va  = (vector signed int *) a;
 8 vector signed int *vaT = (vector signed int *) aT;
 9 vector signed int vtmp[4];
10 vector unsigned char vpat1 = (vector unsigned char) { 0x00, 0x01, 0x02, 0x03,
11                                                       0x04, 0x05, 0x06, 0x07,
12                                                       0x10, 0x11, 0x12, 0x13,
13                                                       0x14, 0x15, 0x16, 0x17 };
14 vector unsigned char vpat2 = (vector unsigned char) { 0x08, 0x09, 0x0a, 0x0b,
15                                                       0x0c, 0x0d, 0x0e, 0x0f,
16                                                       0x18, 0x19, 0x1a, 0x1b,
17                                                       0x1c, 0x1d, 0x1e, 0x1f };
18 vector unsigned char vpat3 = (vector unsigned char) { 0x00, 0x01, 0x02, 0x03,
19                                                       0x10, 0x11, 0x12, 0x13,
20                                                       0x08, 0x09, 0x0a, 0x0b,
21                                                       0x18, 0x19, 0x1a, 0x1b };
22 vector unsigned char vpat4 = (vector unsigned char) { 0x04, 0x05, 0x06, 0x07,
23                                                       0x14, 0x15, 0x16, 0x17,
24                                                       0x0c, 0x0d, 0x0e, 0x0f,
25                                                       0x1c, 0x1d, 0x1e, 0x1f };
26
27 vtmp[0] = vec_perm(va[0], va[2], vpat1);
28 vtmp[1] = vec_perm(va[1], va[3], vpat1);
29 vtmp[2] = vec_perm(va[0], va[2], vpat2);
30 vtmp[3] = vec_perm(va[1], va[3], vpat2);
31
32 vaT[0] = vec_perm(vtmp[0], vtmp[1], vpat3);
33 vaT[1] = vec_perm(vtmp[0], vtmp[1], vpat4);
34 vaT[2] = vec_perm(vtmp[2], vtmp[3], vpat3);
35 vaT[3] = vec_perm(vtmp[2], vtmp[3], vpat4);

　この例題では、行列の1要素が4バイト長であるため、4つのベクタ変数を使用します。vec_perm()関数では、1回の実行で2つのベクタ変数の並べ替えしかおこなえないため、図 2.11に示すように2段階の並べ替えを行い、合計で8回の並べ替え操作が必要になります。

図 2.11 並べ替え操作

　それでは、例題プログラムのソースコード全体について解説します。

例題プログラム (2-2) 転置行列プログラム全体

 1 #include <stdio.h>
 2 #include <altivec.h>
 3
 4 int a[16]  __attribute__((aligned(16))) = {  1,  2,  3,  4,
 5                                              5,  6,  7,  8,
 6                                              9, 10, 11, 12,
 7                                             13, 14, 15, 16 };
 8 int aT[16] __attribute__((aligned(16)));
 9
10 void print_matrix(int *matrix)
11 {
12     int i, j;
13
14     for (i = 0; i < 4; i++) {
15         for (j = 0; j < 4; j++) {
16             printf("%2d ", matrix[i * 4 + j]);
17         }
18         printf("\n");
19     }
20     return;
21 }
22
23 int main(int argc, char **argv)
24 {
25     vector signed int *va  = (vector signed int *) a;
26     vector signed int *vaT = (vector signed int *) aT;
27     vector signed int vtmp[4];
28     vector unsigned char vpat1 = (vector unsigned char) { 0x00, 0x01, 0x02, 0x03,
29                                                           0x04, 0x05, 0x06, 0x07,
30                                                           0x10, 0x11, 0x12, 0x13,
31                                                           0x14, 0x15, 0x16, 0x17 };
32     vector unsigned char vpat2 = (vector unsigned char) { 0x08, 0x09, 0x0a, 0x0b,
33                                                           0x0c, 0x0d, 0x0e, 0x0f,
34                                                           0x18, 0x19, 0x1a, 0x1b,
35                                                           0x1c, 0x1d, 0x1e, 0x1f };
36     vector unsigned char vpat3 = (vector unsigned char) { 0x00, 0x01, 0x02, 0x03,
37                                                           0x10, 0x11, 0x12, 0x13,
38                                                           0x08, 0x09, 0x0a, 0x0b,
39                                                           0x18, 0x19, 0x1a, 0x1b };
40     vector unsigned char vpat4 = (vector unsigned char) { 0x04, 0x05, 0x06, 0x07,
41                                                           0x14, 0x15, 0x16, 0x17,
42                                                           0x0c, 0x0d, 0x0e, 0x0f,
43                                                           0x1c, 0x1d, 0x1e, 0x1f };
44
45     printf("--- original matrix ---\n");
46     print_matrix(a);
47
48     /* vec_perm() part 1 */
49     vtmp[0] = vec_perm(va[0], va[2], vpat1);
50     vtmp[1] = vec_perm(va[1], va[3], vpat1);
51     vtmp[2] = vec_perm(va[0], va[2], vpat2);
52     vtmp[3] = vec_perm(va[1], va[3], vpat2);
53
54     printf("--- transform 1 ---\n");
55     print_matrix((int *) vtmp);
56
57     /* vec_perm() part 2 */
58     vaT[0] = vec_perm(vtmp[0], vtmp[1], vpat3);
59     vaT[1] = vec_perm(vtmp[0], vtmp[1], vpat4);
60     vaT[2] = vec_perm(vtmp[2], vtmp[3], vpat3);
61     vaT[3] = vec_perm(vtmp[2], vtmp[3], vpat4);
62
63     printf("--- transform 2 ---\n");
64     print_matrix(aT);
65
66     return 0;
67 }