# MATLAB - Fortran interface

**URL:** <https://fortran-lang.discourse.group/t/matlab-fortran-interface/1322>\
**Category:** Help\
**Created:** [May 29, 2021, 4:15pm UTC](https://fortran-lang.discourse.group/t/matlab-fortran-interface/1322 "2021-05-29T16:15:04Z")\
**Posts on this page:** 1\
**Showing post:** 9

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [May 31, 2021, 12:59am UTC](https://fortran-lang.discourse.group/t/matlab-fortran-interface/1322/9 "2021-05-31T00:59:37Z")

</div>

## Using a Fortran callback in the MATLAB ODE solvers

MATLAB provides multiple [external language interfaces](https://uk.mathworks.com/help/matlab/external-language-interfaces.html?s_tid=CRUX_lftnav). 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](https://uk.mathworks.com/help/matlab/cc-mx-matrix-library.html).

The [C++ MEX API](https://uk.mathworks.com/help/matlab/cpp-mex-file-applications.html) 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`](https://uk.mathworks.com/help/matlab/apiref/matlab.data.array.html) class defined in the [MATLAB Data API](https://uk.mathworks.com/help/matlab/matlab-data-array.html). 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](https://uk.mathworks.com/help/matlab/ref/ode23s.html#bu654ce) from MATLAB’s [`ode23s`](https://uk.mathworks.com/help/matlab/ref/ode23s.html) solver documentation. The van der Pol ODE’s in Fortran are defined by the procedure `f_vanderpol`:

```plaintext
! 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:

```cpp
// 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_vanderpol`has 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](https://gcc.gnu.org/onlinedocs/gfortran/Interoperability-Options.html). With `gfortran` you can use the command:

```auto
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`:

```auto
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:

```cpp
// 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](https://google.github.io/styleguide/cppguide.html#Names_and_Order_of_Includes) we place these before the MATLAB headers:

```cpp
// 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

```cpp
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`](https://uk.mathworks.com/help/matlab/matlab_external/structure-of-c-mex-function.html) 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`](https://uk.mathworks.com/help/matlab/ref/ode23s.html#bu5ypsh_sep_shared-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:

```cpp
        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`](https://uk.mathworks.com/help/matlab/matlab_external/cpp-mex-api.html#mw_bfc7456f-1d43-48d4-9225-f0b6e9060beb) container we used the [`operator[]`](https://uk.mathworks.com/help/matlab/matlab_external/cpp-mex-api.html#mw_bfc7456f-1d43-48d4-9225-f0b6e9060beb). Each element of the input container is an array of parent class [`matlab::data::Array`](https://uk.mathworks.com/help/matlab/apiref/matlab.data.array.html). 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](https://uk.mathworks.com/help/matlab/matlab_external/copy-cpp-api-matlab-arrays.html) by assigning `inputs[1]` to an instance of type [`matlab::data::TypedArray<T>`](https://uk.mathworks.com/help/matlab/apiref/matlab.data.typedarray.html) (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`:

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

```

Arrays are created using the [`matlab::data::ArrayFactory`](https://uk.mathworks.com/help/matlab/apiref/matlab.data.arrayfactory.html) 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:

```cpp
        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:

```cpp
        outputs[0] = res;

```

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

> **Complete example in C++**
>
> ```auto
> /* 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:

```matlab
mex vanderpol.cpp

```

To test the function run the following script:

```matlab
[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:

```cpp
        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:

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

```

The [`release`](https://uk.mathworks.com/help/matlab/apiref/matlab.data.typedarray.html#bvmn12r-2) member function returns a smart pointer object (`buffer_ptr_t<T>` is an alias of [`std::unique_ptr<T[]>`](http://www.cplusplus.com/reference/memory/unique_ptr/) 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()`](http://www.cplusplus.com/reference/memory/unique_ptr/get/) returns the stored pointer. For brevity we could also chain the two commands as follows

```cpp
        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`:

```cpp
        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):

```cpp
        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:

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

            outputs[0] = res;

```

When creating the result array, we need to use [`std::move`](http://www.cplusplus.com/reference/utility/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**
>
> ```auto
> /* 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:

```plaintext
$ 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

```matlab
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

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

```

showing the expected output

 ![vanderpol](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/db05851a8981838ef3c64bbf6144c32d999bbdf7.png)

---

_[View the full topic](https://fortran-lang.discourse.group/t/matlab-fortran-interface/1322)._