Linking with Fortran encoded subroutines

A great deal of legacy code used in computational science and engineering is written in Fortran, especially in the form of subroutine packages. With a few simple rules one can call Fortran encoded subroutines. A detailed example is given in notes on Fortran binding. The main consideration is that Fortran matrices are stored by columns with both subscripts starting at 1:

   A(1,1) A(2,1) ... A(50,1) A(1,2) A(2,2) ... A(1,50) ... A(50,50)

As long as the precision is the same and we match up the beginning of the array in C++ and Fortran we establish an equivalence between a[0][0] and A(1,1), a[0][1] and A(2,1), etc. We just have to remember that the subscripts are reversed and our subscript numbering is one less than Fortran's. Another consideration is that Fortran 77 (a very common version for legacy Fortran code) expects all subprogram parameters (arguments) to be pointers. Of course, we must be sure to match the Fortran precision for all numeric types. And finally, unlike C++, Fortran allows the array dimensions in formal subprogram parameters to be variable, so this dimension is often passed as a parameter as we now illustrate.

Suppose our mulmatvec subroutine were written in Fortran, and the author adopted the usual convention that the first Fortran matrix subscript is the row index and the second, the column. You might find the following description of this subroutine

      SUBROUTINE MULMATVEC(M, N, MAX, A, X, B)
      REAL*8 A(MAX,*), X(*), B(*)
      INTEGER M, N, MAX
C
C  INPUT PARAMETERS
C
C  M     - NUMBER OF ROWS IN A AND NUMBER OF COMPONENTS IN B
C  N     - NUMBER OF COLUMNS IN A AND NUMBER OF COMPONENTS IN X
C  MAX   - LEADING DIMENSION OF A
C  A     - REAL*8 TWO-DIMENSIONAL ARRAY CONTAINING THE MATRIX
C  X     - REAL*8 ARRAY CONTAINING THE VECTOR
C
C  OUTPUT PARAMETERS
C
C  B     - REAL*8 ARRAY CONTAINING THE RESULT OF MATRIX A TIMES VECTOR X

Here is how we would modify our code to call this subroutine. The actual Fortran subroutine is in a separate file, of course. We don't give the Fortran code here.

#include <iostream>
using namespace std;
#define MAX 50

extern "C"{
  void mulmatvec_(int *m, int *n, int *max, 
      double a[MAX][MAX], double x[MAX], double b[MAX]);
}

int main(){
  int i,j,m,n;
  double a[MAX][MAX], x[MAX], b[MAX];

  cout << "Enter the number of rows in the matrix \n";
  cin >> m;
  cout << "Enter the number of columns in the matrix \n";
  cin >> n;

  cout << "Enter the matrix by rows\n";
  for(i = 0; i < m; i++)
    for(j = 0; j < n; j++)
      cin >> a[j][i];

  cout << "Enter the vector\n";
  for(j = 0; j < n; j++)
    cin >> x[j];

  int max = MAX;
  mulmatvec_(&m,&n,&max,a,x,b);

  cout << "\nA*x = \n";
  for(i = 0; i < m; i++)
    cout << b[i] << "\n";
}

Let's examine some notable features of this code.

extern "C"{
  void mulmatvec_(int *m, int *n, int *max, 
      double a[MAX][MAX], double x[MAX], double b[MAX]);
}

A prototype is required. The extern "C" declaration tells the compiler not to modify the function name. The extra underscore character in the name may be needed, because Fortran compilers often add them. The Fortran type REAL*8 is equivalent to C++ double.

  for(i = 0; i < m; i++)
    for(j = 0; j < n; j++)
      cin >> a[j][i];

The order of subscripts has been reversed.

  int max = MAX;
  mulmatvec_(&m,&n,&max,a,x,b);

All parameters are pointers. Pointer have to point to a location in memory where the number is kept. Since MAX is not a variable, but just a macro, we have to store its value somewhere and pass a pointer to that location. That is the reason for introducing the variable max. Another way to do it is to promote MAX to a variable by declaring it to be a global constant. This is done by replacing the #define directive with the C++ declaration

const int MAX = 50;

Declaring it to be constant assures that the compiler won't let us accidentally change it. Then we can pass &MAX as a parameter, as long as we also change int * to const int * in the prototype for that parameter.

For another example and more details, please see the supplementary notes on Fortran binding.