NSC: NAG C Header Files for the SGI

Calling NAG Fortran Library Routines from C Language Programs
Using the NAG C Header Files

Ian Hounam

NAG Ltd, Oxford

© The Numerical Algorithms Group Ltd, Oxford UK. 1998

19 March 1998

Modified by Bo Einarsson 2001-08-10 and 2002-02-28 for the SGI 3800

These header files are for Mark 18 but work also with Mark 19 for all routines which are in both versions. If you use one of the 62 new routines only available in Mark 19 or one of the 95 new routines in Mark 20 contact Bo Einarsson for assistance.

1 Introduction

A great number of systems allow the C programmer to call other language routines. Indeed the ANSI standard definition of C provides a powerful argument checking facility that, given the correct definition of function prototypes, can facilitate cross language communication between C and, say, Fortran.

A header file containing the function prototypes can be included in the user's program to allow the C compiler to check argument passage. Such a header file has been created for the current NAG Fortran Library and the current NAG Graphics Library. This was done automatically from the Library source code in order to ensure its correctness. This document explains how to call Fortran routines from C using the NAG header file.

This document describes how this is commonly achieved in Unix systems.

2 Argument Types

The table below lists the mapping between relevant Fortran and C argument and function types.
   DOUBLE PRECISION D       double d;
   INTEGER I                long int i;
   LOGICAL L                long int l;
   CHARACTER*n S            char s[n];
   DOUBLE COMPLEX Z         struct {double re,im;} z;
All arguments are passed as pointers to a variable; this means that a constant must first be assigned to a variable and then the address of the variable passed.

As Fortran stores multi-dimension arrays in column major order whereas C stores in row major order, either

  1. the C routine must store and manipulate the transpose of the problem matrix, or
  2. the C routine must transpose the matrix before and after calling the Fortran routine.
Fortran CHARACTER type is different to C null terminated strings so that the actual length of the string must be passed as an extra argument. For example the Fortran subroutine -
         SUBROUTINE X00XXF(D1,S,D2)
where the arguments are declared -
         DOUBLE PRECISION D1,D2
         CHARACTER*5 S
would, with an f77 compatible compiler, be called from C thus -
         x00xxf_(&d1,s,&d2,length_of_s);
where the arguments are declared -
         double d1,d2;
         char s[5];
         int length_of_s = 5;
(length_of_s is the length of the string s)

If there is more than one string argument, further extra length arguments are appended, in order, to the end of the argument list.

The Fortran type DOUBLE COMPLEX (or COMPLEX*16) is provided in the NAG header files by the typedef "Complex", which expands to "struct {double re,im;}"

3 Subroutines and Functions

Fortran subroutines are declared as void functions in C.

Procedure arguments, i.e. function or subroutine names, are passed by address in the normal C manner.

4 Complex functions and Silicon Graphics

With the native Fortran compilers on DEC Alpha Unix, IBM AIX and Silicon Graphics, the C programmer wishing to access Fortran DOUBLE COMPLEX FUNCTIONs has real problems. The conventions used by the native Fortran compilers on these systems do not allow the C programmer to access the return value of DOUBLE COMPLEX FUNCTIONs directly because the DOUBLE COMPLEX return value is stored in registers that do not map to any C data type.

As it is impossible to access the return value of a DOUBLE COMPLEX FUNCTION from C, another solution is to write a Fortran "jacket" routine to convert the DOUBLE COMPLEX FUNCTION to a SUBROUTINE with an extra initial argument containing a pointer to the return value. The "jacket" routine then calls S01EAF.

For example to call S01EAF, write a jacket routine called S01EAFJ which would have the following prototype.

extern void s01eafj_(Complex *ret_val, CONST Complex *z,
                     int *ifail);
The jacket routine is then called from the C program rather than S01EAF.
      DOUBLE COMPLEX FUNCTION S01EAFJ(RET_VAL, Z, IFAIL)
      DOUBLE COMPLEX RET_VAL, Z, S01EAF
      INTEGER IFAIL
      RET_VAL = S01EAF(Z, IFAIL)
      END
This routine can be compiled and linked with the C files and NAG Fortran Library in the normal way.

Similar jacket routines can be written for user supplied functions.

The supplied C Header files do not contain prototypes for these jacket functions. Use the provided version of the C Header File and modify the prototypes appropriately.

5 Character functions

A character function is implemented as a void function with two extra arguments added to the start of the argument list, the first is a pointer to the returned string and the second its length.

e.g.

      CHARACTER*10 FUNCTION F(I)
is called thus:
     extern void f_(char *,int,int *);
      char c[10];
      int clen = 10,i;
      f_(&c,clen,&i)

6 Routine Names

Some Unix systems (including SGI) require the addition of an underscore to the name of Fortran routines called from C. e.g. X01AAF becomes x01aaf_.

7 Multi-Dimension Arrays

There is no syntax in the C language to specify an arbitrary array dimension and it is not possible to know in advance the dimensions of the arrays in your program. The best that the C Header Files can do is to declare the multi-dimensional arrays as one dimension arrays with a comment stating the actual number of dimensions, e.g. for 2 dimensions:
int array[] /* 2 dimension */

and for 3 dimensions:

double array[] /* 3 dimension */
The prototype for a hypothetical NAG Fortran routine with a 2 dimensional DOUBLE PRECISION array argument would look like this:
extern void nagsub_(double[] /* 2 dimensions */);
A simple program to call this routine might look like this:
main ()
{
  double p[2][2];
  nagsub_((double *)p);
}
Note that we need to cast the 2 dimensional C array actual argument to (double *).

The example prototype below shows how to call a hypothetical NAG routine that takes a single subroutine argument. This user supplied subroutine takes a 2 dimension DOUBLE PRECISION array and an integer which specifies the leading dimension of the Fortran array.

extern void nagsub_(void (*f) (double[] /* 2 dimension */, int *));

The C code for the user supplied function is listed below. The 2 dimension array is passed as a pointer to double and the code must carry out array indexing using the dimension information passed from Fortran. In this case, the macro P uses the leading dimension of the Fortran array, which is the trailing dimension of the C array, to index into the array p. The array p is only referenced through this macro.

void fun(double p[], int *tdp)
{
#define P(I,J) p[(I)*(*tdp) + (J)]
  P(0,0) = 0.0;
  P(1,1) = 1.0;
}
The main function looks like this:
main ()
{
  void fun(double p[], int *tdp);
  nagsub_(fun);
}
Example 3 below shows a complete example program that illustrates these concepts.

10 Linking your program

It is always necessary to link your program with the Fortran run time library routines. This can most easily be achieved on most Unix systems by linking using the f77 command, not cc or ld. For example,
cc -c myprog.c
f90 -o myprog myprog.o -lnag -L'/opt/nag/flsg618db/lib32/mips4/'

12 Examples

All the examples below are written in ANSI standard C; minor changes to remove function prototypes will be necessary for non-ANSI compilers.

All three examples below are available as downloadable files.

12.1 Example 1

This example shows how to call C05AJF using a procedure parameter, the function "f". The example is a translation of the Fortran Library example program into C. This example applies to both standard Unix and VMS, but on some platforms the name c05ajf does not have an appended underscore. No changes are required on SGI!
#include <math.h>
#include <stdio.h>
#include <nagmk18.h>

main()
{
  double eps, eta, x;
  int ifail, k, nfmax;
  double f(double *);

  printf("C05AJF Example Program Results\n\n");
  for (k=1; k<=2; k++)
    {
      eps = k==1 ? 0.1e-3 : 0.1e-4;
      x = 1.0;
      eta = 0.0;
      nfmax = 200;
      ifail = 1;
      /* Remove underscore as appropriate */
      c05ajf_(&x,&eps,&eta,f,&nfmax,&ifail);
      if (ifail==0)
        printf("With eps = %e   root = %f\n",eps,x);
      else
        {
          printf("ifail = %d\n",ifail);
          if (ifail==3 | ifail==4)
            printf("With eps = %e final value = %f\n",eps,x);
        }
     }
}
double f(double *x)
{
  return exp(-*x) - *x;
}

12.2 Example 2

This is a simplified version of the NAG Fortran Library example program for F01QDF. Array sizes have been chosen to be the same as the problem size to illustrate better the use of two dimensional arrays.

#include <stdio.h>
#include <nagmk18.h>

/* Simplified example program for F01QDF */

main()
{
  int i,ifail,j,m = 5,n = 3,ncolb = 2;
  int lda = m,ldb = m;
  /* Initialise arrays in column major order */
  static double a[3][5] =
    {
      2.0, 2.0, 1.6, 2.0, 1.2,
      2.5, 2.5, -0.4, -0.5, -0.3,
      2.5, 2.5, 2.8, 0.5, -2.9
    };
  static double b[2][5] =
    {
      1.1, 0.9, 0.6, 0.0, -0.8,
      0.0, 0.0, 1.32, 1.1, -0.26
    };
  double work[2],zeta[3];

  printf("F01QDF Example Program Results\n\n");
  ifail = 0;
  f01qcf_(&m,&n,(double *)a,&lda,zeta,&ifail);
  ifail = 0;
  f01qdf_("Transpose","Separate",&m,&n,(double *)a,&lda,zeta,&ncolb,
          (double *)b,&ldb,work,&ifail,9,8);
  printf("Matrix  Q'*B\n");
  for (i=0; i<m; i++)
    {
      for (j=0; j<ncolb; j++)
        {
          printf("%f ",b[j][i]);
        }
      printf("\n");
    }
}

12.3 Example 3

This example illustrates the use of multi-dimension arrays in user supplied functions. It is a translation of the NAG Fortran Library example program for D03PCF. Please refer to the NAG Fortran Library documentation for details of this program. Note that this example is only shown here in the "Standard Unix" version. Please make the appropriate changes for other platforms.
/* D03PCF Example Program
   C Version. NAG Copyright 1997. */

#include <stdio.h>
#include <math.h>

/* Include the appropriate NAG C Header File */
#include <nagmk18.h>

/* Global data */
double alpha;

/* Routine for PDE initial condition */
void uinit(double u[20][2], double *x, int npts)
{
  int i;
  for (i=0; i<npts; i++)
    {
      u[i][0] = 2.0*alpha*x[i];
      u[i][1] = 1.0;
    }
}

/* Define the system of PDEs */
void pdedef(int *npde, double *t, double *x, double *u, double dudx[],
            double p[], double q[], double r[], int *ires)
{
/* Macro to do subscripting into 2 dimensional array p, which
   is passed as double p[]. */
#define P(I,J) p[(I)*(*npde) + (J)]
  q[0] = 4.0*alpha*(u[1] + *x * dudx[1]);
  q[1] = 0.0;
  r[0] = *x * dudx[0];
  r[1] = dudx[1] - u[0]*u[1];
  P(0,0) = 0.0;
  P(0,1) = 0.0;
  P(1,0) = 0.0;
  P(1,1) = 1.0 - (*x)*(*x);
}

void bndary(int *npde, double *t, double u[], double *ux, int *ibnd,
            double beta[], double gamma[], int *ires)
{
  if (*ibnd==0)
    {
      beta[0] = 0.0;
      beta[1] = 1.0;
      gamma[0] = u[0];
      gamma[1] = -u[0]*u[1];
    }
  else
    {
      beta[0] = 1.0;
      beta[1] = 0.0;
      gamma[0] = -u[0];
      gamma[1] = u[1];
    }
}

main()
{
  const int npde=2, npts=20, intpts=6, itype=1, neqn=40, niw=64, nw=1128;
  const double acc = 1.0e-4;
  double hx, pi, piby2, tout, ts;
  const int m = 1, itask = 1, itrace = 0;
  int i, ifail, ind, it, j;
  /* u[npts][npde], uout[itype][intpts][npde], w[nw], x[npts] */
  double u[20][2], uout[1][6][2], w[1128], x[20];
  /* xout[intpts] */
  double xout[] = {0.0,0.4,0.6,0.8,0.9,1.0};
  int iw[64];
  void bndary(int *npde, double *t, double *u, double *ux, int *ibnd,
              double *beta, double *gamma, int *ires);
  void pdedef(int *npde, double *t, double *x, double *u, double *dudx,
              double p[], double *q, double *r, int *ires);
  void uinit(double u[20][2], double *x, int npts);

  printf("D03PCF Example Program Results\n\n");
  alpha = 1.0;
  ind = 0;

  /* Set spatial mesh points */

  piby2 = 0.5*x01aaf_(&pi);
  hx = piby2/(npts-1);
  x[0] = 0.0;
  x[npts-1] = 1.0;
  for (i=1; i<=npts-2; i++)
    x[i] = sin(hx*i);

  /* Set initial conditions */

  ts = 0.0;
  tout = 0.1e-4;
  printf("Accuracy requirement = %lf\n",acc);
  printf("Parameter alpha = %lf\n\n",alpha);
  printf("  T  /  X  ");
  for (i=0; i<=5; i++)
    printf("%8.4lf",xout[i]);
  printf("\n");

  /* Set the initial values */

  uinit(u,x,npts);
  for (i=0; i<=4; i++)
    {
      ifail = -1;
      tout = 10.0*tout;
      /* Cast 2-D array u to double * to match prototype */
      d03pcf_(&npde,&m,&ts,&tout,pdedef,bndary,(double *)u,&npts,x,&acc,w,
              &nw,iw,&niw,&itask,&itrace,&ind,&ifail);

      /* Interpolate at required spatial points */

      /* Cast 2-D array u and 3-D array uout
         to double * to match prototypes */
      d03pzf_(&npde,&m,(double *)u,&npts,x,xout,&intpts,&itype,
              (double *)uout,&ifail);
      printf("%6.4lf U(1)",tout);
      for (j=0; j<intpts; j++)
        printf("%8.4lf",uout[0][j][0]);
      printf("\n       U(2)");
      for (j=0; j<intpts; j++)
        printf("%8.4lf",uout[0][j][1]);
      printf("\n");
    }
  printf("\n");

  /* Print integration statistics */

  printf("Number of integration steps in time                    %d\n",iw[0]);
  printf("Number of residual evaluations of resulting ODE system %d\n",iw[1]);
  printf("Number of Jacobian evaluations                         %d\n",iw[2]);
  printf("Number of iterations of nonlinear solver               %d\n",iw[4]);
}
© The Numerical Algorithms Group Ltd, Oxford UK. 1998