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