Я пытаюсь использовать sseev LAPACK (http://www.netlib.org/lapack/explore-html/d3/d88/group__real_s_yeigen_ga63d8d12aef8f2711d711d9e6bd833e46.html) для вычисления собственных значений и собственных векторов. В качестве теста я использую простую матрицу из 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
Я получаю те же собственные значения (4, 1, 1) обратно, что и в этом вопросе, однако мои собственные векторы разные. Ожидаемые значения
1 1 1
1 0 -2
1 -1 1
но я получаю
0.57735 1 1
0.57735 0.408248 1
0.57735 0.408248 0.707107
Как я могу получить правильные собственные векторы, используя LAPACK ssyev?
код:
#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;
}
и скомпилировал его с
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