Calling Fortran routines from C/C++ code
From my grad school days, I had heard that it was a common practice to mix C/C++ code with Fortran libraries such as LAPACK
or BLAS
. So I set about trying to leverage LAPACK
’s eigenvalue solver routines for blitzdg
’s EigenSolver class.
The process wasn’t overly challenging, but there were some intricacies that needed special attention. I’ll spell those out here in hopes that it may help someone in the future, or that it will serve as a reference for myself.
Declaring an external Fortran routine as a C/C++ function
In my case, I want to call the Fortran routine called dsyevd
, which is a workhorse for solving eigenvalue problems of real, symmetric matrices. Its official API doc is here.
In order to call dsyevd
in my C++ code, I need to declare a function with the same signature as the Fortran routine, but it’s a bit more convoluted than that. Here is what the declaration code ends up looking like:
Questions you might have:

Why is it wrapped in an
extern "C" {}
block? This tells the compiler to give the function Cstyle linkage. C++ compilers will often mangle the name. This is bad, because we need the function name to match up with the Fortran routine, otherwise the linker will not associate the two pieces of code.

Why is there an underscore after the function’s name?
 Fortran compilers like
gfortran
will usually append an underscore to the routine’s name in the object file. There is a good discussion on Fortran/C++ interoperability here. Thus, in order for the linker to find the code associated with our function, we append an underscore to the name.
 Fortran compilers like

Why are all the arguments pointers?
 This is a difference between how C/C++ and Fortran compilers work. C/C++ passes arguments by value. Fortran always passes arguments by reference. Therefore, the C++ declaration of the Fortran routine contains a pointer in every argument.
Calling a LAPACK routine from C/C++ code
The next step is to actually call the C function we declared which links to a Fortran routine in a C++ method. The code looks like this:
First, we set some parameters and convert the blitz matrix A
to a contiguous memory block (a 1D array) in MatrixConverter.fullToPodArray
.
The dsyevd_
function then gets called twice. First to determine the optimal workspace parameters, then to actually solve the eigenvalue problem. The Fortran routine stores the resulting eigenvectors in the storage allocated by the array Apod
(Here, pod
= plainold datatype).
The eigenvectors are then converted from pod
format to 2D blitzarray format, and similarly the eigenvalues are put into the 1D blitz array that the caller passes in.