3.6 富士通スパコン

3.6.1 Fujitsu Supercomputer PRIMEHPC

Fujitsu Supercomputer PRIMEHPC FX700 [14](以下、FX)を用いた高速化について述べます。
図3-6-1にFX700の仕様を示します。
FXのCPUは48コアから成り、それぞれ512ビットのSIMDベクトル演算を行うことができます。
内側のループをコンパイラーの自動ベクトル化によってSIMDベクトル計算し、 外側のループをOpenMPによって並列計算します。
さらに複数のノードで計算するには領域分割を行いMPIを用いて並列計算します。
本プログラムは前節まででOpenMPとMPIによって並列化されていますので、FXに対応することは比較的容易です。
ただしベクトル化するためにvectorモードを使用する必要があります。

リスト3-6-1にFXのコンパイルメッセージを示します。
17行以下のループがOpenMPによって並列化され、 48行以下のループがベクトル長16(=512ビット/単精度32ビット)でSIMDベクトル化されていることがわかります。 ただし46行のように強制的にベクトル化する必要があります。
富士通スパコンの注意点は以下の通りです。

  1. マクロ"__FUJITSU"が定義されている。
  2. ベクトル化するループのindexは局所変数にする。
  3. ベクトル化するループのindexは64ビット整数にする。
  4. ベクトル化できないループを強制的にベクトル化するにはループの直前に制御文"#pragma loop novrec"を代入する。 このときコンパイルオプション"-Kocl"が必要。
  5. 一番内側のループではなるべくif文を使わない。


リスト3-6-1 FXのコンパイルメッセージ(vectorモード)
$ fcc -c -Kfast,openmp,ocl,optmsg=2 -I../include -Nlst=t sor.c
$ cat sor.lst
Fujitsu C/C++ Version 4.0.0   Sat Jun  6 17:49:12 2020
Compilation information
  Source file       : sor.c
(line-no.)(optimize)
        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   p           for (        i = 0; i <= Nx; i++) {
       18   p           for (int64_t j = 0; j <= Ny; j++) {
       19   p                   const int64_t kmin = ((i + j) % 2 == oe) ? 0 : 1;
       20   p                   const int64_t kmax = Nz;
       21   p                   const int64_t n = D(i, j, kmin);
       22   p                   real_t *ptr = &V[n];
       23   p                   real_t *ptr_vxp = &V[n + ni];
       24   p                   real_t *ptr_vxm = &V[n - ni];
       25   p                   real_t *ptr_vyp = &V[n + nj];
       26   p                   real_t *ptr_vym = &V[n - nj];
       27   p                   real_t *ptr_vzp = &V[n + nk];
       28   p                   real_t *ptr_vzm = &V[n - nk];
       29   p                   real_t *ptr_e   = &Epsr_v[n];
       30   p                   real_t *ptr_exp = &Epsr_v[n + ni];
       31   p                   real_t *ptr_exm = &Epsr_v[n - ni];
       32   p                   real_t *ptr_eyp = &Epsr_v[n + nj];
       33   p                   real_t *ptr_eym = &Epsr_v[n - nj];
       34   p                   real_t *ptr_ezp = &Epsr_v[n + nk];
       35   p                   real_t *ptr_ezm = &Epsr_v[n - nk];
       36   p                   real_t rxp = RXp[i];
       37   p                   real_t rxm = RXm[i];
       38   p                   real_t ryp = RYp[j];
       39   p                   real_t rym = RYm[j];
       40   p                   real_t *ptr_rzp = &RZp[kmin];
       41   p                   real_t *ptr_rzm = &RZm[kmin];
       42   p                   int    *ptr_id = &idVolt_v[n];
       43             #ifdef __NEC__
       44             #pragma _NEC ivdep
       45             #elif __FUJITSU
       46             #pragma loop novrec
       47             #endif
                       <<< Loop-information Start >>>
                       <<<  [OPTIMIZATION]
                       <<<    SIMD(VL: 16)
                       <<<    PREFETCH(SOFT) : 4
                       <<<     SEQUENTIAL : 4
                       <<<      (unknown): 4
                       <<<    SPILLS :
                       <<<      GENERAL   : SPILL 0  FILL 10
                       <<<      SIMD&FP   : SPILL 0  FILL 0
                       <<<      SCALABLE  : SPILL 0  FILL 0
                       <<<      PREDICATE : SPILL 0  FILL 0
                       <<< Loop-information  End >>>
       48   p      v            for (int64_t k = kmin; k <= kmax; k += 2) {
       49   p      v                    const real_t e = *ptr_e;
       50
       51   p      v                    const real_t axp = (e + (*ptr_exp)) * rxp;
       52   p      v                    const real_t axm = (e + (*ptr_exm)) * rxm;
       53   p      v                    const real_t ayp = (e + (*ptr_eyp)) * ryp;
       54   p      v                    const real_t aym = (e + (*ptr_eym)) * rym;
       55   p      v                    const real_t azp = (e + (*ptr_ezp)) * (*ptr_rzp);
       56   p                           const real_t azm = (e + (*ptr_ezm)) * (*ptr_rzm);
       57
       58   p      v                    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   p      v                    *ptr += res;
       67   p      v                    res2 += res * res;
       68
       69   p      v                    ptr += 2;
       70   p      v                    ptr_vxp += 2;
       71   p      v                    ptr_vxm += 2;
       72   p      v                    ptr_vyp += 2;
       73   p      v                    ptr_vym += 2;
       74   p      v                    ptr_vzp += 2;
       75   p      v                    ptr_vzm += 2;
       76   p      v                    ptr_e += 2;
       77   p      v                    ptr_exp += 2;
       78   p      v                    ptr_exm += 2;
       79   p      v                    ptr_eyp += 2;
       80   p      v                    ptr_eym += 2;
       81   p      v                    ptr_ezp += 2;
       82   p      v                    ptr_ezm += 2;
       83   p      v                    ptr_id += 2;
       84   p      v                    ptr_rzp += 2;
       85   p      v                    ptr_rzm += 2;
       86   p      v            }
       87   p           }
       88   p           }
       89
       90               return res2;
       91             }
Total prefetch num: 4
Optimization messages
  jwd6004s-i  "sor.c", line 48: リダクション演算を含むループ制御変数'k'のループをSIMD化しました。
  jwd8222o-i  "sor.c", line 48: このループで必要なプリフェッチの数が、ハードウェアプリフェッチの許容数を超えたため、prefetch命令を生成しました。
  jwd8664o-i  "sor.c", line 48: ループ内に関数呼出しなどの最適化対象外の命令があるため、ソフトウェアパイプライニングを適用できません。
  jwd8209o-i  "sor.c", line 58: 多項式の演算順序を変更しました。
Statistics information
  Option information
    Command line options : -c -Kfast,openmp,ocl,optmsg=2 -I../include -Nlst=t
    Effective options    : -g0 -mt -Qy -std=gnu11 -x- -O3 -Knoalias_const -Kalign_loops


3.6.2 FXの計算時間

表3-6-1にFOCUSスパコンXシステムで計算条件を変えたときの計算時間を示します。
すべてのコアでOpenMPで並列計算するよりも MPIを用いて4プロセスに分割するほうが2倍以上速いことがわかります。
これは図3-6-2のブロック図と関係すると思われます。

表3-6-1 FXの計算時間(FOCUSスパコンXシステム, benchmark600, vectorモード)
計算条件 計算時間
OpenMP 48スレッド 102.5秒
MPI 4プロセス x OpenMP 12スレッド 44.9秒


図3-6-1 FX700の仕様 [14]


図3-6-2 A64FXのブロック図 (FX1000の場合、 https://pr.fujitsu.com/jp/news/2018/08/22-1.html )