# Is it possible to use the DGEEV and DSYEV LAPACK subroutines in Mathematica?

Here is my problem: I'm diagonalizing some matrices using Eigensystem[] to obtain eigenvalues and eigenvectors and I'm diagonalizing the same matrices using a Fortran code that uses DGEEV (DGEEV computes the eigenvalues and right eigenvectors for a real nonsymmetric matrix) or DSYEV (DSYEV computes eigenvalues and eigenvectors of a real symmetric matrix) LAPACK subroutines.

These matrices are of course square but can be real nonsymmetric or symmetric matrices. The problem is that when I compare what I obtain using Mathematica and Fortran I have the same eigenvalues (that are degenerate) but not the same eigenvectors.

So I know that for degenerate eigenvalues Mathematica can give non orthogonal eigenvectors (found it here) and I tried stuffs with Orthogonalize[] and Normalize[] but nothing conclusive yet (I'm still working in this direction).

So I was wondering if any of you know if we can use the DGEEV or DSYEV LAPACK subroutines in Mathematica ?

I have done some research and found some stuffs.

I found the Mathematica equivalent functions to the LAPACK functions as stated by Nasser with a table that we can find here.

I also found here that some (or all ?) BLAS and some LAPACK functions can be used directly in Mathematica by using

LinearAlgebraBLAS


or

LinearAlgebraLAPACK


followed by the name and the right arguments of the function.

I have tried some of these functions without specifying the arguments and Mathematica let me know that I don't have the right number of arguments, here an example:

LinearAlgebraLAPACKGETRS[];

LinearAlgebraLAPACKGETRS::argrx: LinearAlgebraLAPACKGETRS called with 0 arguments; 4 arguments are expected.


So I tried the same with the DGEEV function (and some other LAPACK functions like the DSYEV subroutine) and Mathematica said nothing.

I also found here something with SystemModel[] by using:

SystemModel["Modelica.Math.Matrices.LAPACK.dgeev", "ModelicaDisplay"]


But then I don't know if it's possible to use it and if yes, how to use it.

• This answer uses GEEV; you might want to study the syntax in that function as a guide on what is retained and omitted compared to the original FORTRAN (e.g. the dimension arguments are typically omitted). SYEV should be similar. (See also this related question.) Jun 4, 2022 at 12:37
• Thanks for your comment but I still don't have the same answer than the DGEEV subroutine. But I will continue to look at your answer. Thanks ! Jun 4, 2022 at 13:03
• You might want to elaborate when you say you "still don't have the same answer"; otherwise we're all just guessing here. Jun 4, 2022 at 13:05
• I mean I still have different eigenvectors with your proposition and the DGEEV or DSYEV subroutines in Fortran. Jun 4, 2022 at 14:03
• "different eigenvectors" - different, how? Did you account for normalization, or was the usual check identity $\mathbf A\mathbf x=\lambda\mathbf x$ not satisfied for some of the eigenvectors generated? Jun 4, 2022 at 15:03

You can use LibraryLink to call LAPACK directly. It involves however quite a lot of boilerplate code. Here is an example:

Needs["CCompilerDriver"];
name = "cf";
ClearAll["cf"]

file = Export[FileNameJoin[{\$TemporaryDirectory, name, ".cpp"}],
"
#include \"WolframLibrary.h\"
#include <cmath>
#include <lapack.h>

EXTERN_C DLLEXPORT int " <> name <>
"(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
// Tell the LibraryFunction what the types of its arguments are.

MTensor A_in            = MArgument_getMTensor(Args[0]);
MTensor lambda_real_out = MArgument_getMTensor(Args[1]);
MTensor lambda_imag_out = MArgument_getMTensor(Args[2]);
MTensor V_out           = MArgument_getMTensor(Args[3]);

// Get the pointers to the buffers of MTensors, so that we can pass them to LAPACK functions.

mreal * A            = libData->MTensor_getRealData(A_in);
mreal * lambda_real  = libData->MTensor_getRealData(lambda_real_out);
mreal * lambda_imag  = libData->MTensor_getRealData(lambda_imag_out);
mreal * V            = libData->MTensor_getRealData(V_out);

const int n  = static_cast<int>( libData->MTensor_getDimensions(A_in)[0] );

// Request the optimal buffer_size for the scratch buffer.

mreal buffer_dummy = 0.;
int   buffer_size  = -1;
int   info = 0;

dgeev_( \"N\", \"V\", &n, A, &n, lambda_real, lambda_imag, nullptr, &n, V, &n, &buffer_dummy, &buffer_size, &info );

if( info == 0 )
{
buffer_size = round(buffer_dummy);

// Allocate some scratch buffer on which LAPACK can work.

mreal * buffer = (mreal*)malloc( buffer_size * sizeof(mreal) );

// The actual call to LAPACK.

dgeev_( \"N\", \"V\", &n, A, &n, lambda_real, lambda_imag, nullptr, &n, V, &n, buffer, &buffer_size, &info );

// We have to clean up after ourselves...
free(buffer);

}

// Make sure that the MTensors passed by reference are _not_ cleaned up by the LibraryFunction.

libData->MTensor_disown(lambda_real_out);
libData->MTensor_disown(lambda_imag_out);
libData->MTensor_disown(V_out);

return info;
}"
,
"Text"
];
lib = CreateLibrary[{file}, name,
"ShellOutputFunction" -> Print,
"CompileOptions" -> {}
, "IncludeDirectories" -> {
"/opt/local/include"    (*I installed LAPACK via macports, so this is the path where to look for the header file.*)
}
, "LibraryDirectories" -> {
"/opt/local/lib"        (*I installed LAPACK via macports, so this is the path where to look for the library file.*)
}
];

{
{Real, 2, "Constant"}, (*argument is passed as constant reference*)
{Real, 1, "Shared"},   (*argument is passed by mutable reference; can be modified.*)
{Real, 1, "Shared"},   (*argument is passed by mutable reference; can be modified.*)
{Real, 2, "Shared"}    (*argument is passed by mutable reference; can be modified.*)
},
"Void" (*No return value.*)
]


And here is a basic usage example:

n = 4;
A = RandomReal[{-1, 1}, {n, n}];
A = A\[Transpose] . A;

(*Allocate arrays for the return values.*)
\[Lambda]real = ConstantArray[0., n];
\[Lambda]imag = ConstantArray[0., n];
U = ConstantArray[0., {n, n}];

(*Hand A and the arrays for the return values over to the LibraryFunction .*)
cf[A, \[Lambda]real, \[Lambda]imag, U];
(*Return values are now written to \[Lambda]real,\[Lambda]imag and U.*)

{\[Mu], V} = Eigensystem[A];
Max[Abs[A . Transpose[V] - Transpose[V] . DiagonalMatrix[\[Mu]]]];

Max[Abs[Sort[\[Lambda]real] - Sort[\[Mu]]]]


5.60663*10^-15

That cf literally returns the output of dgeev. So the outputs might have to be postprocessed, in particular when nonreal eigenvalues are present. See netlib.org for details.

• Thank you for your answer but I'm not used to headers. I will have a look at it. For the moment I have an error with the #include <lapack.h>. I guess that I have to make a directory with all the headers requested ? Thanks again. Jun 4, 2022 at 13:10
• Yes, you to install an appropriate LAPACK distribution ad to make sure that the C++ compiler can find it. Unfortunately, this depends quite heavily on the system that you use. Jun 4, 2022 at 13:13
• Ok I will look at it, thanks. I'm on MacOs BigSur 11.6.2. Jun 4, 2022 at 14:04
• Good, then I can help you. Make sure that XCode is installed and then install Macports. Afterwards, you can execute sudo port install OpenBLAS` in the command line to install OpenBLAS (which ships a LAPACK implementation). Then the above code should run fine. Jun 4, 2022 at 14:08
• Yes the code runs fine ! Thanks a lot ! Jun 4, 2022 at 14:44