Openmp c++：错误：折叠的循环不完全嵌套答案

【问题标题】：Openmp c++: error: collapsed loops not perfectly nestedOpenmp c++：错误：折叠的循环不完全嵌套
【发布时间】：2021-08-19 06:06:39
【问题描述】：

我有以下想要并行的串行代码。我理解在嵌套循环中使用 collapse 子句时，在 for(i) 循环之前和之后不要有代码是很重要的，因为这是不允许的。那么我如何将嵌套的 for 循环与这样的 if 语句并行：

void foo2D(double Lx, double Ly, double KX[], double KY[], double ksqu[], double ninvksqu[], int FourierMeshType){
    
// Make a variable to account for way wavenumbers are set in FFTW3. 
    
int k_counter = 0;      
    
// kx corresponds to modes [0:n/2-1 , -n/2:-1]. This is why there is an additional step, just due to FFTW3's structuring
    #pragma omp parallel for collapse(2) 
    for(int i = 0; i < nx ; i++){
        for(int j = 0; j < nyk; j++){           
            if (i < nx/2){ // i<nx/2 --> 0 to 127
                KX[j + nyk*i] =2.*M_PI*i/Lx; //K = 2pi/L* [0:nx/2-1 , -n/2:-1]' : creates a column vector with nx/2 elements starting from 0 ( from 0 to 127 then -128 to -1) 256/2 is 128      
            }
            if( i >= nx/2){ // i >= nx/2 --> from -128 to -1
                KX[j + nyk*i] =  2.* M_PI * (-i + 2.*k_counter) / Lx;           
            }
        }
        if( i >= nx/2){ // error here 
            k_counter++;
        }
    }
}

其中 nx, ny, nyk定义为：

static const int nx = 256;
static const int ny = 256;
static const int nyk = ny/2 + 1;

有没有更好的方法来重写这个？

【问题讨论】：

标签： c++ multithreading parallel-processing openmp nested-loops

【解决方案1】：

正如1201ProgramAlarm 在 cmets 中指出的那样，您可以通过消除两个循环之间存在的if 分支来消除错误：

#pragma omp parallel for collapse(2) 
for(int i = 0; i < nx ; i++){
    for(int j = 0; j < nyk; j++){    
       ...       
    }
    if( i >= nx/2){ <== remove this
        k_counter++;
    }
}

这个：

if( i >= nx/2){ <== remove this
    k_counter++;
}

可以预先静态计算。所以您的代码如下所示：

void makeFourierMesh2D(double Lx, double Ly, double KX[], double KY[], double ksqu[], double ninvksqu[], int FourierMeshType){
    
    #pragma omp parallel for collapse(2) 
    for(int i = 0; i < nx ; i++){
        for(int j = 0; j < nyk; j++){           
            if (i < nx/2){
                KX[j + nyk*i] =2.*M_PI*i/Lx; 
            }
            if( i >= nx/2){
                int k_counter = i - nx/2;
                KX[j + nyk*i] =  2.* M_PI * (-i + 2.*k_counter) / Lx;           
            }
        }
    }
}

然后您可以将(-i + 2.*k_counter) 简化为(-i + 2 * ( i - nx/2))，然后将(-i + 2*i - nx) 简化为i - nx。

要进一步删除其他循环，您可以在单独的一对循环中计算每个分支，如下所示：

    #pragma omp parallel for collapse(2) 
    for(int i = 0; i < nx/2 ; i++){
        for(int j = 0; j < nyk; j++){           
             KX[j + nyk*i] =2.*M_PI*i/Lx; 
        }
    }

    #pragma omp parallel for collapse(2) 
    for(int i = nx/2; i < nx ; i++){
        for(int j = 0; j < nyk; j++){           
             KX[j + nyk*i] =  2.* M_PI * ((i - nx) / Lx);    
        }
    }

和以前一样，您也可以将(-i + 2.*k_counter) 简化为i - nx。

【讨论】：