Problem with eigenvectors returned by DSTEGR

Open discussion regarding features, bugs, issues, vendors, etc.

Problem with eigenvectors returned by DSTEGR

Postby snezak » Thu Jan 19, 2012 9:36 pm

Hi,

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);
   
}
snezak
 
Posts: 2
Joined: Thu Jan 19, 2012 9:23 pm

Re: Problem with eigenvectors returned by DSTEGR

Postby Julien Langou » Fri Jan 20, 2012 12:06 pm

I did not review the whole code but, first, the way you allocate you matrix Z to store the eigenvectors is not correct.
The n^2 elements in Z needs to be contiguous. When you have Z[i][j]. You do not necessarily have contiguity.

You need to do:
Code: Select all
   
double *Z;
Z=malloc(N*N*sizeof(double));

Then to access Z(i,j), you get it with Z[(i-1)+(j-1)*N].
(Assuming i and j goes from 1 to N.)

Second, with DSTEGR, you compute the eigenvector of the tridiagonal matrix.
They are not the same as the one of the original matrix. You need to back transform.
(The eigenvalues are the same indeed.)

Your problem is that we do not have DSPEVR driver ... Sorry about that. That would have make your life easier.
It' in the wish list.

For other methods, you could use the drivers: DSPEV, DSPEVD, DSPEVX. This would avoid you quite some lines of code.
This is if you want to work in symmetric packed format but I see that you are allocating the full A, then convert it back to AP a packed matrix ...
So why not using the full A (but do not do A[i][j], do A[i+j*n], see the comment on Z above.)
Then you can readily use the drivers DSYEV, DSYEVD, DSYEVR, DSYEVX.
(And forget about AP.)

Julien
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Problem with eigenvectors returned by DSTEGR

Postby snezak » Sat Jan 21, 2012 8:24 pm

Julien, thanks for you thorough answer!
I'm starting to like LAPACK now that I know how to run it;) Using the full (not tridiagonalized A) with one of the mentioned routines is really handy, DSYEVR was what I was looking for.

Anyway, I found many routines, whose names begin with SORM- that mulitply the orthogonal matrix returned by some algorithm in the form of a product of elementary reflectors with another matrix. Is any of theses suitable to multiply Q returned by DSPTRD (that is stored in arrays AP and TAU) with Z returned by DSTEGR in order to obtain the eigenvalues of the original matrix that was tridiagonalized by DSPTRD?
snezak
 
Posts: 2
Joined: Thu Jan 19, 2012 9:23 pm

Re: Problem with eigenvectors returned by DSTEGR

Postby Julien Langou » Sun Jan 22, 2012 12:56 am

Anyway, I found many routines, whose names begin with SORM- that mulitply the orthogonal matrix returned by some algorithm in the form of a product of elementary reflectors with another matrix. Is any of theses suitable to multiply Q returned by DSPTRD (that is stored in arrays AP and TAU) with Z returned by DSTEGR in order to obtain the eigenvalues of the original matrix that was tridiagonalized by DSPTRD?


Yes definitely. Best is to have a look on how it is done in the drivers themselves. So for an example of use of DORMTR, you can look at DSYEVD. For an example of use of DORGTR, you can look at DSYEV. Julien.
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA


Return to User Discussion

Who is online

Users browsing this forum: Bing [Bot], Yahoo [Bot] and 3 guests

cron