Did anybody come across any resources of interfacing Modern Fortran with Matlab. The documentation provides very outdated API with F77. Have you used the c interop route ?
The resources are scarce on this topic. But if you use CFI, then any doc that applies to C, also applies to Fortran. It took me a year to implement a Fortran-MATLAB interface that I needed. But we are here to not let such delays happen anymore. If you elaborate further on your MFI case, I might be able to help. In general, you need to write a MEX file, which is typically a C code that calls MATLAB-specific C libraries to construct the interface bridge. Inside the Mex file, you so call your Fortran dll file. Then you compile the mex file and call the MEX library (which is really just another dll file but with a different name) and the Mex file calls your Fortran dll.
1 Like
The Fortran Wiki has pages on interoperability, including a page for Matlab, where Fortran 95 Interface to MATLAB API with extras by James Tursa seems relevant.
Alternatively one could write Fortran 77 style wrappers for modern Fortran code and have Matlab call the F77 code.
1 Like
Thanks @shahmoradi.
I have a Fortran code for solving stiff differential equations. Matlab also has it - ode23s, but it is orders of magnitude slower than Fortran implementation.
Just to outline it’s working -
I have 3 stiff set of odes and are to be integrated together wrt time (independent variable).
My solver takes time steps t and t+1 and 3 dependent variables from the previous step (t) as input and solves for 3 dependent variables at t+1 as output.
To be more clear on C-interop, do I just need to define the interface such that on my Fortran side I expose the inputs and outputs to the c-mex file and rest of the workings (variables and subroutines and function) can be in Fortran ? Or should I change some internal workings also ? In other words, should my c-mex file know anything other than inputs and outputs ?
1 Like
So my understanding is that you need to do a callback to MATLAB from Fortran, right? That is, you want MATLAB to load a Fortran shared file that contains the ode solver, then pass a MATLAB callback function odefun to the Fortran library so that Fortran calls your MATLAB function to return the solution. This is the interface of ode23s in MATLAB that I see on the MATLAB website,
[t,y] = ode23s(odefun,tspan,y0)
If this is the problem, then I believe it is possible (though not quite straightforward, all because of MATLAB’s weird interfacing conventions, not because of Fortran or C). But before doing so, are you sure the performance bottleneck is in the solver and not in the odefun? It would be a surprise if MATLAB intrinsics are significantly slower than C/Fortran.
1 Like
Yes ! I think the bottleneck is coming from interp1 function that I am using in odefun. Thanks @shahmoradi. I am not looking for callback to MATLAB from Fortran. My odefun and solver both are in Fortran, I am using this ODE solving as a part of a bigger problem - where I need output from the solver to be included as a part of calculation that I am doing in MATLAB.
Also I am trying to write a robust API for calling Fortran so that I can use the beauty of MATLAB (plotting/GUI/Data-analysis) and speed of Fortran.
1 Like
In such a case, it may actually be more sensible to write the results to an output file and read it in MATLAB for further post-processing. External file IO is typically on the order of 100-1000 microseconds. If you need faster information exchange between Fortran and MATLAB, on the order of 1-10 microseconds, then writing a MEX file would be sensible.
1 Like
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:
- Writing the Fortran procedure
- Defining the C prototype of the Fortran procedure
- Wrapping the Fortran procedure in a C++ MEX function
- 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
7 Likes
Just for comparison, here is a C MEX file for the same callback:
// vanderpol.c
#include "mex.h"
#include "matrix.h"
#include "vanderpol_mod.h"
void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
// Value of independent variable
const double *t = mxGetDoubles(prhs[0]);
// Get number of rows in column vector y
const int m = mxGetM(prhs[1]);
// Get pointer to dependent variables
const double *y = mxGetDoubles(prhs[1]);
// Create output buffer
plhs[0] = mxCreateDoubleMatrix(m,1,mxREAL);
double *dydt = mxGetDoubles(plhs[0]);
f_vanderpol(&m,t,y,dydt);
}
To compile this function use the command
mex -R2018a vanderpol.c vanderpol_mod.o
3 Likes
Thanks a lot @ivanpribec !!! I am exactly looking for this.
Great help !
The c++ side looks little daunting to me, but I will learn it soon.
The C++ API is newer and offers the benefits like better type-safety and automatic memory management (due to use of C++ containers). It is also part of a general trend that more and more software is C++ (as are the developers).
For interfacing to Fortran I was surprised by how much easier it was to write the C MEX function. You just need to follow the steps
- unpack contents of
prhsinto corresponding C types - prepare any output buffers to pass to Fortran
- call your Fortran routine
If you avoid doing any other operations in the C MEX function you should remain pretty safe.
3 Likes
If you need to reduce execution speed, fortran is hard to beat. However, if you just want the functionality of the fortran code available in matlab, then refactor the goto’s, convert to f90 style, and try f2matlab. F2matlab converts fortran source to matlab source and works on smaller codes relatively well, but I can help with the process for larger codes.
Another option is to use system calls from matlab to call fortran executables, then read in the results in matlab.
I have gone the MEX file route before, but it can be painful and tedious. But if you need to retain the execution speed, then stick with fortran.
3 Likes
@Ashok interp1 is indeed very slow. I would suggest you to use interp1qr from Matlab file exchange. Here is the link. Alternatively, I’ve written myself a fast function to do linear interpolation: you can find it in my Github page: Income-Fluctuation-Problem/myinterp1.m at main · aledinola/Income-Fluctuation-Problem · GitHub
On calling Fortran from Matlab: unfortunately the documentation available on Mathworks is only for FORTRAN 77 but if you write the subroutine in modern fortran it works as well. I’ve done it several times: first, you find the bottleneck in your Matlab code. Then, you outsource the slow part to a Fortran subroutine that you call from Matlab using a gateway function. If I have time I’ll post an example on my github.
Thanks a lot @aledinola. Yeah interp1 is slow. I ended up writing everything in Fortran - and also one variant with what Ivan has pointed with C mex file.
I will take a look at your interp1 also.
Thanks.
1 Like
I find it easier to just use Fortran throughout instead of C/C++ in the Matlab gateway function. Yes, you have to use F77 formatting, but most modern features I’ve tried works anyway 
1 Like