将C++数组传递给Fortran子例程会导致nan值出现在结果中

Question

我正在尝试使用C++应用程序中的Fortran代码.具体来说,我正在尝试与SLATEC的drc3jj.f进行交互.但是,Fortran子例程返回一个数组,该数组的大小取决于传递给函数的参数.

如果数组的大小为1,则我打印的C++数组包含适当的值.但是,如果此大小大于1,则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

肉在这里.我们实际上从C++调用Fortran子例程并打印输出:

#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

NaN实际上有正确的标志.

无论如何,是否有另一种(正确的)方法将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