Passing C++ array to Fortran subroutine causes nan values to appear in result

Question

I'm trying to use Fortran code from a C++ application. Specifically, I am trying to interface with SLATEC's drc3jj.f. However, the Fortran subroutine returns an array which size depends on the parameters passed to the function.

If the size of the array is 1, then the C++ array that I print contains the appropriate value. However, if this size is greater than one, the C++ array contains NaNs where there should be output values.

Below is the code I use. This merely links the Fortran subroutine to the C++ application.

#ifndef FORTRANLINKAGE_H
#define FORTRANLINKAGE_H

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

#endif // FORTRANLINKAGE_H

The meat is down here. We actually call the Fortran subroutine from C++ and print the output:

#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;
}

If we call the program with ./myProgram 2 8 2 8, it properly outputs 1/sqrt(21). However, if we try ./myProgram 2 8 2 7, where the size of the array is actually 2, we get this result:

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

The NaNs actually have the proper sign.

Anyway, is there another (proper) way to pass C++ arrays to Fortran? Is this even the issue?

Answer 1

The issue resides not in the interface between C++ and Fortran, but rather in the obsolete Fortran implementation. The file drc3jj.f is part of the SLATEC library, which has utility functions that return constants that depend on the machine it's being run on (machine constants). They are defined in files d1mach.f, i1mach.f and r1mach.f.

However, since Fortran 95, there exists intrinsics such as huge(), tiny(), spacing() and epsilon() that are guaranteed to return the proper values for any machine.

The solution, then, is to remove any reference to d1mach(int) in the drc3jj() subroutine and replace them with the appropriate intrinsics.

Moreover, linking directly to the Fortran subroutines can always be tricky because it's compiler dependent; better is to use iso_c_binding in fortran90 to define the interface to C for you in a typesafe way:

!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

and

// 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

and then you just call

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;
}

and compiling and running gives

$ 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