openMP交互方向暗黙のメソッド

Question

こんにちは私はopenMPでこのコードの計算を並列化しようとしています。有限差分法を用いた流体力学渦度の計算である。私はそうするためにAlternating direction暗黙の方法を使用しています。

私はその実行をスピードアップしたいと思います。 （ここではNx = Ny = 100）

問題は、openMpを使用するとコードをスピードアップするのではなく速度を落とすことです。私は共有変数を指定しようとしましたが、これはあまり役に立ちません。

すべてのベスト

void ADI(double vort[][Ny], double psi[][Ny], double n[][Ny], 
double cls[][Ny],double AAx[], double BBx[], double CCx[], double DDx[], 
double AAy[], double BBy[], double CCy[], double DDy[], 
double cx[][Ny], double cy[][Ny], double epsx[][Ny], double epsy[][Ny], 
double vortx[], double vorty[Ny-2], double dx, double Dxs, double coefMass, 
double coefMasCls) 
{ 
    ////////////calcul sur y//////////// 
    //calcul coef ADI 
    int i=0, j=0; 

    #pragma omp parallel for private(Dxs,i) shared(psi,vort) 
    for (i=0; i<Nx; i++) //Boundary condition sur x 
    { 
     vort[i][0]=(psi[i][0]-psi[i][1])*2/Dxs; 
     vort[i][Ny-1] = (psi[i][Ny-1]-psi[i][Ny-2])*2/Dxs; 
    } 

    #pragma omp parallel for private(Dxs,j) shared(psi,vort) 
    for (j=0; j<Ny; j++) //Boundary condition 
    { 
     vort[0][j] = (psi[0][j]-psi[1][j])*2/Dxs; 
     vort[Nx-1][j] = (psi[Nx-1][j]-psi[Nx-2][j])*2/Dxs; 

    } 

    for (j=1; j<Ny-1; j++) //interior points 
    { 
     #pragma omp parallel for private(coefMasCls,coefMasCls,i) shared(psi,vort,n,cls) 
     for (i=1; i<Nx-1; i++) //interior points 
     { 

      vort[i][j] = vort[i][j] - coefMass * (n[i+1][j]-n[i-1][j])- coefMasCls * (cls[i+1][j]-cls[i-1][j]);; 

     } 
     //i=0; 
     //vort[i][j] = vort[i][j] + coefMass*(n[1][j]-n[1][j]); 
     //i=Nx-1; 
     //vort[i][j] = vort[i][j] + coefMass*(n[Nx-2][j]-n[Nx-2][j]); 

    } 


    for (i=1; i<Nx-1; i++) //interior points 
    { 
     for (j=1; j<Ny-1; j++) //interior points 
     { 


      AAy[j] = -.5 * (.5 * (1 + epsy[i][j]) * cy[i][j-1] + dx); 

      BBy[j] = 1 + dx + .5 * epsy[i][j] * cy[i][j]; 

      CCy[j] = .5 * (.5 * (1 - epsy[i][j]) * cy[i][j+1] - dx); 

      DDy[j] = .5 * (.5 * (1 + epsx[i][j]) * cx[i-1][j] + dx) * vort[i-1][j] 
      + (1 - dx - .5 * epsx[i][j] * cx[i][j]) * vort[i][j] 
      + .5 * (- .5 * (1 - epsx[i][j]) * cx[i+1][j] + dx) * vort[i+1][j]; 


      vorty[j] = vort[i][j]; 
     } 


     DDy[1]=DDy[1] - AAy[1] * vort[i][0]; //the AA[0] are not taken into account in the tridiag methode. Include it in the second hand 
     DDy[Ny-2]=DDy[Ny-2] - CCy[Ny-2]* vort[i][Ny-1]; //moving boundary condition 
     //DDy[Ny-3]= DDy[Ny-3]; //vorticity nul on the free slip boundary condition 


     tridiag(AAy, BBy, CCy, DDy, vorty, Ny-1); //ne calcul pas le point en 0 et en Ny-1 

     for (j=1; j<Ny-1; j++) 
     { 
      vort[i][j]=vorty[j]; 
     } 

    } 

    ////////////calcul sur x ////////// 
    //calcul coef ADI 

    for (j=1; j<Ny-1; j++) 
    { 

     for (i=1; i<Nx-1; i++) 

     { 

      AAx[i] = -.5* (.5 * (1 + epsx[i][j]) * cx[i-1][j] + dx); 
      BBx[i] = 1 + dx + .5 * epsx[i][j] * cx[i][j]; 
      CCx[i] = .5 * (.5 * (1 - epsx[i][j]) * cx[i+1][j] - dx) ; 

      DDx[i]= .5 * (.5 * (1 + epsy[i][j]) * cy[i][j-1] + dx) * vort[i][j-1] 
      + (1 - dx - .5 * epsy[i][j] * cy[i][j]) * vort[i][j] 
      + .5 * (-.5 * (1 - epsy[i][j]) * cy[i][j+1] + dx) * vort[i][j+1]; 

      vortx[i]=vort[i][j]; 

     } 


     DDx[1] = DDx[1] - AAx[1]* vort[0][j]; 
     DDx[Nx-2] = DDx[Nx-2] - CCx[Nx-2] * vort[Nx-1][j]; 

     tridiag(AAx, BBx, CCx, DDx, vortx, Nx-1); //ne calcul pas le point en 0 et en Nx-1 

     for (i=1; i<Nx-1; i++) 
     { 
      vort[i][j]=vortx[i]; 

     } 


    } 

} 

Answer 1

まず最初に行うには、並列化が最も悪い影響を与えることができますが、キャッシュスラッシングを経験されるだろうようにそこの最後のループは非常に見えたループ隔離するために実際にあります。 *

double vort[Nx][Ny]; 
// ... 
for (int j=1; j<Ny-1; ++j) { 
    #pragma omp parallel for 
    for (int i=1; i<Nx-1; ++i) { 
     vort[i][j] -= f(i, j); 
    } 
} 

任意の所与のスレッドは、オフセットJ + kにおけるvortの値を読み取り、順番に更新しようとしているNY、j個+（K + 1）* NY、j個+（K + 2：構造体を少し単純化）* Nyなど。forループがスレッド間でどのようにチャンクされるかによって異なります。これらのアクセスのそれぞれは、8バイトを更新するためにキャッシュライン分のデータを引き出す予定です。外側のループが再び始まると、今アクセスしたデータのどれもがキャッシュに残ることはありません。

すべてが同じであるため、アレイアクセスを配列して最小のストライド方向（C配列の場合は最後のインデックス）に移動できる場合は、キャッシュの動作がはるかに優れています。ディメンションサイズ100の場合、配列はそれほど大きくないので、これは大きな違いになります。例えば、 Nx、Ny = 1000の場合、アレイに「間違った方法」でアクセスすることは壊滅的なことでしょう。

これはシリアルコードでパフォーマンスが低下することがありますが、スレッドを追加するだけでそれはずっと悪くなると思います。つまり、これらの内部ループのそれぞれで実行される計算量は非常に小さくなります。あなたはメモリの帯域幅にかかわらず制約を受ける可能性が高いでしょう。

補遺

ただ、明示的であることを、「正しい」ループのアクセスは次のようになります。

for (int i=1; i<Nx-1; ++i) { 
    for (int j=1; j<Ny-1; ++j) { 
     vort[i][j] -= f(i, j); 
    } 
} 

そして、それを並列化するために、あなたがより良いチャンクにスレッド間でデータをコンパイラに許可することができます

#pragma omp parallel for collapse(2) 
for (int i=1; i<Nx-1; ++i) { 
    for (int j=1; j<Ny-1; ++j) { 
     vort[i][j] -= f(i, j); 
    } 
} 

最後に、（スレッドが上の踏み誤った共有を避けるために：collapseディレクティブを使用することにより、お互いのキャッシュライン）、アレイの2つの隣接する行が同じキャッシュライン内のデータを共有していないことを確認するとよいでしょう。各行がメモリ内のキャッシュラインサイズの倍数に揃えられていることを確認したり、各行の最後にパディングを追加するだけで済みます。

double vort[Nx][Ny+8]; // 8 doubles ~ 64 bytes 

（これは十分である、64バイトのキャッシュラインを仮定。）