c++ - OpenMP 二维拉普拉斯方程代码 : G++ blows up while Intel C converges perfectly

Question

问题:使用有限差分迭代求解拉普拉斯方程 div2(u)=0。

边界条件:u(x,0)=1 u(x,1)=u(0,y)=u(1,y)=0

算法是:

u[i,j]= 0.25 * (u[i,j-1] + u[i,j+1] + u[i-1,j] + u[i+1,j])

环境:Arch Linux + Gcc 4.8.2 & Intel Parallel Studio 2013 SP4。

求解拉普拉斯方程的代码如下:

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <time.h> 
#include <sys/time.h>
using namespace std;
double getTimeElapsed(struct timeval end, struct timeval start)
{
    return (end.tv_sec - start.tv_sec) + (end.tv_usec - start.tv_usec) / 1000000.00;
}
double max(double a,double b)
{
    if (a>b) return a;
    else return b;
}
int main(int argc, char *argv[])
{
    struct timeval t0, t1; 
    double htime;
    double **laplace,prev,qa=0,accu=10e-4;
    int i,j,rowcol,step=0; 
    if (argc != 3) cout << "Usage: hw4 [points at row] [accuracy]" << endl;
    else
    {
        rowcol=atoi(argv[1]);
        accu=atof(argv[2]);
        laplace = new double *[rowcol+1];
        for (i=0;i<=rowcol;i++) 
            laplace[i]=new double [rowcol+1];
        #pragma omp parallel for //MP init
        for (i=0;i<=rowcol;i++)
            for (j=0;j<=rowcol;j++)
            {
                if (j==0) laplace[i][j] = 1.0;
                else laplace[i][j] = 0.0;
            }
        gettimeofday(&t0, NULL);
        while(1)
        {
            #pragma omp parallel for shared(rowcol,laplace,prev) reduction (max:qa) //mp calculation
            for(i=1;i<=rowcol-1;i++)
                for(j=1;j<=rowcol-1;j++)
                {
                    prev=laplace[i][j];
                    laplace[i][j]=0.25*(laplace[i+1][j]+laplace[i-1][j]+laplace[i][j+1]+laplace[i][j-1]);
                    qa=max(qa,abs(laplace[i][j]-prev));
                }
            if (qa < accu)
            {
                gettimeofday(&t1, NULL);
                htime=getTimeElapsed(t1,t0);
                cout << "Done!" << endl << "Time used: " << htime << endl;
                for (i=0;i<=rowcol;i++) delete [] laplace[i];
                delete [] laplace;
                exit(0);
            }
            else 
            {
                step++;
                if (step%80==0) cout << "Iteration = " << step << " ,  Error= " << qa << endl;
                qa=0;
            }
        }
    }
    return 0;
}

我用 100x100 网格和 1e-06 精度测试了程序。

英特尔 C++ 通过 6000 次迭代完美地完成了程序。它产生与顺序代码相同的结果。

但对于 GCC，它未能收敛。

我不知道为什么。

PS:g++编译的程序能运行，但没有像intel版那样收敛。

温斯顿

Answer 1

您在 j 和 prev 中存在竞争条件。将它们设为私有(private)。

#pragma omp parallel for private(j) //MP init
for (i=0;i<=rowcol;i++)
    for (j=0;j<=rowcol;j++)

和

#pragma omp parallel for private(j,prev) reduction (max:qa) //mp calculation
for(i=1;i<=rowcol-1;i++)
    for(j=1;j<=rowcol-1;j++) {
        prev=laplace[i][j];