It may happen that we find a useful subroutine written in ANSI Fortran 77 and want to call it from a C++ program.

Let's start with an illustration. Follow the links to the explanatory notes.

On netlib we find a subprogram by C.S. Duris [published reference ACM
TOMS 6, 92-103 (1980)] for constructing a cubic spline. We put it in a
file **dcsint.f**. At the top of the Fortran file we find
declarations and some related information:

SUBROUTINE DCSINT(IENT, H, N, TNODE, G, END1, ENDN, B, C, D) DIMENSION TNODE(N), G(N), B(N), C(N), D(N) C C INPUT PARAMETERS (NONE OF THESE PARAMETERS C ARE CHANGED BY THIS SUBROUTINE.) C C IENT - SPECIFIES END CONDITIONS WHICH ARE IN EFFECT. C H - THE STEP SIZE USED FOR THE DISCRETE CUBIC SPLINE. C N - NUMBER OF NODES (TNODE) AND DATA VALUES (G). C (N.GE.2) C TNODE - REAL ARRAY CONTAINING THE NODES (TNODE(I).LT. C TNODE(I+1)). C G - REAL ARRAY CONTAINING THE INTERPOLATING DATA. C END1 - END CONDITION VALUE AT TNODE(1). C ENDN - END CONDITION VALUE AT TNODE(N). C C OUTPUT PARAMETERS C C B - REAL ARRAY CONTAINING COEFFICIENTS OF C (T-TNODE(I)),I=1,2,....,N-1. C C - REAL ARRAY CONTAINING COEFFICIENTS OF C (T-TNODE(I))**2,I=1,2,...,N-1. C D - REAL ARRAY CONTAINING COEFFICIENTS OF C (T-TNODE(I))**3,I=1,2,...,N-1. ...Here is how we construct The C++ code for calling this routine:

extern "C"{Here is a Makefile for compiling and linking the code:^{[1]}void dcsint_(int *ient, float *h, int *n, float tnode[], float g[], float *end1, float *endn, float b[], float c[], float d[]);^{[3][4]}} #define MAX 100 ... int main() { int^{[2]}ient = 2; // End conditions (natural spline) float^{[2]}h; // Step size int n; // Number of node values float tnode[MAX]; // Spline nodes x[i] float g[MAX]; // Data values y[i] float end1 = 0; // End condition at node i = 0 (natural spline)^{[5]}float endn = 0; // End condition at node i = n-1 (natural spline)^{[5]}float b[MAX]; // Coefficient of (t - tnode[i]) i = 0,1, ..., n-2^{[5]}float c[MAX]; // Coefficient of (t - tnode[i])^2 i = 0,1, ..., n-2^{[5]}float d[MAX]; // Coefficient of (t - tnode[i])^3 i = 0,1, ..., n-2^{[5]}... // Load tnode[i] and g[i] for i = 0,1, ..., n-1^{[5]}... h = tnode[1] - tnode[0]; // Note that array names g, b, c, and d are already pointers to floats // but we use & to get pointers to ient, h, n, end1, endn dcsint_(&ient, &h, &n, tnode, g, &end1, &endn, b, c, d);^{[3]}^{[4]}...

OBJECTS = dcsint.o main.o main: ${OBJECTS} g++ -O ${OBJECTS} -o mainThe local installation of

**
[1] C++ prototype declarations must be extern "C".
**

You must include a prototype declaration of the Fortran routine you use. In C++ you must suppress "name mangling" by making the prototypeextern "C"as in the example below. (Don't use this "extern" in a C program.)

The variable type in C and C++ must match the Fortran type. Once one has determined the Fortran type, the translation to C and C++ is straightforward, but may depend on local hardware and compiler conventions.A number of LAPACK routines use single characters to specify options. For example the abstract might specifyHere is how to determine the Fortran variable type. All declarations for Fortran variables are found at the top of the code. Fortran variable types are either implicit (undeclared) or explicit (declared). Implicit types are set by the first letter of the variable name. By default variables beginning with

I, J, K, L, M, NareINTEGERtypes. Otherwise they areREAL(single precision). This default may be overridden with anIMPLICITstatement near the top of the code. For exampleIMPLICIT REAL*8 (A-H,O-Z)causes the default real type to be double precisionREAL*8. An array variable declared withDIMENSIONtakes on the implicit type. Any explicit type declaration overrides the implicit type. On our Solaris and Linux x86 systems, translate the most common numeric types as follows:

REAL, REAL*4 float REAL*8, DOUBLE PRECISION double INTEGER int COMPLEX, COMPLEX*8 float[2] (an array of two floats) COMPLEX*16 double[2] (an array of two doubles)

The Fortran compiled object module is usually given a lower case name. With many Fortran compilers by default, an extra suffix underscore character is also attached. If you find this happening, you have to add it to the name you call, as we did, or change the Fortran compilation options, if possible.

All Fortran arguments are passed as pointers, whether they are input or output values.

An array in C and C++ is indexed starting from 0. In Fortran the index starts at 1. This is not a problem for passing simple arrays, because we just hand the Fortran routine the address of the first element, which we think of as element 0 and Fortran thinks of as element 1. But the distinction is important in interpreting the documentation.

Fortran matrix elementFor a discussion of the historical background of language differences and further linkage issues, see Nelson Beebe's notes .A(3,5)translates to C/C++ matrix elementa[4][2]. (Subtract 1 for zero base indexing and reverse the order of the subscripts.)The Fortran declaration

REAL A(10,20)translates to the C/C++ declarationfloat a[20][10]. Note that Fortran documentation referring to the "leading" dimension of an array should be translated as the "last" dimension of an array: 10 in this example.

N (input) INTEGER The order of the matrix A. N >= 0. KD (input) INTEGER The number of superdiagonals of the matrix A if UPLO = 'U', or the number of subdiagonals if UPLO = 'L'. KD >= 0. AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i- j,j) = A(i,j) for j<=i<=min(n,j+kd).Suppose that in your C/C++ code you created a matrix

Write C/C++ code to pack the matrix, following the instructions above,
using the **'U'** convention.

double ab[MAX][LDAB];where

Next we translate the packing instructions. For the requested
**'U'** case we must apply the rule

AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=jfor all

for ( j = 1; j <= n; j++ ) for ( i = max(1,j-kd); i <= j; i++ ) ab[j-1][kd+i-j] = a[j-1][i-1];Note that we reversed the order of the subscripts and reduced each subscript by one to get the zero base indexing in C/C++.

Physics Department Home Page

*This page is maintained by:
*