I am solving an eigensystem for a symmetric real matrix with DSTEGR routine, which is used from a C program. The eigenvalues returned are exactly as expected, but the eigenvectors aren't. I am interpreting the i-th column of the Z matrix as the eigenvector corresponding to eigenvalue in W[i]. I am not sure whether there is something wrong in my code or are the resulting eigenvectors just misinterpreted. The actual problem being solved is that of finding eigenmodes of a rectengular membrane vibrations, so the eigenvectos are well known.
Here is the code:
- Code: Select all
// cc membrane.c -L /usr/lib -llapack -lblas
#include <stdio.h>
#include <stdlib.h>
#define n 30 // size of grid is n x n
#define N n*n // matrix dimension
#define PI 3.14159265358979323846
int main()
{
// **1** fill matrix
int i,j;
double tmp;
double **A=malloc(N*sizeof(double *));
for(i=0; i<N; i++) A[i]=malloc(N*sizeof(double));
for(i=0; i<N; i++)
for(j=0; j<N; j++) A[i][j]=0.0;
for(i=0; i<N; i++) {
if(i==0 || i==n-1 || i==(n-1)*n || i==n*n-1) tmp=-6.0;
else if( (i>0 && i<n-1) || (i>(n-1)*n && i<n*n-1) ) tmp=-5.0;
else if( i%n==0 || i%n==n-1 ) tmp=-5.0;
else tmp=-4.0;
A[i][i]=tmp;
}
for(i=0; i<n; i++)
for(j=0; j<n; j++) {
if(j<n-1) {
A[i*n+j][i*n+j+1]=1.0;
A[i*n+j+1][i*n+j]=1.0;
}
if(i>0) {
A[i*n+j][(i-1)*n+j]=1.0;
A[(i-1)*n+j][i*n+j]=1.0;
}
}
// A is symmetric
// **2** transform the matrix in order to call a Fortran routine from C
double *AP=malloc((N*(N+1)/2)*sizeof(double));
for(j=0; j<N; j++)
for(i=0; i<=j; i++) AP[i+j*(j+1)/2]=A[i][j]; // column major order
// **3** reduction to tridiagonal form
// DSPTRD http://www.netlib.org/lapack/double/dsptrd.f
char UPLO='U';
int DIM=N;
double D[N], E[N-1], TAU[N-1];
int INFO;
dsptrd_(&UPLO,&DIM,AP,D,E,TAU,&INFO); // all arguments are passed as references
// **4** MRRR algorithm to find eigenvalues and eigenvectors
// DSTEGR http://www.netlib.org/lapack/double/dstegr.f
char JOBZ='V'; // Compute eigenvalues and eigenvectors.
char RANGE='A'; // All eigenvalues will be found.
double E2[N]; for(i=0; i<N-1; i++) E2[i]=E[i];
double VL, VU;
int IL, IU;
double ABSTOL;
int M;
double W[N];
double **Z=malloc(N*sizeof(double *));
for(i=0; i<N; i++) Z[i]=malloc(N*sizeof(double));
for(i=0; i<N; i++) for(j=0; j<N; j++) Z[i][j]=0.0;
int LDZ=N;
int ISUPPZ[2*N];
double WORK[18*N];
int LWORK=18*N;
double IWORK[10*N];
int LIWORK=10*N;
dstegr_( &JOBZ, &RANGE, &DIM, D, E2, &VL, &VU, &IL, &IU, // all arguments are passed as references
&ABSTOL, &M, W, *Z, &LDZ, ISUPPZ, WORK, &LWORK,
IWORK, &LIWORK, &INFO );
// print results
// ** print eigenvalues in reverse order
//for(i=1; i<=N; i++) printf("%d \t %.16f \n",i,-W[N-i]*N/PI/PI); // the eigenvalues returned here are as expected
// ** print eigenvectors
// Z contains the orthonormal eigenvectors
// the i-th column of Z holds the eigenvector associated with W(i)
int k=N-1; // select eigenvector
for(i=0; i<N; i++) {
(i%n==0) ? printf("\n") : printf(" ");
printf("%f",Z[i][k]); // the eigenvector returned here is not as expected!?
}
// deallocation
free(AP);
for(i=0;i<N;i++) {
free(A[i]);
//free(Z[i]);// this gives segfault!?
}
free(A);
//free(Z);
}
