Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
7113 Discussions

Diagonalization of symmetric matrices - form of matrix?

Schoenle__Christoph
1,598 Views

I would like to compute the eigenvalues of a symmetric matrix and wanted to use the LAPACKE_dsyev function from the MKL Library in C++ for that.

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 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;
	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];
// Pass over the full matrix
	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 << endl;
// Pass only the upper triagonal
	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 << 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

I guess I am missing something obvious here, but why, if the full matrix has to be given anyway, is it necessary at all to specify 'U' or 'L'? Or am I doing something wrong elsewhere?

Thanks for any help, and apologies if this question is trivial (which I have the feeling it must be), I am rather new to using MKL.

0 Kudos
2 Replies
Schoenle__Christoph
1,598 Views

Ok, I have been given an answer already on Stackoverflow (https://stackoverflow.com/questions/58181635/full-matrix-or-triagonal-part-necessary-in-lapacke-function-for-diagonalization).

Apparently, an array of the full size has to be given, and the flag 'U' or 'L' then determines which triagonal of the matrix is used and which is ignored (so the latter can contain anything like zeros or random data).

So would conclude that everything does make sense here, documentation is just slightly misleading.

0 Kudos
mecej4
Honored Contributor III
1,598 Views

You only need to fill in the lower triangle with values. The strictly upper triangle need not be initialized with values. If eigenvectors are also requested, the full matrix is used to hold the eigenvectors. The way the subroutine is designed, it has two parts which have to be mutually consistent as to the representation of the matrix. If you don't ask for eigenvectors, the second part is skipped, but the subroutine arguments are still the same, because the underlying code was written Fortran 77, which does not permit arguments to be skipped whether they get used or not.

0 Kudos
Reply