MATLAB - Fortran interface

Using a Fortran callback in the MATLAB ODE solvers

MATLAB provides multiple external language interfaces. Mathworks recommends using the C++ interface whenever possible (after preparing this post I found out the C interface might still be the better option). If you are working with MATLAB R2017b or earlier, you can only use the C Matrix API.

The C++ MEX API allows MATLAB users to create MATLAB functions in C++ that integrate seamlessly with MATLAB. The API is built on top of the matlab::data:Array class defined in the MATLAB Data API. Since C++ allows calling functions written in C, and Fortran provides standard tools for interoperability with C, we can use the C++ MEX API to call a procedure in Fortran.

We can break the process of interfacing Fortran and MATLAB into four steps:

  1. Writing the Fortran procedure
  2. Defining the C prototype of the Fortran procedure
  3. Wrapping the Fortran procedure in a C++ MEX function
  4. Compilation and linking

Step 3 will be the hardest part as it requires good knowledge of C++ and some familiarity with the C++ MEX API. It took me circa two hours to produce a working example.

Writing the Fortran procedure

As stated in the title our goal is to solve an ODE defined by a Fortran procedure. We will reproduce the van der Pol equation example from MATLAB’s ode23s solver documentation. The van der Pol ODE’s in Fortran are defined by the procedure f_vanderpol:

! vanderpol_mod.f90

module vanderpol_mod

  use iso_c_binding, only: c_int, c_double
  implicit none

contains

  subroutine f_vanderpol(n,t,y,dydt) bind(c,name="f_vanderpol")
    integer(c_int), intent(in) :: n
      !! Number of ODE's
    real(c_double), intent(in) :: t
      !! Value of independent variable
    real(c_double), intent(in) :: y(n)
      !! Values of dependent variables y_i = y_i(t), i = 1,...,N
    real(c_double), intent(out) :: dydt(n)
      !! RHS of the system of ODE's

    dydt(1) = y(2)
    dydt(2) = 1000*(1 - y(1)**2)*y(2) - y(1)
  end subroutine

end module

Defining the C prototype of the Fortran procedure

This is the easiest step of the entire process. We just create a header file containing the C prototype of the Fortran function:

// vanderpol_mod.h

#ifdef __cplusplus
extern "C" {
#endif

void f_vanderpol(const int *n, const double *t, const double *y, double *dydt);

#ifdef __cplusplus
}
#endif

In agreement with the Fortran/C procedure interface conventions, all arguments are passed as pointers. Dummy variables with the intent(in) attribute are also given the const specifier to make sure that calling code in C knows these pointers are not modified internally. The extern "C" block makes sure the procedure f_vanderpolhas C linkage (compiler does not mangle the name). It is similar to the bind(C) attribute in the Fortran code. By using the preprocessor directives the extern "C" block is only added when the header is included by a C++ compiler. This way we only need to maintain a single header file for calling codes in C or C++.

A simple way to generate the C prototypes defined in a Fortran module is to export them automatically. With gfortran you can use the command:

gfortran -fc-prototypes -fsyntax-only vanderpol_mod.f90 > vanderpol_mod.h

Wrapping the Fortran procedure in a C++ MEX function

MATLAB projects typically follow the convention of defining functions in individual files. E.g. the van der Pol equations could have been placed directly in the file vanderpol.m:

function dydt = vanderpol(t,y)
dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];

Similarly, for each C++ (or Fortran) function we want to expose to MATLAB we will have a file defining a C++ class that overrides the function call operator, operator() to create a function object (also called functor). To define a C++ MEX function we first need to import the necessary MATLAB headers:

// vanderpol.cpp

#include "mex.hpp"
#include "mexAdapter.hpp"

After the MATLAB header we wclude the header file we just created containing the C prototype of the Fortran procedure. In the implementation of the C++ MEX function we will also need some definitions from C++ standard library. By convention we place these before the MATLAB headers:

// vanderpol.cpp

/* C++ standard library headers */
#include <vector>

/* MATLAB related headers */
#include "mex.hpp"
#include "mexAdapter.hpp"

/* Header containing prototypes of Fortran routine */
#include "vanderpol_mod.h"

The MEX function will generally have the following structure

class MexFunction : public matlab::mex::Function {
    matlab::data::ArrayFactory factory;
public:
    void operator()(matlab::mex::ArgumentList outputs, matlab::mex::ArgumentList inputs) {

        // check input arguments
        ...
        // implement function
        ...
        // assign outputs
        ...
    }
}

The class inherits from the base matlab::mex::Function class. The name of the class, MexFunction, should not be modified.

C++ MEX Function

Before we create the Fortran glue code, we will look at the simpler case of a pure C++ MEX function. In this case the vanderpol_mod.h header we just created previously is not needed.

In accordance with the odefun argument expected by the ode23s solver, the inputs to our MEX function will be the scalar t and a column vector y. The function should create a new column vector with the RHS values of the ODE and assign it to the list of outputs.

For simplicity we will skip checking the input arguments. Instead we proceed to extract the values from the input list. For the two input arguments this is done by the following code:

        const double t_scalar = inputs[0][0];
        const matlab::data::TypedArray<double> y_vector = inputs[1];

To retrieve the input arguments from the matlab::mex::ArgumentList container we used the operator[]. Each element of the input container is an array of parent class matlab::data::Array. Since we know the first input argument is a scalar, we simply access the first element, hence input[0][0] (don’t forget C and C++ use zero-based indexing). For the second argument containing the vector of dependent variables we create a shared data copy by assigning inputs[1] to an instance of type matlab::data::TypedArray<T> (a sub-class of matlab::data::Array).

Next we proceed to create an empty result array. The dimensions of this array are equal to the dimensions of the input y_vector:

        auto dims = y_vector.getDimensions();
        auto res = factory.createArray<double>(dims);

Arrays are created using the matlab::data::ArrayFactory instance, which is a private member of our MEX function class. The result res is an instance of TypedArray<double>.

To implement the van der Pol equations in C++ we use the operator[] syntax:

        res[0] = y_vector[1];
        res[1] = 1000*(1 - y_vector[0]*y_vector[0])*y_vector[1] - y_vector[0];

Finally, we assign the result array to the list of outputs:

        outputs[0] = res;

The complete MEX function in C++ can be found by opening the box below:

Complete example in C++
/* vanderpol
 * Evaluate the van der POL ODEs for mu = 1
*/

#include "mex.hpp"
#include "mexAdapter.hpp"

using matlab::mex::ArgumentList;
using namespace matlab::data;

class MexFunction : public matlab::mex::Function {
    ArrayFactory factory;
public:
    void operator()(ArgumentList outputs, ArgumentList inputs) {

        const double t_scalar = inputs[0][0];
        const TypedArray<double> y_vector = inputs[1];

        auto dims = y_vector.getDimensions();
        auto res = factory.createArray<double>(dims);

        res[0] = y_vector[1];
        res[1] = 1000*(1 - y_vector[0]*y_vector[0])*y_vector[1] - y_vector[0];

        outputs[0] = res;
    }
};

To compile this function in MATLAB, use the command:

mex vanderpol.cpp

To test the function run the following script:

[t,y] = ode23s(@vanderpol,[0 3000],[2 0]);
plot(t,y(:,1),'-o')

C++ MEX function calling Fortran

Instead of performing the calculations in C++, we want to replace them with a call to our Fortran procedure. Before we can achieve this, we have to coerce the MATLAB data into types matching the C prototype of our Fortran procedure.

The number of equations is easiest. We just access the first element of the dims object:

        auto dims = y_vector.getDimensions();
        int n = static_cast<int>(dims[0]);

The automatic return value of the TypedArray<T> member function getDimensions() is an object of type matlab::data::ArrayDimensions, which is just an alias for std::vector<size_t>. We add a static_cast just to show awareness of the potential overflow that could occur when assigning a variable of type size_t to a smaller int variable. Unless you are planning to solve more than 2**31 - 1 (upper limit of int) ODE’s this is unlikely to ever become a problem.

For the vector y the Fortran procedure expects a pointer to a buffer of double precision values, i.e a double *. The code needed to do this is:

        buffer_ptr_t<double> y_p = y_vector.release();
        double *y = y_p.get();

The release member function returns a smart pointer object (buffer_ptr_t<T> is an alias of std::unique_ptr<T[]> with a custom deleter routine). Since y_vector is a (non-const) shared array, the release method actually creates a copy of the underlying array buffer. (This turns out to be a disadvantages of using the C++ MEX function route to wrap a Fortran procedure. It is a result of incompatible mechanisms programming languages use to provide type-safety and automatic memory management.) The subsequent call to get() returns the stored pointer. For brevity we could also chain the two commands as follows

        double *y = y_vector.release().get()

Now that the input arguments are ready, we procced to create the double precision buffer of output values dydt:

        size_t nelems = n;
        auto dydt_p = factory.createBuffer<double>(nelems);
        double *dydt = dydt_p.get();

Instead of directly creating an array like in the pure C++ version, we now use the MATLAB array factory to create a buffer first. The buffer must be sized with the number of elements. For 2-D arrays or higher, this will be equal to the product of the array dimensions.

Finally we have reached a point where we can call our Fortran routine (watch out for arguments passed by reference):

        f_vanderpol(&n, &t_scalar, y, dydt);

Before we can assign dydt to the MEX function output, we need to coerce it back into a MATLAB array type. Luckily, the factory class offers a method exactly for this purpose:

            TypedArray<double> res = factory.createArrayFromBuffer(
                dims, std::move(dydt_p));

            outputs[0] = res;

When creating the result array, we need to use std::move to force move semantics (whatever that means ). If you forget to use this helper function the MEX function will fail to compile.

The see the complete C++ MEX function calling Fortran snippet open the following link:

Complete example calling Fortran
/* vanderpol
 * Evaluate the van der POL ODEs for mu = 1
*/

#include <iostream>
#include <vector>

#include "mex.hpp"
#include "mexAdapter.hpp"

#include "vanderpol_mod.h"


using matlab::mex::ArgumentList;
using namespace matlab::data;

class MexFunction : public matlab::mex::Function {
    ArrayFactory factory;
public:
    void operator()(ArgumentList outputs, ArgumentList inputs) {

        // Retrieve inputs
        const double t_scalar = inputs[0][0];
        TypedArray<double> y_vector = inputs[1];

        // Get size of column vector
        auto dims = y_vector.getDimensions();
        int n = static_cast<int>(dims[0]);

        // Prepare pointer to y values
        buffer_ptr_t<double> y_p = y_vector.release();
        double *y = y_p.get();

        // Prepare pointer to dydt values
        size_t nelems = n;
        auto dydt_p = factory.createBuffer<double>(nelems);
        double *dydt = dydt_p.get();

        // Call Fortran procedure
        f_vanderpol(&n,&t_scalar,y,dydt);

        // Create output array
        TypedArray<double> res = factory.createArrayFromBuffer(
            dims, std::move(dydt_p));

        // Assign output
        outputs[0] = res;
    }
};

Compilation and linking

With the procedures ready what remains to be done is compiling the executable code. First we navigate to our work directory and compile the the Fortran code:

$ gfortran -Wall -O2 -fPIC -c vanderpol_mod.f90

This creates the object file vanderpol_mod.o. Next, from the MATLAB command window, run the command

mex -v vanderpol.cpp vanderpol_mod.o

which links the Fortran code statically. For distribution and reuse it is probably better to place the Fortran procedures in a shared library, but I did not explore how to achieve this.

Again we can run the script

[t,y] = ode23s(@vanderpol,[0 3000],[2 0]);
plot(t,y(:,1),'-o')

showing the expected output

vanderpol

7 Likes