I need a C eigenproblem solver under Ubuntu. To this end, I gave LAPACKE from a lapack 3.5.0 shot and was actually able to write a lower sample program, which is an example that I built from an orthogonal system and a diagonal matrix
EV = [
.6, -.8, 0
.8, .6, 0
0, 0, 1
];
D = [
2, 0, 0
0, -3, 0
0, 0, 0
];
producing A: = EV D EV '.
While the bottom program is working fine, the results are strangely inaccurate. Here's the end of the output:
...
Lambda: -3.07386, 0, 1.87386
EV = [
-0.788205 0 -0.615412
0.615412 0 -0.788205
-0 1 0
];
Info: 0
As documentation I used www.physics.orst.edu/~rubin/nacphy/lapack/routines/dsyev.html
My full source:
#include <iostream>
#include <string>
#include <sstream>
#include <cstdlib>
#include <cblas.h>
#include <lapacke.h>
using namespace std;
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)
{
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 };
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;
char jobz = 'V';
char uplo = 'L';
int N = 3;
double ATriL[9] = { A[0], A[1], A[2], A[4], A[5], A[8], 0, 0, 0 };
int lda = 3;
double w[3] = { 0, 0, 0 };
int lwork = 15;
double work[lwork];
int info = 0;
info = LAPACKE_dsyev(LAPACK_COL_MAJOR, jobz, uplo, N, &ATriL[0], lda, &work[0]);
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;
}
Finally, relevant questions: why are the results so inaccurate and what can I do to improve them?