c – LAPACKE特征解决方案不准确.怎么改进呢？
作者：互联网
我需要Ubuntu下的C eigenproblem求解器.为此我从lapack 3.5.0中给了LAPACKE一个镜头并且实际上设法写了下面的示例程序
例子我是从正交系统和对角矩阵构造的
EV = [
.6, -.8, 0
.8, .6, 0
0, 0, 1
];
D = [
2, 0, 0
0, -3, 0
0, 0, 0
];
通过产生A：= 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;
}
最后实际问题：为什么结果如此不准确,我该怎么做才能改进它们？
解决方法:
结果是准确的 – 对于你给lapack的数据.不幸的是,你没有解决你想要的问题.
具有较低与较高三角形部分的部分仅用于其他一些(内部？)算法.在您的简单情况下,您不必担心这一点.如果你通过矩阵A而不是ATril你应该没事.
更详细：
你正在构建双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
其特征向量实际上就是你得到的那些.
标签：c,lapack 来源： https://codeday.me/bug/20190831/1772840.html