1

C++ 関数から Fortran サブルーチンを呼び出し、その C++ 関数が OpenMP の parallel for 構造内で呼び出されると、Fortran サブルーチンは時々異なる値を返します。これは、同じ入力 (50 個の引数) で同じ結果を返すブラックボックス サブルーチンです。何百もの異なる入力の組み合わせに対して実行するために、サブルーチン呼び出しを並列化しました。プログラムを 2 回実行し、各サブルーチン実行の結果を出力すると、結果は同じではありません。

問題の詳細:

  1. シリアル バージョンは一貫性があり、正常に動作しており、常に同じ答えが得られます。
  2. サブルーチンは疑似乱数を使用しません。
  3. サブルーチンは、同じ .F90 ファイル内の他のサブルーチンを呼び出します。
  4. ネストも、openmp プラグマも、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