OctaveFortran

ScaryOctave > CategoryCode > CategoryExternal

---

This is a guide to incorporating Fortran code in Octave, by example. Octave is written in C++, so the approach taken here is to write C++ functions which call Fortran functions or subroutines, and then set up the C++ functions as `dynamically linked' functions for Octave. These can then be called from the Octave command line or in Octave functions just like any other Octave function.

First, a very simple example to get started. It won't be "hello, world", because strings are a bit tricky. Instead, say we have the Fortran subroutine

       SUBROUTINE NINE (R)
       DOUBLE PRECISION R
       R = 9D0
       END SUBROUTINE NINE

we'd write the C++ wrapper

 #include <octave/oct.h>
 #include "f77-fcn.h"
 extern "C"
 { 
   int F77_FUNC (nine, NINE) (double& R);
 }
 DEFUN_DLD (o9, args, ,
 	   "- Loadable Function: [NINE] = o9 ()\n\
      \n\
      Returns the number nine.")
 {
   octave_value_list retval;
   double r;
   F77_XFCN (nine, NINE, (r));
   if (f77_exception_encountered)
     {
       error ("unrecoverable error in o9");
       return retval;
     }
   retval(0) = octave_value (r);
   return retval;
 }

compile them (at the shell, not in Octave, the "$" being the shell's prompt) with

 $ mkoctfile o9.cc nine.f

assuming those are the names of the files containing the C++ wrapper and the Fortran subroutine. Then o9 is available in Octave like any other function:

 octave:1> o9 
 ans = 9

More details

Fortran has different calling conventions from C/C++, and conventions vary from system to system, and perhaps between compilers on the same system. A set of macros hides the details.

• Fortran is not case-sensitive but the linker is, so the compiler changes the name of the function to all uppercase or all lowercase. The F77_FUNC(name,NAME) macro selects the appropriate name given upper and lower case versions of the name. Use it instead of the name when declaring or calling the function.

• If you don't trust your fortran routine with all inputs you should call it with:

F77_XFCN (name, NAME, (args...));

Otherwise, if you are sure it will not generate any signals or run for excessive time, you can use:

F77_FUNC(name,NAME)(args...);

• Use of underscores in the Fortran subroutine name causes some problems. Some compilers add trailing underscores, some add none, some only add an underscore if there is not one there already. If your function name contains a trailing underscore, you need to use F77_FUNC_(name,NAME).

• Fortran subroutines may or may not return a value (functions always do of course). Use F77_RET_T as the return type for a Fortran subroutine in the declaration.

• In Fortran, all parameters are pass-by-reference, which means that you need to declare type& for scalar parameters and type* for vector parameters.

• Fortran strings have an intrinsic length, and unlike C, do not use a trailing \0 to determine the length of the string. When you pass a string or a string constant to a subroutine, you have to also pass the length. Some compilers treat strings as structures (so length and data are passed together) and others pass the length of the strings after every other parameter has been passed.

Rather than using "char*" for string argument declarations and one extra "int" declaration at the end for each string length, you need either

F77_CHAR_ARG_DECL

or

F77_CONST_CHAR_ARG_DECL

depending on the const-ness of the argument.

At the end of the argument list, leave the option for a string length declaration if that's what the compiler uses:

F77_CHAR_ARG_LEN_DECL

To call a function with a string argument, you need to use one of the following in place of the string argument:

F77_CHAR_ARG(x) for non-const arg, len = strlen(x)

F77_CONST_CHAR_ARG(x) for const arg, len = strlen(x)

F77_CXX_STRING_ARG(x) for const arg, len = x.lenth()

F77_CHAR_ARG2(x,n) for non-const arg, len = n

F77_CONST_CHAR_ARG2(x,n) for const arg, len = n

Regardless, you need to use the following macro at the end of the argument list for each string you pass. It may or may not be ignored.

F77_CHAR_ARG_LEN(n)

• There is also a macro F77_RETURN(retval) which is for callbacks to fortran functions.

Examples

There are a number of examples in the Octave sources to begin with, such as dozens of the C++ files in liboctave and a couple (balance.cc[1] and qz.cc[2]) in src/DLD-FUNCTIONS. In liboctave/CColVector.cc[3], for example, a call is made to the BLAS subroutine ZGEMV which performs matrix-vector multiplication.

It's possible to get the basic idea from these existing files, but let's look at a simple but complete example.

Say we had a Fortran subroutine TNINE:

           SUBROUTINE TNINE (IOPT, PARMOD, PS, X, Y, Z, BX, BY, BZ)
           INTEGER IOPT
           DOUBLE PRECISION PARMOD(10), PS, X, Y, Z, BX, BY, BZ    
     C     This is just a test subroutine body, to check connexions.
     C     Put the sum of PARMOD in PS, and X, Y, Z into BX, BY, BZ
           INTEGER I
           PS = 0D0
           DO 1 I=1, 10
              PS = PS + PARMOD (I)
      1    CONTINUE
           BX = X
           BY = Y
           BZ = Z
           END

A minimal C++ wrapper would look like:

    #include <octave/oct.h>
    #include "f77-fcn.h"
    extern "C"
     {
       int F77_FUNC (tnine, TNINE) (const int& IOPT, const double* PARMOD,
 				    double& PS,
 				    const double& X, const double& Y,
 				    const double &Z,
 				    double& BX, double& BY, double& BZ );
     }
     DEFUN_DLD (t96, args, ,
 	       "- Loadable Function: [PS, BX, BY, BZ] = t96 (PM, X, Y, Z)\n\
     \n\
     Returns the sum of PM in PS and X, Y, and Z in BX, BY, and BZ.")
     {
       octave_value_list retval;
       const int dummy_integer = 0;
       Matrix pm;
       const double x = args(1).double_value(), y = args(2).double_value(),
       z = args(3).double_value();
       double ps, bx, by, bz;
       pm = args(0).matrix_value ();
       F77_XFCN (tnine, TNINE,
 	        (dummy_integer, pm.fortran_vec(), ps, x, y, z, bx, by, bz) );
       if (f77_exception_encountered)
         {
           error ("unrecoverable error in t96");
           return retval;
         }
       retval(0) = octave_value (ps);
       retval(1) = octave_value (bx);
       retval(2) = octave_value (by);
       retval(3) = octave_value (bz);
       return retval;
     }

Compile this (in the Bourne Again Shell) with

    $ mkoctfile t96.cc tnine.f

and run it in Octave like:

    octave> [p, x, y, z] = t96 (1:10, sqrt (2), pi, e)
    p = 55
    x = 1.4142
    y = 3.1416
    z = 2.7183

The main thing to be added to this approach to extend it to more complicated Fortran subprograms is the way to handle different Fortran data types. For this, go to the OctaveFortranMex page or refer to the examples in the Octave sources cited above.

---

CategoryCode CategoryExternal