I would like to compute the eigenvalues of a symmetric matrix and wanted to use the LAPACKE_dsyev function from the Intel MKL Library in C++ for that. But I am a bit puzzled regarding the form in which the matrix needs to be passed.
From the documentation https://software.intel.com/en-us/mkl-developer-reference-c-syev, I concluded that I would have to pass only the upper/lower triangular part of the matrix. It says there about the argument that it "is an array containing either upper or lower triangular part of the symmetric matrix A". However, it seems that actually one needs to pass the pointer to the full matrix to the routine.
Say I want to diagonalize the following matrix:
[[-2, 0, 0.5, 0],
[0, 0.5, -2, 0.5],
[0.5, -2, 0.5, 0],
[0, 0.5, 0, -1]],
which has eigenvalues [ 2.56, -2.22, -1.53, -0.81]
Then in the following code, only the first option gives the correct values.
#include <iostream>
#include"mkl_lapacke.h"
using namespace std;
int main(){
MKL_INT N = 4;
// use full matrix
double matrix_ex_full[16] = {-2,0,0.5,0,0,0.5,-2,0.5,0.5, -2, 0.5, 0, 0,0.5,0,-1};
double evals_full[4];
MKL_INT test1 = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'N', 'U', N, matrix_ex_full,N, evals_full);
cout << "success = " <<test1 << endl;
for (MKL_INT i = 0;i<4;i++)
cout << evals_full[i] << endl;
// use upper triagonal only
double matrix_ex_uppertri[10] = {-2, 0, 0.5, 0, 0.5, -2, 0.5, 0.5, 0, -1};
double evals_uppertri[4];
MKL_INT test2 = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'N', 'U', N, matrix_ex_uppertri,N, evals_uppertri);
cout << "success = " <<test2 << endl;
for (MKL_INT i = 0;i<4;i++)
cout << evals_uppertri[i] << endl;
// (Compiled with g++ test.cpp -o main -m64 -I/share/opt/intel/compilers_and_libraries_2018.0.128/linux/mkl/include -L/share/opt/intel/compilers_and_libraries_2018.0.128/linux/mkl/lib/intel64 -Wl,--no-as-needed -lmkl_intel_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl)
}
Why do I only get the correct eigenvalues when passing over the pointer to the full matrix?
I have the feeling this must be a trivial question, but what am I missing? If the full matrix always has to be specified, then it makes no sense to me to specify with 'U' or 'L' which triangle of the matrix is been given. Or am I doing something wrong elsewhere?
Thanks a lot for your help!
According to the documentation of Lapack regarding naming scheme, letters sy in the name of function dsyev() refer to a symmetric matrix, while letters d refers to double precision and ev refers to eigenvalues. Nevertheless, The format sy corresponds to a conventional storage as a 2D array of shape consistent with that of the matrix. Either the upper triangular part or the lower part is used depending on the value of argument UPLO.
To use a packed format, take a look at the function dspev(), where the letters sp corresponds to packed storage of symmetric matrices. The function LAPACKE_dspev() from LAPACKE provides a convenient interface to C.
Here is a sample code compiled by g++ main.cpp -o main -llapacke -llapack -lblas -lm -Wall:
#include <iostream>
#include <math.h>
extern "C" {
#include <lapacke.h>
}
using namespace std;
int main(int argc, char *argv[])
{
lapack_int N = 4;
// use full matrix
double matrix_ex_full[16] = {-2,0,0.5,0, 0,0.5,-2,0.7, 0.5, -2, 0.5, 0, 0,0.7,0,-1};
double evals_full[4];
lapack_int test1 = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'N', 'U', N, matrix_ex_full,N, evals_full);
cout << "success = " <<test1 << endl;
for (int i = 0;i<4;i++)
cout << evals_full[i] << endl;
// use upper triagonal only
double matrix_ex_uppertri[10] = {-2, 0, 0.5, 0, 0.5, -2, 0.7, 0.5, 0, -1};
double evals_uppertri[4];
int test2 = LAPACKE_dspev(LAPACK_COL_MAJOR, 'N', 'L', N, matrix_ex_uppertri, evals_uppertri,NULL,N);
cout << "success = " <<test2 << endl;
for (int i = 0;i<4;i++)
cout << evals_uppertri[i] << endl;
return 0;
}
As shown in the docs of LAPACKE_dspev_work(), making use of LAPACK_COL_MAJOR saves some extra allocatation of temporary arrays. Since the eigenvectors are not computed, the outcome remains consistent with that of LAPACKE_dsyev(LAPACK_ROW_MAJOR,...).