/*
This file has my implementation of the LAPACK routine dsygv
for C++. This program solves for the eigenvalues and, if
desired, the eigenvectors for a symmetric, generalized eigenvalue
problem of the form Ax = lBx, where A and B are two symmetric
matrices, x is an eigenvector, and l is an eigenvalue. The
program assumes that the upper triangles of the two matrices
are stored.
There are two function calls defined in this header, of the
forms
void dsygv(double **A, double **B, int n, double *E)
void dsygv(double **A, double **B, int n, double *E, double **Evecs)
A: the n by n matrix from the left hand side of the equation
B: the n by n matrix from the right hand side of the equation
n: the order of the square matrices A and B
E: an n-element array to hold the eigenvalues l
Evecs: an n by n matrix to hold the eigenvectors of the
problem, if they are requested.
The function call is defined twice, so that whether or not
eigenvectors are called for the proper version is called.
Scot Shaw
28 September 1999
*/
#include
using namespace std;
void dsygv(double **A, double **B, int n, double *E);
void dsygv(double **A, double **B, int n, double *E, double **Evecs);
double *dsygv_ctof(double **in, int rows, int cols);
void dsygv_ftoc(double *in, double **out, int rows, int cols);
void dsygv_normalize(double **Evecs, int N);
extern "C" void dsygv_(int *itype, char *jobz, char *uplo, int *n,
double *a, int *lda, double *b, int *ldb,
double *w, double *work, int *lwork,
int *info);
void dsygv(double **A, double **B, int n, double *E)
{
char jobz, uplo;
int itype, lda, ldb, lwork, info;
double *a, *b, *work;
itype = 1; /* This sets the type of generalized eigenvalue problem
that we are solving. We have the possible values
1: Ax = lBx
2: ABx = lx
3: BAx = lx */
jobz = 'N'; /* V/N indicates that eigenvectors should/should not
be calculated. */
uplo = 'U'; /* U/L indicated that the upper/lower triangle of the
symmetric matrix is stored. */
lda = n; // The leading dimension of the matrix A
ldb = n; // The leading dimension of the matrix B
lwork = 3*n-1;
work = new double[lwork]; /* The work array to be used by dsygv and
its size. */
a = dsygv_ctof(A, n, lda);
b = dsygv_ctof(B, n, ldb); /* Here we convert the incoming arrays, assumed
to be in double index C form, to Fortran
style matrices. */
dsygv_(&itype, &jobz, &uplo, &n, a, &lda, b, &ldb, E, work, &lwork, &info);
if ((info==0)&&(work[0]>lwork))
cout << "The pre-set lwork value was sub-optimal for the job that\n"
<< "you gave dsygv. The used value was " << lwork
<< " whereas " << work[0] << " is optimal.\n";
delete a;
delete b;
delete work;
}
void dsygv(double **A, double **B, int n, double *E, double **Evecs)
{
char jobz, uplo;
int itype, lda, ldb, lwork, info;
double *a, *b, *work;
itype = 1;
jobz = 'V';
uplo = 'U';
lda = n;
ldb = n;
lwork = 3*n-1;
work = new double[lwork];
a = dsygv_ctof(A, n, lda);
b = dsygv_ctof(B, n, ldb);
// printf(" ref d1 ");
// lwork=-1;
// dsygv_(&itype, &jobz, &uplo, &n, a, &lda, b, &ldb, E, work, &lwork, &info);
// lwork=int(work[0]);
// long ML = 10000;
// if(lwork>ML)
// {
// lwork=ML;
// }
// printf("\n lwork=%d \n",lwork);
dsygv_(&itype, &jobz, &uplo, &n, a, &lda, b, &ldb, E, work, &lwork, &info);
printf("\n dsygv completed \n");
/*
if ((info==0)&&(work[0]>lwork))
cout << "The pre-set lwork value was sub-optimal for the job that\n"
<< "you gave dsygv. The used value was " << lwork
<< " whereas " << work[0] << " is optimal.\n";
*/
dsygv_ftoc(a, Evecs, n, lda);
dsygv_normalize(Evecs, n);
delete a;
delete b;
delete work;
}
double* dsygv_ctof(double **in, int rows, int cols)
{
double *out;
int i, j;
out = new double[rows*cols];
for (i=0; i