3.7 NECスパコン

3.7.1 NEC SX-Aurora TSUBASA

NEC SX-Aurora TSUBASA [23](以下、SX)について述べます。
SXは1VE(ベクトルエンジン)にベクトル長256のベクトル演算器を8コア持っています。
各VEは32GB(HBM2)のメモリーを持っています。

3.7.2 SXにおける高速化

SXの性能を生かすには以下の高速化技術が必要です。

  1. ベクトル化:一番内側のループをコンパイラーの自動ベクトル化機能によってベクトル化します。 必要なら強制ベクトル化を行います。
  2. スレッド並列化(共有メモリー):OpenMPまたはコンパイラーの自動並列化機能によってスレッド並列を行います。 通常、粒度の一番大きい一番外側のループを並列化します。
  3. MPI並列化(分散メモリー):問題を分割し、MPIを使用して並列化します。 複数のVEを使用するときはMPI並列化が必須です。
本プログラムは前節まででOpenMPとMPIによって並列化されているので、 SXに対応するにはベクトル化が必要になります。

3.7.3 SXの計算時間

表3-7-1にWシステム(A300-8VE)の1VEと8VEの計算時間を示します。
表より8VEでは1VEの最大6.9倍速くなることがわかります。
またmatrixモードはnomatrixモードの2〜3倍速くなります。
8VEでは使用可能なメモリー容量が8倍になり、より大きな問題が計算できるようになります。
セル数が小さいときは計算の速いmatrixモードを使用し、 メモリーが足りなくなったら計算時間は増えますがnomatrixモードを使用してください。
実行コマンドは以下の通りです。

$ ./oth_ncc -n 8 -matrix data/benchmark/benchmark100.oth (1VEのとき)
$ mpiexec -ve 0-7 -n 8 ./oth_ncc_mpi -n 8 -matrix data/benchmark/benchmark100.oth (8VEのとき)

表3-7-1 SXの計算時間(Wシステム、()内は1VEとの速度比)
VE数benchmark100benchmark200benchmark300benchmark400
nomatrixmatrixnomatrixmatrixnomatrixmatrixnomatrixmatrix
120.0秒 (1.0)3.9秒 (1.0)89.0秒 (1.0)26.6秒 (1.0)329.0秒 (1.0)92.9秒 (1.0)625.5秒 (1.0)メモリー不足
8 3.9秒 (5.1)1.5秒 (2.6)16.0秒 (5.6) 6.5秒 (4.1) 47.6秒 (6.9)17.1秒 (5.4) 93.1秒 (6.7)36.8秒

3.7.4 コンパイルメッセージ

(1) nomatrixモード

3.1からnomatrixモードの主要部はmatrixE{x|y|z}.cとなります。
リスト3-7-1にコンパイルメッセージを示します。 一番外側のループ(25行目)がOpenMPによって並列化され、 一番内側のループ(34行目)は複雑ですが、 28行目のように強制ベクトル化を行うとベクトル化されることがわかります。

リスト3-7-1 nomatrixモードの主要部(matrixEx.c)のコンパイルメッセージ

$ cat -n matrixEx.c | head -n 50
     1  /*
     2  matrixEx.c
     3  u = Av (Ex)
     4  */
     5
     6  #include "oth.h"
     7  #include "complex.h"
     8  #include "complex13.h"
     9
    10  void matrixEx(int pre, int ifreq, d_complex_t a[], const d_complex_t v[], d_complex_t u[], d_complex_t b[])
    11  {
    12          const int m0 = ifreq * NMaterial;
    13
    14          int64_t i0 = i0Ex;
    15          int64_t i1 = i1Ex;
    16          int64_t j0 = j0Ex;
    17          int64_t j1 = j1Ex;
    18          int64_t k0 = k0Ex;
    19          int64_t k1 = k1Ex;
    20
    21          int64_t i;
    22  #ifdef _OPENMP
    23  #pragma omp parallel for
    24  #endif
    25          for (        i = i0; i <= i1; i++) {
    26          for (int64_t j = j0; j <= j1; j++) {
    27  #if defined(__NEC__)
    28  #pragma _NEC ivdep
    29  #elif defined(__FUJITSU)
    30  #pragma loop novrec
    31  #elif defined(__CLANG_FUJITSU)
    32  #pragma clang loop vectorize(assume_safety)
    33  #endif
    34          for (int64_t k = k0; k <= k1; k++) {
    35                  d_complex_t c[NNZ];
    36                  d_set13(c, d_complex(0, 0));
    37
    38                  if (MEX(i, j, k) == PEC) {
    39                          c[0] = d_complex(1, 0);
    40                  }
    41                  else {
    42                          if ((Abc == MUR) && (
    43                                  (!Pbc.ypbc && ((j == 0) || (j == Ny))) ||
    44                                  (!Pbc.zpbc && ((k == 0) || (k == Nz))))) {
    45                                  // Mur
    46                                  c[0] = d_complex(1, 0);
    47                                  d_complex_t cv = Material[MEX(i, j, k) + m0].cv;
    48                                  d_complex_t cel = d_complex(0, 0);
    49                                  if      (!Pbc.ypbc && (j ==  0)) {
    50                                          cel = d_rmul(KDYn[ 0], cv);
(以下略)
$ ncc -c -O2 -I../include -fopenmp -fdiag-vector=2 -fdiag-parallel=2 matrixEx.c
ncc: par(1801): matrixEx.c, line 23: Parallel routine generated.: matrixEx$1
ncc: par(1803): matrixEx.c, line 25: Parallelized by "for".
ncc: opt(3008): matrixEx.c, line 25: Reference within a conditional branch moved outside loop.
ncc: opt(3008): matrixEx.c, line 26: Reference within a conditional branch moved outside loop.
ncc: opt(3008): matrixEx.c, line 34: Reference within a conditional branch moved outside loop.
ncc: vec( 101): matrixEx.c, line 34: Vectorized loop.
ncc: vec( 128): matrixEx.c, line 83: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 84: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 85: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 86: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 87: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 88: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 92: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 122: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 127: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 131: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 134: Fused multiply-add operation applied.
ncc: vec( 128): matrixEx.c, line 152: Fused multiply-add operation applied.

(2) matrixモード

3.1からmatrixモードの主要部はprodmvE{x|y|z}.cとなります。
リスト3-7-2にコンパイルメッセージを示します。 一番外側のループ(21行目)がOpenMPによって並列化され、 一番内側のループ(43行目)は強制ベクトル化を行うことなくベクトル化されることがわかります。

リスト3-7-2 matrixモードの主要部(prodmvEx.c)のコンパイルメッセージ

$ cat -n prodmvEx.c
     1  /*
     2  prodmvEx.c
     3  matrix mode
     4  */
     5
     6  #include "oth.h"
     7
     8  void prodmvEx(d_complex_t a[], d_complex_t v[], d_complex_t u[])
     9  {
    10          const int64_t i0 = i0Ex;
    11          const int64_t i1 = i1Ex;
    12          const int64_t j0 = j0Ex;
    13          const int64_t j1 = j1Ex;
    14          const int64_t k0 = k0Ex;
    15          const int64_t k1 = k1Ex;
    16
    17          int64_t i;
    18  #ifdef _OPENMP
    19  #pragma omp parallel for
    20  #endif
    21          for (        i = i0; i <= i1; i++) {
    22          for (int64_t j = j0; j <= j1; j++) {
    23                  int64_t n = NEX(i, j, k0);  // NEXK=1
    24                  d_complex_t *aptr = &a[n * NNZ];
    25                  d_complex_t *v00p = &v[NEX(i,     j,     k0    )];
    26                  d_complex_t *v01p = &v[NEX(i,     j + 1, k0    )];
    27                  d_complex_t *v02p = &v[NEX(i,     j - 1, k0    )];
    28                  d_complex_t *v03p = &v[NEX(i,     j,     k0 + 1)];
    29                  d_complex_t *v04p = &v[NEX(i,     j,     k0 - 1)];
    30                  d_complex_t *v05p = &v[NEY(i + 1, j,     k0    )];
    31                  d_complex_t *v06p = &v[NEY(i,     j,     k0    )];
    32                  d_complex_t *v07p = &v[NEY(i + 1, j - 1, k0    )];
    33                  d_complex_t *v08p = &v[NEY(i,     j - 1, k0    )];
    34                  d_complex_t *v09p = &v[NEZ(i + 1, j,     k0    )];
    35                  d_complex_t *v10p = &v[NEZ(i,     j,     k0    )];
    36                  d_complex_t *v11p = &v[NEZ(i + 1, j,     k0 - 1)];
    37                  d_complex_t *v12p = &v[NEZ(i,     j,     k0 - 1)];
    38  #if defined(__FUJITSU)
    39  #pragma loop novrec
    40  #elif defined(__CLANG_FUJITSU)
    41  #pragma clang loop vectorize(assume_safety)
    42  #endif
    43          for (int64_t k = k0; k <= k1; k++, n++) {
    44                  const d_complex_t a00 = *aptr++;
    45                  const d_complex_t a01 = *aptr++;
    46                  const d_complex_t a02 = *aptr++;
    47                  const d_complex_t a03 = *aptr++;
    48                  const d_complex_t a04 = *aptr++;
    49                  const d_complex_t a05 = *aptr++;
    50                  const d_complex_t a06 = *aptr++;
    51                  const d_complex_t a07 = *aptr++;
    52                  const d_complex_t a08 = *aptr++;
    53                  const d_complex_t a09 = *aptr++;
    54                  const d_complex_t a10 = *aptr++;
    55                  const d_complex_t a11 = *aptr++;
    56                  const d_complex_t a12 = *aptr++;
    57
    58                  const d_complex_t v00 = *v00p++;
    59                  const d_complex_t v01 = *v01p++;
    60                  const d_complex_t v02 = *v02p++;
    61                  const d_complex_t v03 = *v03p++;
    62                  const d_complex_t v04 = *v04p++;
    63                  const d_complex_t v05 = *v05p++;
    64                  const d_complex_t v06 = *v06p++;
    65                  const d_complex_t v07 = *v07p++;
    66                  const d_complex_t v08 = *v08p++;
    67                  const d_complex_t v09 = *v09p++;
    68                  const d_complex_t v10 = *v10p++;
    69                  const d_complex_t v11 = *v11p++;
    70                  const d_complex_t v12 = *v12p++;
    71
    72                  u[n].r =
    73                          (a00.r * v00.r) - (a00.i * v00.i) +
    74                          (a01.r * v01.r) - (a01.i * v01.i) +
    75                          (a02.r * v02.r) - (a02.i * v02.i) +
    76                          (a03.r * v03.r) - (a03.i * v03.i) +
    77                          (a04.r * v04.r) - (a04.i * v04.i) +
    78                          (a05.r * v05.r) - (a05.i * v05.i) +
    79                          (a06.r * v06.r) - (a06.i * v06.i) +
    80                          (a07.r * v07.r) - (a07.i * v07.i) +
    81                          (a08.r * v08.r) - (a08.i * v08.i) +
    82                          (a09.r * v09.r) - (a09.i * v09.i) +
    83                          (a10.r * v10.r) - (a10.i * v10.i) +
    84                          (a11.r * v11.r) - (a11.i * v11.i) +
    85                          (a12.r * v12.r) - (a12.i * v12.i);
    86
    87                  u[n].i =
    88                          (a00.r * v00.i) + (a00.i * v00.r) +
    89                          (a01.r * v01.i) + (a01.i * v01.r) +
    90                          (a02.r * v02.i) + (a02.i * v02.r) +
    91                          (a03.r * v03.i) + (a03.i * v03.r) +
    92                          (a04.r * v04.i) + (a04.i * v04.r) +
    93                          (a05.r * v05.i) + (a05.i * v05.r) +
    94                          (a06.r * v06.i) + (a06.i * v06.r) +
    95                          (a07.r * v07.i) + (a07.i * v07.r) +
    96                          (a08.r * v08.i) + (a08.i * v08.r) +
    97                          (a09.r * v09.i) + (a09.i * v09.r) +
    98                          (a10.r * v10.i) + (a10.i * v10.r) +
    99                          (a11.r * v11.i) + (a11.i * v11.r) +
   100                          (a12.r * v12.i) + (a12.i * v12.r);
   101          }
   102          }
   103          }
   104  }
$ ncc -c -O2 -I../include -fopenmp -fdiag-vector=2 -fdiag-parallel=2 prodmvEx.c
ncc: par(1801): prodmvEx.c, line 19: Parallel routine generated.: prodmvEx$1
ncc: par(1803): prodmvEx.c, line 21: Parallelized by "for".
ncc: vec( 101): prodmvEx.c, line 43: Vectorized loop.
ncc: vec( 128): prodmvEx.c, line 72: Fused multiply-add operation applied.
ncc: vec( 128): prodmvEx.c, line 87: Fused multiply-add operation applied.