3

I am trying to use LAPACK's ssyev (http://www.netlib.org/lapack/explore-html/d3/d88/group__real_s_yeigen_ga63d8d12aef8f2711d711d9e6bd833e46.html) to calculate eigenvalues and eigenvectors. As a test, I'm using the simple matrix from https://math.stackexchange.com/questions/1271788/finding-eigenvectors-of-symmetric-matrices-when-the-general-solutions-of-eigenva:

2 1 1
1 2 1
1 1 2

I get the same eigenvalues (4, 1, 1) back as in that question, however, my eigenvectors are different. The expected values are

1  1  1
1  0 -2
1 -1  1

but I get

0.57735 1           1
0.57735 0.408248    1
0.57735 0.408248    0.707107

How can I get the correct eigenvectors using LAPACK ssyev?

The code:

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <random>
#include <chrono>
#include <mkl.h>


typedef float T;

void (*gemm)(const CBLAS_LAYOUT, const CBLAS_TRANSPOSE, const CBLAS_TRANSPOSE, const long long, const long long, const long long, const float, const float *, const long long, const float *, const long long, const float, float *, const long long) = &cblas_sgemm;
long long (*syev)(int, char, char, long long, float *, long long, float *) = &LAPACKE_ssyev;

int (*ssNewTask)(VSLSSTaskPtr *, const long long *, const long long *, const long long *, const float *, const float *, const long long *) = &vslsSSNewTask;
int (*ssEditTask)(VSLSSTaskPtr, const long long, const float *) = &vslsSSEditTask;
int (*ssEditMoments)(VSLSSTaskPtr, float *, float *, float *,  float *,  float *, float *,  float *) = &vslsSSEditMoments;
int (*ssEditCovCor)(VSLSSTaskPtr,  float *,  float *, const long long *,  float *, const long long *) = &vslsSSEditCovCor;
int (*ssCompute)(VSLSSTaskPtr, const unsigned long long, const long long) = &vslsSSCompute;

#ifdef SS_FAST

int ss_method = VSL_SS_METHOD_FAST;

#else

int ss_method = VSL_SS_METHOD_1PASS;

#endif

int main(int argc, char* argv[]) {
    MKL_INT nRows = 3;

    T* MATRIX = new T[nRows * nRows]; 
    // Input matrix:
    //  2 1 1
    //  1 2 1
    //  1 1 2
    MATRIX[0] = 2;
    MATRIX[1] = 1;
    MATRIX[2] = 1;
    MATRIX[3] = 1;
    MATRIX[4] = 2;
    MATRIX[5] = 1;
    MATRIX[6] = 1;
    MATRIX[7] = 1;
    MATRIX[8] = 2;

    T* eigenvalues = new T[nRows];

    syev(LAPACK_ROW_MAJOR, 'V', 'U', nRows, MATRIX, nRows, eigenvalues);
    // eigenvalues and vectors are outputted in reverse. This gives:
    // eigenvalue 1: 4
    // eigenvalue 2: 1
    // eigenvalue 3: 1
    std::cout << "eigenvalue 1: " << eigenvalues[2] << std::endl;
    std::cout << "eigenvalue 2: " << eigenvalues[1] << std::endl;
    std::cout << "eigenvalue 3: " << eigenvalues[0] << std::endl<< std::endl;


    // Expected output matrix:
    // 1  1  1
    // 1  0 -2
    // 1 -1  1
    //
    // But get
    // 0.57735  1           1
    // 0.57735  0.408248    1
    // 0.57735  0.408248    0.707107
    std::cout << MATRIX[8] << "\t" << MATRIX[7] << "\t" << MATRIX[6] << std::endl;
    std::cout << MATRIX[5] << "\t" << MATRIX[4] << "\t" << MATRIX[3] << std::endl;
    std::cout << MATRIX[2] << "\t" << MATRIX[1] << "\t" << MATRIX[0] << std::endl;
    std::cout<<std::endl;
}

and compiled it with

clang++ -Wall -std=c++11 -stdlib=libc++ -DMKL_ILP64 -m64 -DSINGLE_PRECISION -I/path/to/mkl/compilers_and_libraries_2019.3.199/linux/mkl/include -c pca.debug.cc -o pca.debug.o
clang++ -fopenmp=libiomp5 -DMKL_ILP64 -Wall -std=c++11 -stdlib=libc++ -DMKL_ILP64 -m64 -DSINGLE_PRECISION -L/path/to/mkl/compilers_and_libraries_2019.3.199/linux/mkl/lib/ \
-I/path/to/mkl/compilers_and_libraries_2019.3.199/linux/mkl/include/ -lmkl_intel_ilp64 -lmkl_core -lmkl_intel_thread -ldl -lpthread -lm -o pcadebug pca.debug.o
Niek de Klein
  • 8,524
  • 20
  • 72
  • 143
  • Two eigenvalues are identical (1). This means that any linear combination of two eigenvectors (corresponding to this eigenvalue) will be a valid eigenvector! In other words, there is no unicity – Damien Jul 05 '19 at 14:15
  • That being said, effectively, the vectors you get are not correct – Damien Jul 05 '19 at 14:22
  • Related: https://stackoverflow.com/questions/32327760/how-to-use-dsyev-routine-to-calculate-eigenvalues – Damien Jul 05 '19 at 14:41

0 Answers0