1

当从 C++ 函数调用 Fortran 子例程并且在 OpenMP 并行 for 构造中调用 C++ 函数时，Fortran 子例程会不时返回不同的值。它是一个黑盒子例程，应该返回相同的结果和相同的输入（50 个参数）。我并行化了子例程调用，以便为数百种不同的输入组合运行它。如果我运行程序两次并打印每个子程序执行的结果，结果是不一样的。

有关问题的详细信息：

  1. 串行版本一致且工作正常，始终给出相同的答案；
  2. 该子程序不使用伪随机生成数；
  3. 子程序调用同一个.F90文件中的其他子程序；
  4. 没有嵌套，没有 openmp pragma 或包含在 fortran 子例程中；
  5. 如果我尝试在 Fortran 子例程中使用 OpenMP API 函数，它们会返回乱码；
  6. 我在使用 gfortran 编译时使用 -fautomatic、-fopenmp 和 -frecursive（顺便说一句，我使用的是 gcc 5.2.0）并且所有子例程都设置为 RECURSIVE。一切都很好地编译和链接，当我运行.exe时问题确实出现了。
  7. Fortran 子例程不访问 I/O。所有变量都通过参数传递。没有 COMMON 或 SAVED 块。所有子程序都使用虚拟参数，并且输出变量在每个子程序内显式初始化；
  8. 我没有使用任何带有#pragma omp parallel for 的OpenMP 子句。
  9. 如果线程数小于可用处理器数，则结果之间的差异数会减少。将线程与处理器绑定并不能解决问题。

代码很大，但我设法在一个示例中对其进行了简化以说明问题：

//StdAfx.h
#include "other.h"
#include <omp.h>
//...many other includes 

//first.cpp
#include "StdAfx.h"
typedef struct
{
    float x[51];
    float result;
} A;
typedef A *B;
B old=NULL;
int size;
float weight;
int var;

int main()
{
    size = 100;
    old = new (nothrow) A[size];
    long* control=NULL;
    control = new long[size];
    int kk;
    //...
    //changing some control[] values to -1
    var = 5; weight = 0.7;
    //...
    #pragma omp parallel for
    for (kk=0; kk<=size-1; kk++)
    {
        if (control[kk]>-1) old[kk].result = calcresult(old[kk].x,kk);
    }
    ...
    delete [] old;
    oldpop = NULL;
    delete [] control;
    control = NULL;
}
float calcresult(float *x, int k)
{
    int dev=0;
    double kresult;
    dev = 10;
    kresult = othercalcresult(&x[0],k);
    kresult += (weight*dev*double(1.0/var));
    return(kresult);
}

//other.h
float othercalcresult(float *x, int anyk=0);

//other.cpp
extern "C" {
void _stdcall fdlf_(int VET[93],int *N, double *extresult);
}
double anothercalcresult(float *x, int *iVet)
{
    int iN=1;
    double extresult=0.0;
    //stuff here
    //original fortran subroutine has 50 arguments
    fdlf_(iVet,&iN,&extresult);
    return(extresult);
}
float othercalcresult(float *x, int anyk=0)
{
    unsigned int i,ii;
    float otherresult=0.0;
    int ulimit;
    //ulimit = 10;
    //iVet is a two dimensional array iVet
    int** iVet = new int*[numcenarios_anaprog_local];
    for (ii=0; ii<ulimit; ii++) iVet[ii]=new int[93];
    //initialize new vector
    for (i=0; i<ulimit; i++) 
        for (ii=0; ii<93; ii++) 
            iVet[i][ii]=(100*i)+ii;
    double* partialresult=NULL;
    partialresult= new double[ulimit];
    for (int jj=0;jj<ulimit;jj++) partialresult[jj] = 0.0;
    //stuff here
    for (i=0;i<ulimit;i++) partialresult[i] = anothercalcresult(x,iVet[i])
    for (i=0;i<ulimit;i++) otherresult+=float(partialresult[i]);
    return(otherresult);
}

//EXT.F90
RECURSIVE SUBROUTINE AUXSUB1(N,VALUE1)
    INTEGER N
    REAL*8 VALUE1
    VALUE1 = 1 / (2 ** N)
    RETURN
END SUBROUTINE AUXSUB1

RECURSIVE SUBROUTINE AUXSUB2(N,VALUE2)
    INTEGER N
    REAL*8 VALUE2
    VALUE2 = 1 / (3 ** N)
    RETURN
END SUBROUTINE AUXSUB2

RECURSIVE SUBROUTINE FDLF(VET,N,EXTRESULT)
    INTEGER VET(93),N
    REAL*8 VALUE1, VALUE2, EXTRESULT
    VALUE1 = 0.
    VALUE2 = 0.
    EXTRESULT = 0.0
    CALL AUXSUB1(N,VALUE1)
    CALL AUXSUB2(N,VALUE2)
    DO I=1,93
        IF I.LT.47 THEN 
            EXTRESULT = EXTRESULT + VALUE1
        ELSE
            EXTRESULT = EXTRESULT + VALUE2
        END IF
    END DO
    EXTRESULT = 1 / EXTRESULT
    RETURN
END SUBROUTINE FDLF
4

0 回答 0