многопоточность — вызов подпрограммы Fortran из C ++ для OpenMP должен работать параллельно?

При вызове подпрограммы Fortran из функции C ++, а функция C ++ вызывается внутри OpenMP-параллели для конструкции, подпрограмма Fortran время от времени возвращает разные значения. Это подпрограмма черного ящика, которая должна возвращать тот же результат с тем же вводом (50 аргументов). Я распараллелил вызов подпрограммы, чтобы запустить его для сотен различных комбинаций ввода. Если я запускаю программу дважды и печатаю результаты каждого выполнения подпрограммы, результаты не совпадают.

Подробности о проблеме:

1. Серийная версия последовательна и работает нормально, постоянно дает один и тот же ответ;
2. Подпрограмма не использует псевдослучайные числа;
3. Подпрограмма вызывает другие подпрограммы в том же файле .F90;
4. Нет вложенности и нет открытых прагм или нет внутри фортрановых подпрограмм;
5. Если я пытаюсь использовать функции OpenMP API внутри подпрограмм Фортрана, они возвращают бессмысленную информацию;
6. Я использую -fautomatic, -fopenmp и -frecursive при компиляции с gfortran (кстати, я использую gcc 5.2.0), и все подпрограммы были сделаны RECURSIVE. Все компилируется и связывается нормально, и проблема действительно появляется, когда я запускаю .exe.
7. Подпрограммы Fortran не имеют доступа к вводу / выводу. Все переменные передаются через аргументы. Там нет общих или сохраненных блоков. Все подпрограммы используют фиктивные аргументы, а выходные переменные явно инициализируются внутри каждой подпрограммы;
8. Я не использую какие-либо предложения OpenMP с параллельным для #pragma omp.
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