Как я могу получить правильные собственные векторы из LAPACK ssyev? - PullRequest
2 голосов
/ 05 июля 2019

Я пытаюсь использовать 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
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...