Передача массива C++ в подпрограмму Fortran вызывает появление значений nan в результате

Question

Я пытаюсь использовать код Fortran из приложения C++. В частности, я пытаюсь взаимодействовать с SLATEC drc3jj.f. Однако подпрограмма Fortran возвращает массив, размер которого зависит от параметров, переданных функции.

Если размер массива равен 1, тогда массив C++, который я печатаю, содержит соответствующее значение. Однако, если этот размер больше единицы, массив C++ содержит NaN, где должны быть выходные значения.

Ниже приведен код, который я использую. Это просто связывает подпрограмму Fortran с приложением C++.

#ifndef FORTRANLINKAGE_H 
#define FORTRANLINKAGE_H 

extern "C" 
{ 
    extern void drc3jj_(double*,double*,double*,double*,double*, 
         double*,double [],int*,int*); 
} 

#endif // FORTRANLINKAGE_H 

Мясо здесь нет. На самом деле мы называем Fortran подпрограмму из C++ и напечатать результат:

#include "fortranLinkage.h" 

#include <iostream> 
#include <stdlib.h> 

using namespace std; 

void wigner3j(double l2, double l3, double m2, double m3, double coeff []) 
{ 
    double l1min,l1max; 
    int ierr,size(3); 

    drc3jj_(&l2,&l3,&m2,&m3,&l1min,&l1max,coeff,&size,&ierr); 

    cout << "Min: " << l1min << "\t Max: " << l1max << "\t Err: " << ierr << endl; 
} 

int main(int argc, char const *argv[]) 
{ 
    int l1(atoi(argv[1])),l2(atoi(argv[2])),m2(atoi(argv[3])),m3(atoi(argv[4])); 
    double coeff [3]; 

    wigner3j(l1,l2,m2,m3,coeff); 

    for (int i=0;i<3;i++) 
    { 
      cout << coeff[i] << endl; 
    } 
    return 0; 
} 

Если мы называем программу с ./myProgram 2 8 2 8, он правильно выводит 1/SQRT (21). Однако, если мы попытаемся ./myProgram 2 8 2 7, где размер массива на самом деле 2, мы получаем следующий результат:

Min: 9 Max: 10  Err: 0 
-nan 
-nan 
2.08175e-317 

NaNs на самом деле имеют надлежащий знак.

Во всяком случае, есть ли другой (правильный) способ передать массивы C++ в Fortran? Это даже проблема?

Answer 1

Проблема заключается не в интерфейсе между C++ и Fortran, а в устаревшей реализации Fortran. Файл drc3jj.f является частью библиотеки SLATEC, которая имеет функции утилиты, которые возвращают константы, зависящие от машины, на которой она выполняется (машинные константы). Они определены в файлах d1mach.f, i1mach.f и r1mach.f.

Однако, поскольку Fortran 95, существует, например, встроенные функции huge(), tiny(), spacing() и epsilon(), которые гарантированно вернуть правильные значения для любой машины.

Решение заключается в том, чтобы удалить любые ссылки на d1mach(int) в подпрограмму drc3jj() и заменить их соответствующими характеристиками.

Кроме того, привязка непосредственно к подпрограммам Fortran всегда может быть сложной, поскольку она зависит от компилятора; лучше использовать iso_c_binding в Fortran90 определить интерфейс C для вас в типизированного образом:

!wrapper.f90: 

subroutine drc3jj_wrap(l2, l3, m2, m3, l1min, l1max, thrcof, ndim, ier) bind(C) 

    use iso_c_binding 
    implicit none 

    real(c_double), value, intent(in)   :: l2, l3, m2, m3 
    real(c_double), intent(out)     :: l1min, l1max 
    real(c_double), dimension(ndim), intent(out):: thrcof 
    integer (c_int), value, intent(in)   :: ndim 
    integer (c_int), intent(out)    :: ier 

    interface 
      SUBROUTINE DRC3JJ (L2, L3, M2, M3, L1MIN, L1MAX, THRCOF, NDIM, IER) 
       INTEGER NDIM, IER 
       DOUBLE PRECISION L2, L3, M2, M3, L1MIN, L1MAX, THRCOF(NDIM) 
      end SUBROUTINE DRC3JJ 
    end interface 

    call DRC3JJ(l2, l3, m2, m3, l1min, l1max, thrcof, ndim, ier) 

end subroutine drc3jj_wrap 

и

// fortranLinkage2.h 
#ifndef FORTRANLINKAGE_H 
#define FORTRANLINKAGE_H 

extern "C" 
{ 
    extern void drc3jj_wrap(double l2, double l3, double m2, double m3, 
          double *l1max, double *l2max, double *thrcof, 
          int ndim, int *ier); 
} 

#endif // FORTRANLINKAGE_H 

, а затем вы просто позвоните

void wigner3j(double l2, double l3, double m2, double m3, double coeff []) 
{ 
    double l1min,l1max; 
    int ierr,size(3); 

    drc3jj_wrap(l2,l3,m2,m3,&l1min,&l1max,coeff,size,&ierr); 

    cout << "Min: " << l1min << "\t Max: " << l1max << "\t Err: " << ierr << endl; 
} 

и компиляции и работает

$ g++ -c foo2.cc 
$ gfortran -c wrapper.f90 
$ gfortran -c drc3jj.f 
$ g++ -o foo2 foo2.o wrapper.o drc3jj.o -lgfortran 
$ ./foo2 2 8 2 8 
Min: 10 Max: 10  Err: 0 
0.218218 
2.07738e-317 
0 
gpc-f103n084-$ ./foo2 2 8 2 7 
Min: 9 Max: 10  Err: 0 
-0.102598 
-0.19518 
0