3.5 NECスパコン

3.5.1 NEC SX-Aurora TSUBASA

NEC SX-Aurora TSUBASA [13](以下、SX)を用いた高速化について述べます。
SXは1台にベクトル長256のベクトル演算器を8コア持っています。
内側のループをコンパイラーの自動ベクトル化によってベクトル計算し、 外側のループをOpenMPによって並列計算します。
さらに複数のノードで計算するには領域分割を行いMPIを用いて並列計算します。
本プログラムは前節まででOpenMPとMPIによって並列化されていますので、SXに対応することは比較的容易です。
ただしベクトル化するためにvectorモードを使用する必要があります。

リスト3-5-1にSXのコンパイルメッセージを示します。
48行のループがベクトル化されていることがわかります。 ただし、44行のように強制的にベクトル化する必要があります。


リスト3-5-1 SXのコンパイルメッセージ(vectorモード)
$ cat -n sor.c
     1  #include "ost.h"
     2  #include "ost_prototype.h"
     3
     4  real_t update_rb_vector(int oe)
     5  {
     6          const real_t omega = (real_t)Solver.omega;
     7
     8          const int64_t ni = D(1, 0, 0) - D(0, 0, 0);
     9          const int64_t nj = D(0, 1, 0) - D(0, 0, 0);
    10          const int64_t nk = D(0, 0, 1) - D(0, 0, 0);
    11
    12          real_t res2 = 0;
    13          int64_t i;
    14  #ifdef _OPENMP
    15  #pragma omp parallel for reduction (+: res2)
    16  #endif
    17          for (        i = 0; i <= Nx; i++) {
    18          for (int64_t j = 0; j <= Ny; j++) {
    19                  const int64_t kmin = ((i + j) % 2 == oe) ? 0 : 1;
    20                  const int64_t kmax = Nz;
    21                  const int64_t n = D(i, j, kmin);
    22                  real_t *ptr = &V[n];
    23                  real_t *ptr_vxp = &V[n + ni];
    24                  real_t *ptr_vxm = &V[n - ni];
    25                  real_t *ptr_vyp = &V[n + nj];
    26                  real_t *ptr_vym = &V[n - nj];
    27                  real_t *ptr_vzp = &V[n + nk];
    28                  real_t *ptr_vzm = &V[n - nk];
    29                  real_t *ptr_e   = &Epsr_v[n];
    30                  real_t *ptr_exp = &Epsr_v[n + ni];
    31                  real_t *ptr_exm = &Epsr_v[n - ni];
    32                  real_t *ptr_eyp = &Epsr_v[n + nj];
    33                  real_t *ptr_eym = &Epsr_v[n - nj];
    34                  real_t *ptr_ezp = &Epsr_v[n + nk];
    35                  real_t *ptr_ezm = &Epsr_v[n - nk];
    36                  real_t rxp = RXp[i];
    37                  real_t rxm = RXm[i];
    38                  real_t ryp = RYp[j];
    39                  real_t rym = RYm[j];
    40                  real_t *ptr_rzp = &RZp[kmin];
    41                  real_t *ptr_rzm = &RZm[kmin];
    42                  int    *ptr_id = &idVolt_v[n];
    43  #ifdef __NEC__
    44  #pragma _NEC ivdep
    45  #elif __FUJITSU
    46  #pragma loop novrec
    47  #endif
    48                  for (int64_t k = kmin; k <= kmax; k += 2) {
    49                          const real_t e = *ptr_e;
    50
    51                          const real_t axp = (e + (*ptr_exp)) * rxp;
    52                          const real_t axm = (e + (*ptr_exm)) * rxm;
    53                          const real_t ayp = (e + (*ptr_eyp)) * ryp;
    54                          const real_t aym = (e + (*ptr_eym)) * rym;
    55                          const real_t azp = (e + (*ptr_ezp)) * (*ptr_rzp);
    56                          const real_t azm = (e + (*ptr_ezm)) * (*ptr_rzm);
    57
    58                          const real_t res = (*ptr_id ? 0 : 1) * omega * ((
    59                                  axp * (*ptr_vxp) +
    60                                  axm * (*ptr_vxm) +
    61                                  ayp * (*ptr_vyp) +
    62                                  aym * (*ptr_vym) +
    63                                  azp * (*ptr_vzp) +
    64                                  azm * (*ptr_vzm)) / (axp + axm + ayp + aym + azp + azm) - (*ptr));
    65
    66                          *ptr += res;
    67                          res2 += res * res;
    68
    69                          ptr += 2;
    70                          ptr_vxp += 2;
    71                          ptr_vxm += 2;
    72                          ptr_vyp += 2;
    73                          ptr_vym += 2;
    74                          ptr_vzp += 2;
    75                          ptr_vzm += 2;
    76                          ptr_e += 2;
    77                          ptr_exp += 2;
    78                          ptr_exm += 2;
    79                          ptr_eyp += 2;
    80                          ptr_eym += 2;
    81                          ptr_ezp += 2;
    82                          ptr_ezm += 2;
    83                          ptr_id += 2;
    84                          ptr_rzp += 2;
    85                          ptr_rzm += 2;
    86                  }
    87          }
    88          }
    89
    90          return res2;
    91  }
$ ncc -c -O2 -fopenmp -I../include -fdiag-vector=2 -fdiag-parallel=2 -Wall sor.c
ncc: par(1801): sor.c, line 15: Parallel routine generated.: update_rb_vector$1
ncc: par(1803): sor.c, line 17: Parallelized by "for".
ncc: vec( 101): sor.c, line 18: Vectorized loop.
ncc: vec( 101): sor.c, line 48: Vectorized loop.
ncc: opt(1097): sor.c, line 49: This statement prevents loop optimization.
ncc: vec( 128): sor.c, line 58: Fused multiply-add operation applied.
ncc: vec( 128): sor.c, line 67: Fused multiply-add operation applied.
ncc: vec( 126): sor.c, line 67: Idiom detected.: Sum
ncc: par(1809): sor.c, line 88: Barrier synchronization.


3.5.2 SXの計算時間

表3-5-1にFOCUSスパコンWシステムでVE数を変えたときの計算時間を示します。
8VEは1VEの6.5倍速くなります。
実行コマンドは以下の通りです。

$ ost_ncc -n 8 data/benchmark/benchmark600.ost (1VEのとき)
$ mpiexec -ve 0-7 -n 8 ost_ncc_mpi -n 8 data/benchmark/benchmark600.ost (8VEのとき)

表3-5-1 SXの計算時間(FOCUSスパコンWシステム, benchmark600, vectorモード)
VE数計算時間(速度比)
162.5秒 (1.0)
8 9.6秒 (6.5)