別解 (tune-1a)

ディレクティブを用いずにほぼ同じ i-loop に対するコンパイラ SIMD 化を実現するには、次のようにループを分割しループの順序を入れ替えるという書き方もあります。コードの変更量は大きくなりますが他に応用の効くかもしれない手法でもありここで紹介しておきます。

axhelm-1a.cpp
 79        for(int k = 0; k < p_Nq; k++) for(int j = 0; j < p_Nq; ++j)
 80        {
 81            for(int i = 0; i < p_Nq; ++i) {
 82                qr[i] = 0.f;
 83                qs[i] = 0.f;
 84                qt[i] = 0.f;
 85            }
 86            for(int m = 0; m < p_Nq; m++) {
 87                for(int i = 0; i < p_Nq; ++i) {
 88                    qr[i] += s_D[i][m] * s_q[k][j][m];
 89                    qs[i] += s_D[j][m] * s_q[k][m][i];
 90                    qt[i] += s_D[k][m] * s_q[m][j][i];
 91                }
 92            }
 93
 94            for(int i = 0; i < p_Nq; ++i) {
 95                // const dlong id = i + j * p_Nq + k * p_Nq * p_Nq + element * p_Np;
 96                const dlong gbase = element * p_Nggeo * p_Np + k * p_Nq * p_Nq + j * p_Nq + i;
 97                const dfloat r_G00 = ggeo[gbase + p_G00ID * p_Np];
 98                const dfloat r_G01 = ggeo[gbase + p_G01ID * p_Np];
 99                const dfloat r_G11 = ggeo[gbase + p_G11ID * p_Np];
100                const dfloat r_G12 = ggeo[gbase + p_G12ID * p_Np];
101                const dfloat r_G02 = ggeo[gbase + p_G02ID * p_Np];
102                const dfloat r_G22 = ggeo[gbase + p_G22ID * p_Np];
103
104                s_Gqr[k][j][i] = r_G00 * qr[i] + r_G01 * qs[i] + r_G02 * qt[i];
105                s_Gqs[k][j][i] = r_G01 * qr[i] + r_G11 * qs[i] + r_G12 * qt[i];
106                s_Gqt[k][j][i] = r_G02 * qr[i] + r_G12 * qs[i] + r_G22 * qt[i];
107            }
108        }
109
110        for(int k = 0; k < p_Nq; ++k) for(int j = 0; j < p_Nq; ++j)
111        {
112            for(int i = 0; i < p_Nq; ++i) {
113                qr[i] = 0.f;
114                qs[i] = 0.f;
115                qt[i] = 0.f;
116            }
117            for(int m = 0; m < p_Nq; m++) {
118                for(int i = 0; i < p_Nq; ++i) {
119                    qr[i] += s_D[m][i] * s_Gqr[k][j][m];
120                    qs[i] += s_D[m][j] * s_Gqs[k][m][i];
121                    qt[i] += s_D[m][k] * s_Gqt[m][j][i];
122                }
123            }
124            for(int i = 0; i < p_Nq; ++i) {
125                const dlong id = element * p_Np + k * p_Nq * p_Nq + j * p_Nq + i;
126                Aq[id] = qr[i] + qs[i] + qt[i];
127            }
128        }
129    }
130}

i-loop が 3 分割されその真中のブロックが m-loop の内側に入っていること、スカラ変数だった qr, qs, qt が長さ8の配列に変更されていることに注意してください。

このバージョンでは上述の方法と完全に同一の結果とまではならなかったのですが 321 GFLOP/s が得られています。