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"{[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]
...
Here is a Makefile for compiling and linking the code:
OBJECTS = dcsint.o main.o
main: ${OBJECTS}
g++ -O ${OBJECTS} -o main
The local installation of make should have default rules for
compiling the Fortran module.
[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 prototype extern "C" as in the example below. (Don't use this "extern" in a C program.)[2] Fortran types must match C++ types.
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 specify CHARACTER JOBZ, where JOBZ is 'N' or 'V'. In this case a C/C++ constant character string "N" or "V" works. Do not use single quotes, however, since Fortran needs the pointer that is generated with double quotes.Here 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, N are INTEGER types. Otherwise they are REAL (single precision). This default may be overridden with an IMPLICIT statement near the top of the code. For example IMPLICIT REAL*8 (A-H,O-Z) causes the default real type to be double precision REAL*8. An array variable declared with DIMENSION takes 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.[4] Fortran 77 arguments are pointers.
All Fortran arguments are passed as pointers, whether they are input or output values.[5] Fortran array indexing is base 1.
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.[6] Fortran presents multiple array subscripts backwards from C/C++
Fortran matrix element A(3,5) translates to C/C++ matrix element a[4][2]. (Subtract 1 for zero base indexing and reverse the order of the subscripts.)For a discussion of the historical background of language differences and further linkage issues, see Nelson Beebe's notes .The Fortran declaration REAL A(10,20) translates to the C/C++ declaration float 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
a[MAX][MAX] containing an n by n symmetric banded
matrix with kd bands above the diagonal, and you want to pack
it into a new matrix ab for feeding to the subprogram.
Write C/C++ code to pack the matrix, following the instructions above, using the 'U' convention.
double ab[MAX][LDAB];
where MAX and LDAB are predefined macros.
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<=j
for all i,j in the range 1,2,...,n.
Here is a straightforward solution:
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++.
This page is maintained by: Carleton DeTar Mail Form Last modified 27 October 2012