Next: Compiling and Linking with the Packages
Up: Scientific Subroutine Packages
Previous: Introduction

The Care and Feeding of a Fortran Subroutine

The five packages, IMSL, ESSL, LINPACK, EISPACK, and LAPACK are all written in Fortran, so require the Fortran calling procedure. If the calling program is written in Fortran, there should be no special difficulty in setting up the calling sequence. If you are writing C, the procedure is also very straightforward, if you follow simple guidelines. Here is an example, taken from the EISPACK library for the subroutine IMTQL2, which finds the eigenvalues of a symmetric tridiagonal matrix using the QL algorithm (quite similar to the QR algorithm). If you search the documentation file for the string SYMMETRIC TRIDIAGONAL you find the Fortran header

       SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR)
 C
       INTEGER I,J,K,L,M,N,II,NM,MML,IERR
       DOUBLE PRECISION D(N),E(N),Z(NM,N)
       DOUBLE PRECISION B,C,F,G,P,R,S,TST1,TST2,PYTHAG
Further along, you read that the vector D is the list of diagonal elements and E is the list of off-diagonal elements of the symmetric, tridiagonal matrix. The eigenvalues are returned in D and the eigenvectors in the two-dimensional matrix Z.

Fortran example

To set up the call in Fortran, you must match all arguments in position and type. The names can be different. Suppose your array has the diagonal elements in the vector A and the offdiagonal elements in B. You decide to use Z for the eigenvectors. You might set up the call as follows:

       INTEGER N,IERR
       PARAMETER (MAX1 = 20, MAX2 = 20)
       DOUBLE PRECISION A(MAX1),B(MAX1),Z(MAX1,MAX2)
           ...
           <load values into A and B>
           <load values into Z>
       CALL IMTQL2(MAX1,N,A,B,Z,IERR)
       IF(IERR.NE.0)THEN
           <error procedure>
       END IF
The first eigenvector is in Z(I,1), for I = 1,...,N. Notice that N could be variable, up to MAX1. It is essential that the storage dimension MAX1 be reported for NM and not the working dimension N.

C example, general considerations

To set up the call in C, you must also match all variables in position and type. Here are some important pointers for getting the C call to conform to the Fortran conventions:

  1. All library subprogram names are in lowercase.
  2. All arguments are passed by reference, rather than by value.
  3. The array subscript 1 in Fortran corresponds to the subscript 0 in C.
  4. In two and higher dimensional arrays, the order of the subscripts is reversed. Thus in C, a[i][j] corresponds to Fortran A(j+1,i+1).

C example: zero array offset

Taking the same example as in Fortran above, suppose your array has the diagonal elements in the vector a and the offdiagonal elements in b. You decide to use z for the eigenvectors. If you are using the convention that C arrays have zero subscript offset, i.e., your first vector element is a[0], etc, you might set up the call as follows:

/*         Example using zero array subscript offset   */
#define MAX1 20
#define MAX2 20
           ...
       int n, ierr, nm = MAX2;
       double a[MAX1],b[MAX1],z[MAX1][MAX2];
       void imtql2(int *nm, int *n, double a[], double b[], 
             double z[][], int *ierr);
           ...
           <load values into a and b:>
           <Note that a[0] is Fortran A(1), etc.>
           <load values into z>

       imtql2(&nm,&n,a,b,z,&ierr);
       if(ierr<>0)
          {
           <error procedure>
          }
The first eigenvector is in z[0][i] for i = 1...,n. Note that the scalars nm, n, ierr are listed by address, through the use of the symbol &. The array references are automatically passed as addresses according to C convention. Thus the argument a is equivalent to the argument &a[0]. Note also that because of the subscript reversal, the array dimension NM required by the subprogram must be the second dimension MAX2 in the C declaration, where it is first in the Fortran DIMENSION statement. That is why we have initialized nm to MAX2. (Of course, in this example it is natural for both dimensions to be the same.)

C example: zero array offset

If you are using a unit offset for your arrays, i.e., your first vector element is a[1], etc, you would do it this way:

/*         Example using unit array subscript offset   */
#define MAX1 21
#define MAX2 21
           ...
       int n, ierr, nm = MAX2;
       double a[MAX1],b[MAX1],z[MAX1][MAX2];
       void imtql2(int *nm, int *n, double a[], double b[], 
             double z[][], int *ierr);
           ...
           <load values into a and b:>
           <Note that a[1] is Fortran A(1), etc.>
           <load values into z>

       imtql2(&nm,&n,&a[1],&b[1],&z[1][1],&ierr);
       if(ierr<>0)
          {
           <error procedure>
          }
By passing the addresses of the first elements of the arrays, we establish the correspondences a[i] = A(i), z[i][j] = Z(j,i), etc, for all the elements of the arrays that we actually use. The unused elements z[i][0] for i > 1 are harmlessly interpreted as the unused elements Z(MAX2,i-1).

The ESSL library has a C header essl.h that modifies the calling sequence for the routines in the ESSL library, mainly by dispensing with ampersands in front of some of the scalars. The C calling sequence is given in the documentation. To use the documented C style calling sequence, you must #include <essl.h> as a header to your code.


Next: Compiling and Linking with the Packages
Up: Scientific Subroutine Packages
Previous: Introduction

Carleton DeTar