3

私はCFDコードを使用しています（計算流体力学用）。私は最近、ループの 1 つでインテル® コンパイラーが SSE を使用しているのを見る機会があり、このループの計算パフォーマンスがほぼ 2 倍になりました。ただし、SSE および SIMD 命令の使用は、運がよかったようです。ほとんどの場合、コンパイラは何もしません。

AVX命令が近い将来この側面を強化することを考慮して、私はSSEの使用を強制しようとしています.

簡単な一次元伝熱コードを作りました。これは、もう一方の結果を使用する 2 つのフェーズで構成されます (U0 -> U1、U1 -> U0、U0 -> U1 など)。反復すると、安定した解に収束します。メイン コードのループのほとんどは、同じ種類の計算を使用しています。(有限差分)。

ただし、私のコードは通常のループよりも 2 倍遅くなります。結果は同じなので、計算は一貫しています。

私は間違いを犯しましたか？スーパー コンピューター (Westmer を使用) でテストする前に、Core 2 を使用してループをテストしています。

SSEループと参照ループを含むコードは次のとおりです。

#include <stdio.h>
#include <emmintrin.h>
#include <time.h>
//#include <vector>
#define n1 1004
#define niter 200000

int i,j,t;

double U0[n1] __attribute__ ((aligned(16)));
double U1[n1] __attribute__ ((aligned(16)));
double Dx,Dy,Lx,Ly,InvDxDx,Dt,alpha,totaltime,Stab,DtAlpha,DxDx;

__m128d vmmx00;
__m128d vmmx01;
__m128d vmmx02;
__m128d vmmx10;

__m128d va;
__m128d vb;
__m128d vc;
__m128d vd;

clock_t time0,time1;
FILE *f1;
int main()
{
/* ---- GENERAL ---- */
   alpha = 0.4;
   totaltime = 1.0/100.0;
   Dt = totaltime/((niter-1)*1.0);
   Lx = 1.0;
   Dx = Lx/((n1-1)*1.0);
   InvDxDx = 1.0/(Dx*Dx);
   DxDx = Dx*Dx;
   Stab = alpha*Dt*(InvDxDx);
   DtAlpha = Dt*alpha;
/* Stability if result <= 0.5 */
   printf("Stability factor : %f \n",Stab);


   for( i = 0; i < n1; i++){U0[i] = 0.0;}
   U0[1] = 1.0;
   U0[2] = 1.0;
   U0[3] = 1.0;
   U0[n1-2] = 2.0;

//    for ( i = 0; i < n1; i++) {
 //       for ( j = i + 1; j < n2; j++) {
  //          std::swap(U0[i][j], U0[j][i]);
   //     }
    //}

    va = _mm_set1_pd(-2.0);
    vb = _mm_set1_pd(InvDxDx);
    vd = _mm_set1_pd(DtAlpha);

 time0=clock();

 for( t = 0; t < niter; t++)
 {

    for( i = 2; i < n1-2; i+=2)
    {
      //printf("%d  %d  \n",i,j);
      //fflush(stdout);
      vmmx00 = _mm_load_pd(&U0[i]);
      vmmx01 = _mm_loadu_pd(&U0[i+1]);
      vmmx02 = _mm_loadu_pd(&U0[i-1]);

      vmmx10 = _mm_mul_pd(va,vmmx00);   // U1[i][j] = -2.0*U0[i][j];
      vmmx10 = _mm_add_pd(vmmx10,vmmx01);   // U1[i][j] = U1[i][j] + U0[i+1][j];
      vmmx10 = _mm_add_pd(vmmx10,vmmx02);   // U1[i][j] = U1[i][j] + U0[i-1][j];
      vmmx10 = _mm_mul_pd(vb,vmmx10);   // U1[i][j] = U1[i][j] * InvDxDx;

      vmmx10 = _mm_mul_pd(vd,vmmx10);   // U1[i][j] = U1[i][j] * DtAlpha;
      vmmx10 = _mm_add_pd(vmmx10,vmmx00);   // U1[i][j] = U1[i][j] + U0[i][j];

      _mm_store_pd(&U1[i],vmmx10);

      // U1[i][j] = U0[i][j] + DtAlpha*( (U0[i+1][j]-2.0*U0[i][j]+U0[i-1][j])*InvDxDx
    }

    for( i = 2; i < n1-2; i+=2)
    {
      //printf("%d  %d  \n",i,j);
      //fflush(stdout);
      vmmx00 = _mm_load_pd(&U1[i]);
      vmmx01 = _mm_loadu_pd(&U1[i+1]);
      vmmx02 = _mm_loadu_pd(&U1[i-1]);

      vmmx10 = _mm_mul_pd(va,vmmx00);   // U0[i][j] = -2.0*U1[i][j];
      vmmx10 = _mm_add_pd(vmmx10,vmmx01);   // U0[i][j] = U0[i][j] + U1[i+1][j];
      vmmx10 = _mm_add_pd(vmmx10,vmmx02);   // U0[i][j] = U0[i][j] + U1[i-1][j];
      vmmx10 = _mm_mul_pd(vb,vmmx10);   // U0[i][j] = U0[i][j] * InvDxDx;

      vmmx10 = _mm_mul_pd(vd,vmmx10);   // U0[i][j] = U0[i][j] * DtAlpha;
      vmmx10 = _mm_add_pd(vmmx10,vmmx00);   // U0[i][j] = U0[i][j] + U1[i][j];

      _mm_store_pd(&U0[i],vmmx10);

      // U1[i][j] = U0[i][j] + DtAlpha*( (U0[i+1][j]-2.0*U0[i][j]+U0[i-1][j])*InvDxDx
    }

 }

 time1=clock();
 printf("Loop 0, total time : %f \n", (double) time1-time0);
 f1 = fopen ("out0.dat", "wt");
 for( i = 1; i < n1-1; i++)
 {
       fprintf (f1, "%d\t%f\n", i, U0[i]);
 }


// REF

   for( i = 0; i < n1; i++){U0[i] = 0.0;}
   U0[1] = 1.0;
   U0[2] = 1.0;
   U0[3] = 1.0;
   U0[n1-2] = 2.0;

 time0=clock();

 for( t = 0; t < niter; t++)
 {

    for( i = 2; i < n1-2; i++)
    {
       U1[i] = U0[i] + DtAlpha* (U0[i+1]-2.0*U0[i]+U0[i-1])*InvDxDx;
    }

    for( i = 2; i < n1-2; i++)
    {
       U0[i] = U1[i] + DtAlpha* (U1[i+1]-2.0*U1[i]+U1[i-1])*InvDxDx;
    }

 }

 time1=clock();
 printf("Loop 0, total time : %f \n", (double) time1-time0);
 f1 = fopen ("outref.dat", "wt");
 for( i = 1; i < n1-1; i++)
 {
       fprintf (f1, "%d\t%f\n", i, U0[i]);
 }



}

++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++

編集 ：

++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++

あなたの答えを考慮して、これについて議論するのに適切な場所を見つけたので、トピックを拡大して私の目的を説明します. 受け入れていただければ、すべてのループについて順を追って説明します。これは長くなる可能性がありますが、私のドメイン内の多くの人や、OpenFoam などのオープンソース ソルバーにとって非常に役立つ可能性があります。エネルギー消費への影響を考慮していません (私たちは皆、大規模なスーパーカルキュレーターを使用しています)。

私が使用している CFD コードは、512 の Westmer コアで実行するのに 1 か月以上かかります。MPI (Message Passing Interface) を使用して proc 間の通信を行っています。物理フィールドはメッシュと見なすことができるため、シミュレーションのタイプに応じて、1D、2D、または 3D 配列になります。しかし、3D はあなたが想像できる最高のものです。

完全なコードは Fortran 95 であり、実際には単純化された C です。C とのインターフェイスは簡単で、C ルーチンは Fortran から直接呼び出すことができます。型は同じです (int、double、long など)。しかし、Fortran はそのような最適化を許可していません。単純になるように設計されています。そのため、C 命令を調査しています。

すべての CFD コードで、3 種類のループと MPI メモリ分散という同じ問題に直面しています。最初にループについて説明しましょう:

  • 空間微分 (有限差分と呼ばれる) : ループは、すべてのケースの 1D 畳み込みで構成されます(1D、2D、3D、一度に 1 つの軸でのみ微分します) (DF = F[i-1]*A + F[i ]*B + F[i+1]*C)。ただし、1D 以上を使用する場合、メモリアクセスは次のようになります。

    // x1 derivative
for i 1 -> n1
    for j 1 ->n2
        DF_x1[i][j] = F[i-1][j]*A + F[i][j]*B + F[i+1][j]*C

// x2 derivative
for i 1 -> n1
    for j 1 ->n2
        DF_x2[i][j] = F[i][j-1]*D + F[i][j]*E + F[i][j+1]*G

    最初のループでは、メモリ アクセスが継続されません (Fortran での反転、メモリが反転されます)。これが最初の問題です。3D 配列を使用する場合は同上。

  • ポアソン方程式の解決、つまり行列の乗算: ループは、シミュレーションに応じて、1D、2D、または 3D の畳み込みで構成されます。これは実際には 2 次導関数 (DDF = D(DF)) です。

    for i 1 -> n1
    for j 1 ->n2
        DDF[i][j] = F[i-1][j]*A + F[i][j]*B + F[i+1][j]*C + F[i][j-1]*D + F[i][j]*E + F[i][j+1]*G

    このループは、最初に示したループと同じですが、偶数でも奇数でもなく、直接計算されます。

  • 加重ガウス ザイデル解像度、つまり以下と同じループですが、依存関係があります。

    // even
for i 1 -> n1
    for j 1 ->n2
        F1[i][j] = F0[i-1][j]*A + F0[i][j]*B + F0[i+1][j]*C + F0[i][j-1]*D + F0[i][j]*E + F0[i][j+1]*G

//odd
for i 1 -> n1
    for j 1 ->n2
        F0[i][j] = F1[i-1][j]*A + F1[i][j]*B + F1[i+1][j]*C + F1[i][j-1]*D + F1[i][j]*E + F1[i][j+1]*G

    これは、以前に調査したループです。

次に、別の問題に直面します。メモリの分散です。各コアには独自のメモリがあり、他のコアと共有する必要があります。単純化された最後のループを考えてみましょう:

    for t 1 -> niter

        // even
        for i 1 -> n1-2
                F1[i] = F0[i-1]*A + F0[i]*B + F0[i+1]*C

        //odd
        for i 1 -> n1-2
                F0[i] = F1[i-1]*A + F1[i]*B + F1[i+1]*C

n1=512 としますが、RAM 容量が少ないためローカルメモリに格納できません。メモリは、同じコンピューターではなくネットワーク上にある core0 (1->255) と core1 (256-512) に分散されます。この場合、i=256 の導関数は点 i=255 を認識する必要がありますが、この値は他のプロシージャにあります。他のプロセッサの値を含むメモリは、GHOST メモリと呼ばれます。したがって、ループは次のとおりです。

    ! update boundary memory :
    Share to ghost : core0 : F0[255] -> Network -> F0[0] : core1 (don't forget that for core1, the array restart from 0)
    Share to ghost : core1 : F0[1] -> Network -> F0[256] : core0 (you understand that F0[256] is the ghost for core0, and F0[0] is the ghost for core1)
    // even, each core do this loop.
    for i 1 -> n1-2
            F1[i] = F0[i-1]*A + F0[i]*B + F0[i+1]*C

    ! update boundary memory :
    Share to ghost : core0 : F1[255] -> Network -> F1[0] : core1
    Share to ghost : core1 : F1[1] -> Network -> F1[256] : core0

    //odd, each core do this loop.
    for i 1 -> n1-2
            F0[i] = F1[i-1]*A + F1[i]*B + F1[i+1]*C

これに対処する必要があります。Mysticial、あなたは今私がどこに向かっているのかわかります: ループインターレースはこれを考慮に入れる必要があります. これは次のようにできると思います：

        ! update boundary memory :
        Share to ghost : core0 : F0[255] -> Network -> F0[0] : core1 
        Share to ghost : core1 : F0[1] -> Network -> F0[256] : core0

        for t 1 -> niter

        ! compute borders in advance :
        core0 only : F1[255] = F0[254]*A + F0[255]*B + F0[256]*C
        core1 only : F1[1] = F0[0]*A + F0[1]*B + F0[2]*C

        Launch Share to ghost asynchronous : core0 : F1[255] -> Network -> F1[0] : core1 
        Launch Share to ghost asynchronous : core1 : F1[1] -> Network -> F1[256] : core0

        During the same time (this can be done at the same time because MPI support asynchronous communications)
        // even
        for i 2 -> n1-3 (note the reduced domain)
                F1[i] = F0[i-1]*A + F0[i]*B + F0[i+1]*C

        Check that communications are done.


        ! compute borders in advance :
        core0 only : F0[255] = F1[254]*A + F1[255]*B + F1[256]*C
        core1 only : F0[1] = F1[0]*A + F1[1]*B + F1[2]*C

        Launch Share to ghost asynchronous : core0 : F0[255] -> Network -> F0[0] : core1 
        Launch Share to ghost asynchronous : core1 : F0[1] -> Network -> F0[256] : core0

        //odd, each core do this loop.
        for i 2 -> n1-3
                F0[i] = F1[i-1]*A + F1[i]*B + F1[i+1]*C

        Check that communications are done.

どこかでインデックスを間違えていないことを願っています。ここでは、最も単純な最初のタイプのループを考えてみましょう。その後、同様のループ 2 と 3 を実行できます。目的はこれを行うことです（これは画像処理に似ています）：

    // x1 derivative
    for i 1 -> n1
        for j 1 ->n2
            DF_x1[i][j] = F[i-1][j]*A + F[i][j]*B + F[i+1][j]*C

    // x2 derivative
    for i 1 -> n1
        for j 1 ->n2
            DF_x2[i][j] = F[i][j-1]*D + F[i][j]*E + F[i][j+1]*G

私はそれに取り組んでおり、あなたの推奨事項を考慮して、数時間以内に結果コードを投稿します.

4

1 に答える 1

1

コード内の Mysticial の賢明な言葉:

  xmm7 = InvDxDx * dtAlpha; // precalculate
  xmm6 = -2*xmm7;
  xmm0 = *ptr++;   // has values d0, d1
  xmm1 = *ptr++;   // has values d2, d3
  while (loop--) {
     xmm2  = xmm0;  // take a copy of d0,d1
     xmm0 += xmm1;  // d0+d2, d1+d3
     xmm2 = shufps (xmm2,xmm1, 0x47);  // the middle elements d1,d2 ??
     xmm0 *= xmm7;  // sum of outer elements * factor
     xmm2 *= xmm6;  // -2 * center element * factor
     // here's still a nasty dependency
     xmm2 += xmm0;
     xmm0 = xmm1;    // shift registers
     *out_ptr ++= xmm2;  // flush
     xmm1 = *ptr++;  // read in new values

  }

これは、ループを 2 つのセグメントに分割し、それぞれが 502 エントリ離れて、または異なるタスクから動作し、結果をインターリーブすることで改善できます。これにより、依存関係の連鎖が断ち切られます。

また、シフト レジスタ アプローチ xmm0 <- xmm1, xmm1 <- new は、ループを 2 回展開し、1 つおきに xmm0 と xmm1 の意味を変更することで回避できます。

これは別の大きな問題を指摘しています: 組み込み関数を使用すると、レジスタが常にメモリにスピルされます。

于 2012-11-20T17:47:42.513 に答える