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

Question

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

LinearAlgebra`BLAS`

or

LinearAlgebra`LAPACK`

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:

LinearAlgebra`LAPACK`GETRS[];
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.

Answer 1

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";
Quiet[LibraryFunctionUnload[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-&gt;MTensor_getRealData(A_in);
mreal * lambda_real  = libData-&gt;MTensor_getRealData(lambda_real_out);
mreal * lambda_imag  = libData-&gt;MTensor_getRealData(lambda_imag_out);
mreal * V            = libData-&gt;MTensor_getRealData(V_out);

const int n  = static_cast&lt;int&gt;( libData-&gt;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_( \&quot;N\&quot;, \&quot;V\&quot;, &amp;n, A, &amp;n, lambda_real, lambda_imag, nullptr, &amp;n, V, &amp;n, &amp;buffer_dummy, &amp;buffer_size, &amp;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_( \&quot;N\&quot;, \&quot;V\&quot;, &amp;n, A, &amp;n, lambda_real, lambda_imag, nullptr, &amp;n, V, &amp;n, buffer, &amp;buffer_size, &amp;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-&gt;MTensor_disown(lambda_real_out); 
libData-&gt;MTensor_disown(lambda_imag_out);
libData-&gt;MTensor_disown(V_out);

return info;

}"
   ,
   "Text"
   ];
lib = CreateLibrary[{file}, name,
   "ShellOutputFunction" -> Print,
   "CompileOptions" -> {}
   , "LinkerOptions" -> {"-llapack"}
   , "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.)
     }
   ];
cf = LibraryFunctionLoad[lib, name,
  {
   {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.