c++ - 将 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 之间的接口(interface)，而在于过时的 Fortran 实现。文件 drc3jj.f 是 SLATEC 库的一部分，它具有返回常量的实用函数，这些常量取决于运行它的机器(机器常量)。它们在文件 d1mach.f、i1mach.f 和 r1mach.f 中定义。

但是，自 Fortran 95 以来，存在诸如 huge()、tiny()、spacing() 和 等内部函数epsilon() 保证为任何机器返回正确的值。

然后，解决方案是删除 drc3jj() 子例程中对 d1mach(int) 的任何引用，并用适当的内部函数替换它们。

此外，直接链接到 Fortran 子例程总是很棘手，因为它依赖于编译器；更好的方法是在 fortran90 中使用 iso_c_binding 以类型安全的方式为您定义 C 的接口(interface):

!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