openMP 대체 방향 암시 적 방법

Question

안녕하세요, 저는이 코드의 계산을 OpenMP로 병렬 처리하려고합니다. 그것은 유한 차분 implicite 방법으로 유체 역학 와도를 계산합니다. 나는 대체 방향 암시 적 방법을 사용하고 있습니다.

나는 그 실행을 빠르게하고 싶다. (여기서 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

제일 먼저 할

루프 parallelizations 가장 나쁜 영향을 미칠 수있는 분리 사실이지만 캐시 스 래싱을 경험하는 것 같은 마지막 루프는 매우이 보인다. 구조를 단순화 비트 :

상관 주어진 스레드 * 옵셋 J + K * NY, J의 + (K + 1)에서 NY, J의 +를 vort의 값을 읽고 차례로 업데이트하려고

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); 
    } 
} 

(K + 2) * Ny 등이 있습니다. for 루프가 쓰레드에서 어떻게 chunk되는지에 따라 달라집니다. 이러한 각각의 액세스는 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); 
    } 
} 

그리고 병렬화하기 위해, 당신은 스레드에서 데이터를 더 청크 컴파일러를 허용 할 수 있습니다 collapse 지시어를 사용하여 :

#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); 
    } 
} 

마지막으로, 거짓 공유를 방지하기 위해 (스레드에 속 서로의 캐시 라인), 배열의 인접한 두 행이 동일한 캐시 라인에서 데이터를 공유하지 않도록하는 것이 좋습니다. 각 행을 메모리의 캐시 라인 크기의 배수로 정렬하거나 각 행의 끝에 패딩을 추가하는 것만으로도 충분합니다.

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