openmp를 사용하여 eigen C++을 사용하는 Jacobi 알고리즘의 병렬화

Question

난 Numerical Recipes 책에 설명 된 루틴을 기반으로 Jacobi 알고리즘을 구현했지만 openmp를 사용하여 병렬화하려고하므로 매우 큰 행렬을 사용하려고합니다.

//#pragma omp parallel for private (ip,g,h,t,theta,c,s,tau) 
//#pragma omp parallel for private (g,h,t,theta,c,s,tau) 

가에서 내 시도입니다

void ROTATE(MatrixXd &a, int i, int j, int k, int l, double s, double tau) 
{ 
double g,h; 
g=a(i,j); 
h=a(k,l); 
a(i,j)=g-s*(h+g*tau); 
a(k,l)=h+s*(g-h*tau); 

} 

void jacobi(int n, MatrixXd &a, MatrixXd &v, VectorXd &d) 
{ 
int j,iq,ip,i; 
double tresh,theta,tau,t,sm,s,h,g,c; 

VectorXd b(n); 
VectorXd z(n); 

v.setIdentity();  
z.setZero(); 

#pragma omp parallel for 
for (ip=0;ip<n;ip++) 
{ 
    d(ip)=a(ip,ip); 
    b(ip)=d(ip); 
} 

for (i=0;i<50;i++) 
{ 
    sm=0.0; 
    for (ip=0;ip<n-1;ip++) 
    { 
     #pragma omp parallel for reduction (+:sm) 
     for (iq=ip+1;iq<n;iq++) 
      sm += fabs(a(ip,iq)); 
    } 
    if (sm == 0.0) { 
     break; 
    } 

    if (i < 3) 
    tresh=0.2*sm/(n*n); 
    else 
    tresh=0.0; 

    #pragma omp parallel for private (ip,g,h,t,theta,c,s,tau) 
    for (ip=0;ip<n-1;ip++) 
    { 
    //#pragma omp parallel for private (g,h,t,theta,c,s,tau) 
     for (iq=ip+1;iq<n;iq++) 
     { 
      g=100.0*fabs(a(ip,iq)); 
      if (i > 3 && (fabs(d(ip))+g) == fabs(d[ip]) && (fabs(d[iq])+g) == fabs(d[iq])) 
      a(ip,iq)=0.0; 
      else if (fabs(a(ip,iq)) > tresh) 
      { 
       h=d(iq)-d(ip); 
       if ((fabs(h)+g) == fabs(h)) 
       { 
        t=(a(ip,iq))/h; 
       } 
       else 
       { 
        theta=0.5*h/(a(ip,iq)); 
        t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); 
        if (theta < 0.0) 
        { 
         t = -t; 
        } 
        c=1.0/sqrt(1+t*t); 
        s=t*c; 
        tau=s/(1.0+c); 
        h=t*a(ip,iq); 

        #pragma omp critical 
        { 
        z(ip)=z(ip)-h; 
        z(iq)=z(iq)+h; 
        d(ip)=d(ip)-h; 
        d(iq)=d(iq)+h; 
        a(ip,iq)=0.0; 


        for (j=0;j<ip;j++) 
         ROTATE(a,j,ip,j,iq,s,tau); 
        for (j=ip+1;j<iq;j++) 
         ROTATE(a,ip,j,j,iq,s,tau); 
        for (j=iq+1;j<n;j++) 
         ROTATE(a,ip,j,iq,j,s,tau); 
        for (j=0;j<n;j++) 
         ROTATE(v,j,ip,j,iq,s,tau); 
        } 

       } 
      } 
     } 
    } 


} 

}

나는 계산의 대부분을 수행하는 루프 코드에 삽입 모두 의견을 병렬화하고 싶었다. 불행히도 둘 다 잘못된 결과를 산출하게됩니다. 나는 문제가이 블록에있을 수 있습니다 의심 :

z(ip)=z(ip)-h; 
z(iq)=z(iq)+h; 
d(ip)=d(ip)-h; 
d(iq)=d(iq)+h; 

축적의 일반적으로 이런 종류는 감소를 필요로하기 때문하지만, 각 스레드가 배열의 다른 부분에 접근하기 때문에,이 특정 아닙니다.

OpenMP로 작업하기 시작한 이래로 올바른 병렬 처리를 수행하고 있는지 확실하지 않습니다. 따라서 제안이나 권장 사항도 환영받을 것입니다.

Sidenote : Eigen의 SelfAdjointEigenSolver를 포함하여 eigenvalue 및 eigenvector 결정에 더 빠른 알고리즘이 있지만 고유 벡터에 필요한 정밀도를 제공하지 못하는 점을 알고 있으며이 알고리즘은 있습니다.

미리 감사드립니다.

편집 : 나는 4096x4096 크기의 시스템에 대한 계산 시간을 줄이지 못하기 때문에 양자 물리학자가 제공 한 답을 정정하는 것으로 생각했습니다. 어쨌든 코드를 수정하여 제대로 작동하도록 만들었습니다. 그리고 아마도 충분히 큰 시스템을 사용하기 위해 사용할 수있을 것입니다.

의

의 #pragma OMP가 실제로 계산 시간을 줄일 수 있을지 테스트하는 타이머의 사용을 권합니다.

Answer 1

도와 드리겠습니다. 그러나 이것이 귀하의 질문에 대한 답변인지 확신 할 수 없습니다.

코드에 많은 문제가 있습니다. 당신을위한 나의 친절한 조언은 : 당신이하는 일의 함의를 이해하지 못한다면 평행을하지 말라.

어떤 이유에서든 모든 것을 병렬로 생각하는 것처럼 보입니다. #pragma for을 사용하면 더 빨리 처리 할 수 ​​있습니다. 이것은 매우 잘못되었습니다. 산란 스레드는 수행하기에 비용이 많이 드는 작업이므로 상대적으로 많은 메모리와 시간이 소요됩니다. 따라서 모든 루프에 대해 #pragma for을 다시 실행하면 모든 루프마다 스레드가 다시 생성되므로 프로그램 속도가 크게 저하됩니다. 단, 행렬이 실제로 크고 계산 시간이 산란 비용보다 >> 큽니다. 그들.

거대한 행렬을 곱하기를 원할 때 요소 현명하게 (그리고 나는 양자 역학에서 기대 값의 합계가 필요했습니다) 비슷한 문제가되었습니다. OpenMP를 사용하기 위해서는 행렬을 선형 배열로 평면화 한 다음 어레이 청크를 각각 스레드에 배포 한 다음 for 루프를 실행해야합니다. 모든 루프 반복에서는 다른 요소와 독립적 인 요소가 사용됩니다. 그들 모두가 독립적으로 진화하게 만들었다. 이것은 꽤 빠르다. 왜? 왜냐하면 나는 스레드를 두 번 다시 생성 할 필요가 없었기 때문입니다.

왜 잘못된 결과가 나옵니까?그 이유는 공유 메모리 규칙을 존중하지 않기 때문이라고 생각합니다. 개의 변수가 동시에 여러 스레드에서 수정되고 있습니다. 그것은 어딘가에 숨어 있고 당신은 그것을 찾아야 만합니다! 예를 들어, z의 기능은 무엇입니까? 그것은 참고로 물건을 소요합니까? 나는 여기에서 볼 무엇 :

z(ip)=z(ip)-h; 
z(iq)=z(iq)+h; 
d(ip)=d(ip)-h; 
d(iq)=d(iq)+h; 

가 안전하지 아주 멀티 스레딩 같은데, 나는 당신이 무슨 일을하는지 이해가 안 돼요. 수정해야 할 참조를 반환하고 있습니까? 이는 스레드가 안전하지 않은 방법입니다. 왜 당신은 깨끗한 배열을 만들고 이것 대신에 그것들을 다루지 않습니까?

디버그 방법 : 작은 예제 (2x2 매트릭스, 아마도)로 시작하고 2 개의 스레드 만 사용하고 무슨 일이 일어나는지 이해하려고하십시오. 디버거를 사용하고 중단 점을 정의하고 스레드간에 공유되는 정보를 확인하십시오.

뮤텍스를 사용하여 공유 할 때 어떤 데이터가 망가 졌는지 확인하는 것도 고려하십시오. Here하는 방법입니다.

내 추천 : 스레드를 한 번만 생성하려는 경우가 아니면 OpenMP를 사용하지 마십시오. C++ 11 때문에 OpenMP가 곧 사라질 것이라고 저는 믿습니다. OpenMP는 C++이 네이티브 멀티 스레딩 구현을하지 않았을 때 아름다웠습니다. 따라서 std::thread을 사용하는 방법을 배우고 그것을 사용하고 스레드에서 많은 것을 실행해야한다면 std::thread으로 스레드 풀을 만드는 법을 배웁니다. This은 멀티 스레딩을 배우기에 좋은 책입니다.