• Fortran External Calls
  • Introduction
• C External Calls

Wrapperless External Calls to C and Fortran Routines

by Waterloo Maple

Fortran External Calls

Introduction

The Fortran external API (given in libmaplefortran.so) provides routines that translate Maple data type to Fortran data type and vice versa. This allows you to invoke C and Fortran routines from the Maple environment, including from within your Maple procedures. Here are the steps to create C and Fortran external calls from Maple (Fortran 77 is assumed). Note that for Maple 7 and higher, no wrapper is needed.   Maple will do all the data type translations automatically and hence eliminates the need for a wrapper program. This method is cleaner and requires less overhead than the wrapper method because array input does not have to be copied to a "container".

As an example, the following Maple worksheet will invoke a Fortran external call fftwrap  that computes the FFT of a rectangular real signal. FFT is the fundamental operation of signal processing applications. FFT transforms a signal from the time domain to its spectrum in the frequency domain. The Fortran routines to compute the FFT are taken from SLATEC Common Mathematical Library, Version 4.1. Their source code can be downloaded from the website http://netlib2.cs.utk.edu/slatec  

Below is the Fortran routine realfft  (the source file is realfft.f). It is the front end to the FFT routines of the SLATEC package.

C***FFT routine - it is the front end to the routines of the SLATEC library
C X  - real input signal
C RY - real part of FFT(X)
C IY - imaginary part of FFT(X)
C N  - length of input signal
      SUBROUTINE realfft( X, RY, IY, N )
      REAL X(*)
      DOUBLEPRECISION RY(*)
      DOUBLEPRECISION IY(*)
      INTEGER N

      PARAMETER( MAXFFT = 4096 )
      REAL W( 2*MAXFFT + 15 )

C Call the SLATEC routines to compute the FFT of the real signal X
      CALL RFFTI( N, W )
      CALL RFFTF( N, X, W )

C Since the input is a real signal, the FFT routine only computes the first
C half of the spectrum. The other half can be obtained by taking the
C conjugate of the first half.

C Real part of the FFT
      RY(1) = X(1)
      DO 10, I = 1, N/2-1
          RY(I+1) = X(2*I)
          RY(N-I+1) = X(2*I)
10    CONTINUE
      RY(N/2+1) = X(N)

C Imaginary part of the FFT
      IY(1) = 0
      DO 20, I = 1, N/2-1
          IY(I+1) = X(2*I+1)
          IY(N-I+1) = -X(2*I+1)
20    CONTINUE
      IY(N/2+1) = 0    

      END

 

>	restart: N := 512:         # FFT length PULSE_LEN := 16:   # Length of rectangular input signal # Build rectangular input signal X := Vector(N, datatype=float[4]): for i from 1 to PULSE_LEN/2 do     X[i] := 1:     X[N-i+1] := 1: end do:

>	# Centralize and plot input signal Xplot := [seq(X[i], i=1..N)]: Xplot := [op(Xplot[N/2+1..N]), op(Xplot[1..N/2])]: T := [seq(i, i=1..N)]: Xplot := zip((x,y)->[x,y], T, Xplot): plot(Xplot);

>	# Real and Imaginary part of FFT(X) RY := Vector(N, datatype=float[8]): IY := Vector(N, datatype=float[8]):

>

#### EXTERNAL CALL TO FIND FFT OF X #### realfft := define_external( 'realfft_', AUTO, FORTRAN,                             x::ARRAY(float[4]),                             ry::ARRAY(float[8]),                             iy::ARRAY(float[8]),                             n::integer[4],                             LIB="librealfft.so" ): realfft( X, RY, IY, N ):

>	# Centralize and plot the magnitude of FFT(X) Yplot := zip( (u,v)->sqrt(u^2+v^2), RY, IY ): Yplot := [seq(Yplot[i], i=1..N)]: Yplot := [op(Yplot[N/2+1..N]), op(Yplot[1..N/2])]: F := [seq(i, i=1..N)]: Yplot := zip((x,y)->[x,y], F, Yplot): plot(Yplot);

>

C External Calls

Below is the code of  a C routine matmult  to compute matrix multiplication of two real matrices (the source file is matmult.c)

void matmult( double *A, double *B, double *C, int I, int J, int K )
{
    int i, j, k;
    double t;

    for( i = 0; i < I; ++i )
    {
        for( k = 0; k < K; ++k )
        {
            t = 0.0;
            for( j = 0; j < J; ++j )
            {
                t += A[i*J+j] * B[j*K+k];
            }
            C[i*K+k] = t;
        }
    }
}

The following Maple code will directly invoke the C routine matmult  to compute matrix multiplication.

>	A := Matrix(2, 2, datatype=float[8], order=C_order ): A[1,1]:=10 : A[1,2]:=20: A[2,1]:=30 : A[2,2]:=40: A;

>	B := Matrix(2, 3, datatype=float[8], order=C_order ): B[1,1]:=1.5 : B[1,2]:=2.5 : B[1,3]:=3.5: B[2,1]:=4.5 : B[2,2]:=5.5 : B[2,3]:=6.5: B;

The empty matrix that will be filled with the values of AB

>	C := Matrix(2, 3, datatype=float[8], order=C_order ): C;

>	matmult := define_external('matmult',AUTO,                   a::ARRAY(float[8]),                   b::ARRAY(float[8]),                   c::ARRAY(float[8]),                   i::integer[4],                   j::integer[4],                   k::integer[4],                   LIB="libmatmult.so"):

>	matmult(A, B, C, 2, 2, 3); C;

>