how to call fortran routines from C++?

Question

I wish to call fortran routine cbesj.f from my C++ code and how do I achieve this?

Here are steps I have done:

1. Download cbesj.f plus dependencies from netlib amos webpage, http://www.netlib.org/cgi-bin/netlibfiles.pl?filename=/amos/cbesj.f

2. In the source dir,

   f2c -C++PR *.f

   g++ -c *.c

   ar cr libmydemo.a *.o

3. [test_cbesj.cpp][1] and [mydemo.h][2] are used to call the subroutine in this way,

   g++ test_cbesj.cpp -lf2c -lm -L. -lmydemo it returns bug:

test_cbesj.cpp:(.text+0xd6): undefined reference to `cbesj_(complex*, float*, long*, long*, complex*, long*, long*)'

What shall be the proper way to refer to cbesj_ subroutine in my problem? Thanks!

Thanks for casey: I think your approach is the best. But I still have set fault, why? Here we go:

f77 -c *.f    

in modemo.h

//File mydemo.h                                                                                                                              
#ifndef MYDEMO_H
#define MYDEMO_H
#include <stdio.h>      /* Standard Library of Input and Output */
#include "f2c.h"

extern"C" int cacai_(complex *z__, real *fnu, integer *kode, integer *mr,         integer *n, complex *y, integer *nz, real *rl, real *tol, real *el\
im, real *alim);
extern"C" int cairy_(complex *z__, integer *id, integer *kode, complex *ai,         integer *nz, integer *ierr);
extern"C" int casyi_(complex *z__, real *fnu, integer *kode, integer *n,     complex *y, integer *nz, real *rl, real *tol, real *elim, real     \
  *alim);
extern"C" int cbesj_(complex *z__, real *fnu, integer *kode, integer *n,     complex *cy, integer *nz, integer *ierr);
extern"C" int cbinu_(complex *z__, real *fnu, integer *kode, integer *n,     complex *cy, integer *nz, real *rl, real *fnul, real *tol, real *el\
im, real *alim);
extern"C" int cbknu_(complex *z__, real *fnu, integer *kode, integer *n, complex *y, integer *nz, real *tol, real *elim, real *alim);
extern"C" int cbuni_(complex *z__, real *fnu, integer *kode, integer *n,     complex *y, integer *nz, integer *nui, integer *nlast, real *fnul, \
real *tol, real *elim, real *alim);
extern"C" int ckscl_(complex *zr, real *fnu, integer *n, complex *y, integer *nz, complex *rz, real *ascle, real *tol, real *elim);
extern"C" int cmlri_(complex *z__, real *fnu, integer *kode, integer *n,     complex *y, integer *nz, real *tol);
extern"C" int crati_(complex *z__, real *fnu, integer *n, complex *cy, real *tol);
extern"C" int cs1s2_(complex *zr, complex *s1, complex *s2, integer *nz, real *ascle, real *alim, integer *iuf);
extern"C" int cseri_(complex *z__, real *fnu, integer *kode, integer *n, complex *y, integer *nz, real *tol, real *elim, real *alim);
extern"C" int cshch_(complex *z__, complex *csh, complex *cch);
extern"C" int cuchk_(complex *y, integer *nz, real *ascle, real *tol);
extern"C" int cunhj_(complex *z__, real *fnu, integer *ipmtr, real *tol, complex *phi, complex *arg, complex *zeta1, complex *zeta2, complex\
 *asum, complex *bsum);
extern"C" int cuni1_(complex *z__, real *fnu, integer *kode, integer *n, complex *y, integer *nz, integer *nlast, real *fnul, real *tol     \
  , real *elim, real *alim);
extern"C" int cuni2_(complex *z__, real *fnu, integer *kode, integer *n, complex *y, integer *nz, integer *nlast, real *fnul, real *tol     \
  , real *elim, real *alim);
extern"C" int cunik_(complex *zr, real *fnu, integer *ikflg, integer *ipmtr, real *tol, integer *init, complex *phi, complex *zeta1, complex\
 *zeta2, complex *sum, complex *cwrk);
extern"C" int cuoik_(complex *z__, real *fnu, integer *kode, integer *ikflg,     integer *n, complex *y, integer *nuf, real *tol, real *elim, re\
al *alim);
extern"C" int cwrsk_(complex *zr, real *fnu, integer *kode, integer *n,     complex *y, integer *nz, complex *cw, real *tol, real *elim, real *a\
lim);
extern"C" real gamln_(real *z__, integer *ierr);
extern"C" integer i1mach_(integer *i__);
extern"C" real r1mach_(integer *i__);
extern"C" int xerror_(char *mess, integer *nmess, integer *l1, integer *l2,     ftnlen mess_len);
    #endif

in test_cbesj.cpp,

#include "mydemo.h"
#include "f2c.h"
#include <math.h>
#include <iostream>
#include <stdio.h>
#include <assert.h>
using namespace std;
int main(void)
{
  //  double x=86840.;                                                                                                                       
  //int nu=46431,j, err;                                                                                                                     
  complex *z,z__;
  z__.r = 3.0;z__.i = 2.0;z = &z__;
  cout << z->r << '\t' << z->i << endl;
  real *fnu;float fnu__ = 3.0;fnu = &fnu__;
  integer *kode ;long int kode__=1;kode=&kode__;
  integer *n    ;long int n__=1;n=&n__;
  complex *cy;
  integer *nz;
  integer *ierr;
  cbesj_(z, fnu, kode, n, cy, nz, ierr);
  cout << cy->r << '\t' << cy->i << endl;

  return 0;
}

Then,

g++ -c -g test_cbesj.cpp
g++ -o test *.o -lg2c
./test
3   2
Segmentation fault (core dumped)
gdb test 
GNU gdb (Ubuntu/Linaro 7.4-2012.04-0ubuntu2.1) 7.4-2012.04
Copyright (C) 2012 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "i686-linux-gnu".
For bug reporting instructions, please see:
<http://bugs.launchpad.net/gdb-linaro/>...
Reading symbols from /media/Downloads/amos-4/test...done.
(gdb) run
Starting program: /media/Downloads/amos-4/test 
3   2

Program received signal SIGSEGV, Segmentation fault.
0x0804b355 in cbesj_ ()
(gdb) frame 0
#0  0x0804b355 in cbesj_ ()
(gdb) frame 1
#1  0x0805a3ca in main () at test_cbesj.cpp:21
21    cbesj_(z, fnu, kode, n, cy, nz, ierr);

Thanks for roygvib's reply! Good suggestions actually. Here is the changed test_cbesj.cpp:

complex z, cy;
  float fnu;
  long int kode, n, nz, ierr;

  z.r = 3.0; z.i = 2.0;
  fnu = 3.0;
  n = 1; kode = 1;
  cout.precision(16);
  cbesj_( &z, &fnu, &kode, &n, &cy, &nz, &ierr );
  cout << cy.r << '\t' << cy.i << endl;
  cout << "nz=" << nz << endl;
  cout << "ierr=" << ierr << Lendl;

No more seg fault. But out of some reasons, the code does not work as expected:

./test
-1.343533039093018  -1.343533992767334
nz=0
ierr=4

and answers are wrong and ierr also says so from the source code:

C           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
C                    NZ= 0   , NORMAL RETURN
C                    NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO
C                              DUE TO UNDERFLOW, CY(I)=CMPLX(0.0,0.0),
C                              I = N-NZ+1,...,N
C           IERR   - ERROR FLAG
C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
C                    IERR=1, INPUT ERROR   - NO COMPUTATION
C                    IERR=2, OVERFLOW      - NO COMPUTATION, AIMAG(Z)
C                            TOO LARGE ON KODE=1
C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
C                            ACCURACY
C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
C                            CANCE BY ARGUMENT REDUCTION
C                    IERR=5, ERROR              - NO COMPUTATION,
C                            ALGORITHM TERMINATION CONDITION NOT MET

Answer 1

I downloaded cbesj (or zbesj) and related files from netlib, but for some reason they did not work with gfortran4.8 (on Linux x86_64). More precisely, cbesj always gives ierr=4, while I could not compile zbesj because zabs has too many arguments (but please see the update below).

---

So I gave up using the original AMOS codes and tried openspecfun that is also based on AMOS. I compiled the latter simply by typing make (with gfortran), which generated libopenspecfun.a etc. Then I made the following code to test zbesj:

#include <cstdio>

extern "C" {
    void zbesj_( double *zr, double *zi, double *fnu, int *kode, int *n,
                 double *Jr, double *Ji, int *nz, int *ierr );
}

int main()
{
    int    n, nz, ierr, kode;
    double fnu, zr, zi, Jr, Ji;

    fnu = 2.5;
    zr = 1.0; zi = 2.0;
    n = 1; kode = 1;

    zbesj_( &zr, &zi, &fnu, &kode, &n, &Jr, &Ji, &nz, &ierr );

    printf( "fnu = %20.15f\n", fnu );
    printf( "z   = %20.15f %20.15f\n", zr, zi );
    printf( "J   = %20.15f %20.15f\n", Jr, Ji );
    printf( "nz, ierr = %d %d\n", nz, ierr );

    return 0;
}

and compiled as

g++ test_zbesj.cpp libopenspecfun.a -lgfortran

which gives

fnu =    2.500000000000000
z   =    1.000000000000000    2.000000000000000
J   =   -0.394517225893988    0.297628229902939
nz, ierr = 0 0

Because this site says the answer is -0.394517...+ 0.297628...i, the result of zbesj seems correct.

---

[Update]

By reading the original code more carefully and also comparing with openspecfun, it was found that the user needs to uncomment appropriate sections of i1mach.f and r1mach.f (or d1mach.f) depending on the machine environment. For Linux x86_64, it seems sufficient to uncomment the section tagged with

C     MACHINE CONSTANTS FOR THE IBM PC

With this modification, the cbesj worked correctly with ierr=0; otherwise it gives ierr=4 because some constants default to 0. For the double-precision version, we also need to deal with the user-defined ZABS, whose definition conflicts with the intrinsic one. The Intel Fortran (ifort) recognizes the user-defined ZABS as such and compiles them successfully (although with a lot of warnings), while gfortran stops compiling with syntax error (assuming it to be the intrinsic one). openspecfunc avoids this problem by rewriting all the ZABS with the intrinsic one. Anyway, with the above modifications both cbesj and zbesj worked correctly (as expected).

---

[Update 2]

It turned out that using the BLAS version of d1mach.f, r1mach.f, and i1mach.f, things become even simpler because they determine machine-dependent constants automatically and we don't need to modify the codes manually. Please see the Netlib/FAQ page for details. The same trick applies to SLATEC also (e.g., this page).

wget http://www.netlib.org/blas/i1mach.f 
wget http://www.netlib.org/blas/r1mach.f
wget http://www.netlib.org/blas/d1mach.f