HPC
C/C++
Fortran
MIC

メモリーレイアウト変換

同カテゴリーの次の記事

インテル® コンパイラーのプラグマ/ディレクティブ

この記事は、インテル® ソフトウェア・サイトに掲載されている「Memory Layout Transformations」の日本語参考訳です。

はじめに

この記事では、よく使用される構造体配列 (AoS) から配列構造体 (SoA) へのデータ構造の変更について説明します。この変更により、プロセッサー上のデータに効率良くアクセスできるようなコードをコンパイラーが生成できるようになります。

トピック

効率を改善するユーザーコードの変更の 1 つに、メモリーのギャザー (読み取り) /スキャッター (書き込み) 操作の排除があります。このような不規則なメモリー操作は、レイテンシーと帯域幅の両方を増加させるだけでなく、コンパイラーによるベクトル化を制限します。一部のアプリケーションでは、構造体配列 (AOS) から配列構造体 (SOA) へデータ構造を変更することで、複数の配列要素にわたる構造体の 1 つのフィールドにアクセスするギャザー操作を回避し、配列を反復するループがコンパイラーによってベクトル化されるようにできます。

このようなデータ構造の変更は、そのデータ構造が使用されるすべての位置を考慮して、プログラムレベルで (プログラマーが) 行う必要があります。この変更をループレベルで行った場合、ループの前後でデータ構造の変換を行わなければならないため、コストがかかります。

AOS から SOA への変換: ベクトル化されるコードにおいてギャザー/スキャッター操作を回避する一般的な最適化は、構造体配列 (AOS) から配列構造体 (SOA) へデータ構造を変更することです。構造体フィールドごとに個別の配列を保持することで、連続したメモリーアクセスになり、構造体要素がベクトル化されます。AOS 構造はギャザーとスキャッターを必要とするため、SIMD の効率性が低下するだけでなく、メモリーアクセスにかかわる帯域幅とレイテンシーが増加します。ハードウェアでギャザー/スキャッター・メカニズムが利用できる場合であっても、このデータ構造の変更は必要です。一般にギャザー/スキャッター操作は、連続したロードと比べて非常に多くの帯域幅とレイテンシーが要求されます。

SOA 形式は、オリジナルの構造体要素の複数のフィールドへのアクセスに際して局所性が低下するという欠点があります。そのため、例えば、TLB (トランスレーション・ルックアサイド・バッファー) の負担が増加します。オリジナルのソースコードのデータアクセス・パターン (ループ内で構造体の複数のフィールドにアクセスしているかどうか、あるいはすべての構造体要素にアクセスするか、サブセットにアクセスするか) によっては、代わりに配列構造体配列 (AOSOA) に変換したほうが良いケースもあります。

AOSOA へ変換することで、外側では局所性の利点を活用し、最も内側ではユニットストライドにできます。ユニットストライドでベクトル化をフル活用するため、内側の配列はベクトル長の小さな倍数にすることができます。各構造体要素の複数のフィールドは、このレイアウト (小さな配列) になります。外側は、(より大きな) 構造体配列にすることができます。このレイアウトでは、フィールドアクセスがまだ十分に近いため、オリジナルのコードの隣接する構造体要素へのアクセスにおいてページの局所性が得られます。

以下は、各形式の違いを示すサンプルコードです。

AOS 形式 (AOS):

struct node {
  float x, y, z;
};
struct node NODES[1024];

float dist[1024];
for(i=0; i<1024; i+=16){
   float x[16],y[16],z[16],d[16];
   x[:] = NODES[i:16].x;
   y[:] = NODES[i:16].y;
   z[:] = NODES[i:16].z;
   d[:] = sqrtf(x[:]*x[:] + y[:]*y[:] + z[:]*z[:]);
   dist[i:16] = d[:];
}

SOA 形式 (SOA):

struct node1 {
  float x[1024], y[1024], z[1024];
}
struct node1 NODES1;

float dist[1024];
for(i=0; i<1024; i+=16){
  float x[16],y[16],z[16],d[16];
  x[:] = NODES1.x[i:16];
  y[:] = NODES1.y[i:16];
  z[:] = NODES1.z[i:16];
  d[:] = sqrtf(x[:]*x[:] + y[:]*y[:] + z[:]*z[:]);
  dist[i:16] = d[:];
}

AOSOA 形式 (AOSOA またはタイルされた AOS):

struct node2 {
  float x[16], y[16], z[16];
}
struct nodes2 NODES2[64];

float dist[1024];
for(i=0; i<64; i++){
   float x[16],y[16],z[16],d[16];
   x[:] = NODES2[i].x[:];
   y[:] = NODES2[i].y[:];
   z[:] = NODES2[i].z[:];
   d[:] = sqrtf(x[:]*x[:] + y[:]*y[:] + z[:]*z[:]);
   dist[i*16:16] = d[:];
}

配列表記による AOS データ構造から SOA データ構造へのソース変更例

以下は、配列表記を使用して AOS 形式 (オリジナルのコード) から SOA 形式 (変更後のコード) へ変換する、LIBOR* ベンチマークの OpenMP* のサンプルコードです。

オリジナルのサンプルコード (AOS 形式を使用する LIBOR):

long tid = (long)(tid_);
long long int chunk = npath/nthreads;
long long int beg = tid*chunk;
long long int end = min(npath, (tid+1)*chunk);
long long int chunk_inner = PATH_BLOCK_SIZE/nthreads // JSP
printf("chunk = %d, chunk_inner = %dn", chunk, chunk_inner);

 #pragma omp parallel
 for(long long int path = beg; path < end; path += chunk_inner){ 

      #pragma omp for nowait
      for (long long int path_=0; path_ < chunk_inner; path_+=1) {
             int    i, j, k;
             __declspec(align(64)) FPPREC B[nmat], S[nmat], L[n];
             FPPREC sqez, lam, con1, v_scal, vrat;
             FPPREC b, s, swapval;
             FPPREC *zlocal=z+nmat*(tid*chunk_inner+path_);

             for(i=0;i < n;i++) {
                     L[i] = L0[i];
             }
             for(j=0; j < nmat; j++){
                     sqez = SQRT(delta)*zlocal[j];
                     v_scal = ZERO;

                     for (i=j+1; i < n; i++) {
                             lam  = lambda[i-j-1];
                             con1 = delta*lam;
                             v_scal   += con1*L[i]/(ONE+delta*L[i]);
                             vrat = EXP(con1*v_scal + lam*(sqez-HALF*con1));
                             L[i] = L[i]*vrat;
                     }
             }
             b = ONE;
             s = ZERO;
             for (j=nmat; j < n; j++) {
                     b = b/(ONE+delta*L[j]);
                     s = s + delta*b;
                     B[j-nmat] = b;
                     S[j-nmat] = s;
             }

             v_scal = ZERO;
             for (i=0; i < nopt; i++){
                     k = maturities[i] -1;
                     swapval = B[k] + swaprates[i]*S[k] - ONE;
                     if (swapval < ZERO)
                             v_scal += MHUND*swapval;
             }

             // 割引率を適用 //
             for (j=0; j < nmat; j++){
                     v_scal = v_scal/(ONE+delta*L[j]);
             }
             v[tid*chunk_inner+path_]+=v_scal;
      }
}

変更後のサンプルコード (SOA 形式を使用する LIBOR):

int tid = (int)(tid_);
long long int chunk = (npath%nthreads)?(1 + npath/nthreads):(npath/nthreads);
long long int beg = tid*chunk;
long long int end = min(npath, (tid+1)*chunk);
long long int chunk_inner = ((PATH_BLOCK_SIZE%(nthreads*SIMDVLEN))?(1 + PATH_BLOCK_SIZE/(nthreads*SIMDVLEN)):(PATH_BLOCK_SIZE/(nthreads*SIMDVLEN)))*SIMDVLEN; 

 #pragma omp parallel
 for(long long int path = beg; path < end; path += chunk_inner){  

      #pragma omp for nowait
      for (long long int path_=0; path_ < chunk_inner; path_+=SIMDVLEN) {
             int    i, j, k;
              __declspec(align(64)) FPPREC B[nmat][SIMDVLEN], S[nmat][SIMDVLEN], L[n][SIMDVLEN];
             __declspec(align(64)) FPPREC sqez[SIMDVLEN], lam[SIMDVLEN], con1[SIMDVLEN], v_scal[SIMDVLEN], vrat[SIMDVLEN], vrat_tmp[SIMDVLEN];
             __declspec(align(64)) FPPREC b[SIMDVLEN], s[SIMDVLEN], swapval[SIMDVLEN];
              FPPREC *zlocal=z+nmat*(tid*chunk_inner + path_);
             __declspec(align(64)) FPPREC old_output[SIMDVLEN]; 

             old_output[:] = v[tid*chunk_inner+path_:SIMDVLEN];  

             for(i=0;i < n;i++) {
                     L[i][:] = L0[i];
             } 

             for(j=0; j < nmat; j++) {
                     sqez[:] = SQRT(delta)*zlocal[0:SIMDVLEN];
                     v_scal[:] = ZERO;
                     zlocal += SIMDVLEN;
                     for (i=j+1; i < n; i++) {
                             lam[:]  = lambda[i-j-1];
                             con1[:] = delta*lam[:];
                             v_scal[:]   += con1[:]*L[i][:]/(ONE+delta*L[i][:]);
                             vrat_tmp[:] = (con1[:]*v_scal[:] + lam[:]*(sqez[:]-HALF*con1[:]));
                             vrat[:] = EXP(vrat_tmp[:]);
                             L[i][:] = L[i][:]*vrat[:];
                     }
             } 

             b[:] = ONE;
             s[:] = ZERO; 

             for (j=nmat; j < n; j++) {
                     b[:] = b[:]/(ONE+delta*L[j][:]);
                     s[:] = s[:] + delta*b[:];
                     B[j-nmat][:] = b[:];
                     S[j-nmat][:] = s[:];
             } 

             v_scal[:] = ZERO; 

             for (i=0; i < nopt; i++){
                     k = maturities[i] -1;
                     swapval[:] = B[k][:] + swaprates[i]*S[k][:] - ONE;
                     if (swapval[:] < ZERO)
                             v_scal[:] += MHUND*swapval[:];
             } 

             // 割引率を適用 //
             for (j=0; j < nmat; j++){
                     v_scal[:] = v_scal[:]/(ONE+delta*L[j][:]);
             }
             v[tid*chunk_inner+path_:SIMDVLEN] = old_output[:] + v_scal[:];
      }
}

まとめ

データ構造は、論理的には構造体配列 (AoS) のほうが適しているように見えますが、実際にはメモリーの読み取り (ギャザー) と書き込み (スキャッター) が非常に困難です。また、コンパイラーによる多くの最適化を阻み、特に効率良いベクトル化の妨げとなります。このデータ構造を使用するアルゴリズムは、多くの場合、インテル® Xeon プロセッサーとインテル® メニー・インテグレーテッド・コア (インテル® MIC) アーキテクチャーにおいてパフォーマンスの低下の原因となります。

可能な限り、配列構造体 (SoA) に変更することが望ましいでしょう。実際のアプリケーションでは、これは大変な労力を要するかもしれませんが、そうすることで非常に大きな利点が得られる可能性があります。

次のステップ

この記事は、「Programming and Compiling for Intel® Many Integrated Core Architecture」(英語) の一部「Memory Layout Transformations」の翻訳です。インテル® Xeon Phi™ コプロセッサー上にアプリケーションを移植し、チューニングを行うには、各リンクのトピックを参照してください。アプリケーションのパフォーマンスを最大限に引き出すために必要なステップを紹介しています。

インテル® MIC アーキテクチャーの準備に戻る

コンパイラーの最適化に関する詳細は、最適化に関する注意事項を参照してください。

関連記事