LAPACKE eigen解決方案不準確。如何改進？

Question

我需要Ubuntu下的C eigenproblem解析器。爲此我給LAPACKE從LAPACK 3.5.0一個鏡頭和實際上設法寫低示例程序，它是我從正交系統和對角矩陣

EV = [ 
    .6, -.8, 0 
    .8, .6, 0 
    0, 0, 1 
]; 
D = [ 
    2, 0, 0 
    0, -3, 0 
    0, 0, 0 
]; 

通過產生構造 例如：= EV d EV」 。

雖然較低的程序運行良好，但結果奇怪地不準確。 這裏輸出的末尾：

... 
Lambda: -3.07386, 0, 1.87386 
EV = [ 
    -0.788205  0  -0.615412 
    0.615412  0  -0.788205 
    -0    1  0 
]; 
Info: 0 

爲我所用www.physics.orst.edu/~rubin/nacphy/lapack/routines/dsyev.html

文檔我的完整的源：

/** 
* eigen.cpp 
* 
* Given that you put version 3.5.0 into /opt/lapack/ compile this with: 
* g++ eigen.cpp -I"/opt/lapack/lapack-3.5.0/lapacke/include" \ 
    -L"/opt/lapack/lapack-3.5.0" -llapacke -llapack -lblas -lcblas 
* The order of included libraries is important! 
*/ 

#include <iostream> 
#include <string> 
#include <sstream> 
// cstdlib needs including before clapack! 
#include <cstdlib> 
#include <cblas.h> 
#include <lapacke.h> 

using namespace std; 

/** Column major style! */ 
string matrix2string(int m, int n, const double* A) 
{ 
    ostringstream oss; 
    for (int j=0;j<m;j++) 
    { 
    for (int k=0;k<n;k++) 
    { 
     oss << A[j+k*m] << "\t"; 
    } 
    oss << endl; 
    } 
    return oss.str(); 
} 

int main(int argc, char** argv) 
{ 
    //> Preliminaries. ------------------------------------------------- 
    // Column Major Matrices for setting up the problem. 
    double D_orig[9] = { 
    2, 0, 0, 
    0, -3, 0, 
    0, 0, 0 
    }; 
    double EV_orig[9] = { 
    3./5., 4./5., 0, 
    -4./5., 3./5., 0, 
    0, 0, 1 
    }; 
    double A[9] = { 0,0,0,0,0,0,0,0,0 }; 
    double dummy[9] = { 0,0,0,0,0,0,0,0,0 }; 
    // A = EV D EV' 
    // dummy := D EV' 
    // A := EV dummy 
    cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,3,3,3,1,&D_orig[0],3,&EV_orig[0],3,0,&dummy[0],3); 
    cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,3,3,3,1,&EV_orig[0],3,&dummy[0],3,0,&A[0],3); 
    cout << "Set up the problem building A = EV D EV'" << endl << 
    "EV = [" << endl << matrix2string(3,3,&EV_orig[0]).c_str() << "];" << endl << 
    "D = [" << endl << matrix2string(3,3,&D_orig[0]).c_str() << "];" << endl << 
    "A = [" << endl << matrix2string(3,3,&A[0]).c_str() << "];" << endl; 
    //< ---------------------------------------------------------------- 
    //> Actual eigen value problem. ------------------------------------ 
    char jobz = 'V'; // We want both vectors and values. 
    char uplo = 'L'; // 'L' means lower triangular input A as opposed to 'U'. 
    int N = 3; // Matrix dimension, or as they call it, 'order'. 
    // As stated by example ATriL is unnecessary. Just replace all of its 
    // occurences with plain A and all is well. 
    double ATriL[9] = { A[0], A[1], A[2], A[4], A[5], A[8], 0, 0, 0 }; // Lower Triangle of symmetric A. 
    // Note that it is larger than necessary. It will contain the eigenvectors at the end. 
    int lda = 3; 
    double w[3] = { 0, 0, 0 }; // Container for the eigenvalues. 
    int lwork = 15; // Size of the worker array. Set to (NB+2)*N where NB here is the largest blocksize. 
    // Note, however, that the definition of NB is more complex. 
    // Compare http://ftp.mcs.anl.gov/pub/MINPACK-2/lapack/ilaenv.f 
    double work[lwork]; 
    int info = 0; 

    // "double symmetric eigenvalue problem" I presume. 
    // lapack_int LAPACKE_dsyev(int matrix_order, char jobz, char uplo, lapack_int n, 
    //       double* a, lapack_int lda, double* w); 
    info = LAPACKE_dsyev(LAPACK_COL_MAJOR, jobz, uplo, N, &ATriL[0], lda, &work[0]); 
    // Note that the function takes no parameters lwork and w and that the 
    // eigenvalues appear to be written into work. 
    cout << "Ran dsyev(..) -- presumably 'double symmetric eigenvalue'." << endl << 
    "Lambda: " << work[0] << ", " << work[1] << ", " << work[2] << endl << 
    "EV = [" << endl << matrix2string(3,3,&ATriL[0]) << "];" << endl << 
    "Info: " << info << endl; 
    //< ---------------------------------------------------------------- 
    return EXIT_SUCCESS; 
} 

最後實際的問題：爲什麼結果如此不準確和我能改善他們嗎？

Answer 1

結果是準確的 - 您提供給lapack的數據。不幸的是，你沒有解決你想要的問題。

具有較低與較高三角形部分的部分僅適用於其他（內部？）算法。在你的簡單情況下，你不需要擔心這一點。如果你通過你的矩陣A而不是ATril你應該沒問題。

的詳細信息： 您正在構建double ATriL[9]，使得看上去如

A[0], A[4], 0, 
A[1], A[5], 0, 
A[2], A[8], 0 

到LAPACK。現在當你告訴它使用這個矩陣對稱輸入（char uplo = 'L';）的下三角部分，LAPACK將有效地看到矩陣

A[0], A[1], A[2],  -1.2 2.4 0 
A[1], A[5], A[8], == 2.4 0 0 
A[2], A[8], 0    0 0 0 

，其特徵向量實際上是你得到的人。