Full matrix or triagonal part necessary in LAPACKE function for diagonalization?

Question

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!

What is the output of the code? – stephematician Oct 01 '19 at 10:09 — stephematician, Oct 01 '19 at 10:09

francis · Accepted Answer · 2019-10-02T18:34:47.430

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,...).

score 1 · Answer 2 · answered Oct 01 '19 at 10:47

According to this example, it seems that LAPACKE_dsyev requires an upper or lower triangular part to be stored in the full storage scheme. Even documentation says, that parameter a should be an array of size max(1, lda * n), which is 16 in your case.

When I switched to this definition of the matrix:

double matrix_ex_uppertri[16] =
   { -2.0,  0.0,  0.5,  0.0,
      0.0,  0.5, -2.0,  0.5,
      0.0,  0.0,  0.5,  0.0,
      0.0,  0.0,  0.0, -1.0 };

then, I got the the correct eigenvalues.

Need to admit that the documentation is very misleading here. Don't understand why they do not specify the storage scheme clearly.

Ok I see, thank you very much for the clarification, that does make sense! — cschoenle, Oct 01 '19 at 11:01

Full matrix or triagonal part necessary in LAPACKE function for diagonalization?

2 Answers2