LAPACK Archives

[Lapack] DSYEVR fails to compute eigenvalues and eigenvectors for symmet


Hello,

I can not reproduce your problem on my machine. A C code with your matrix 
embedded is in attachment. The output is copy-pasted below. I have used 
the reference BLAS and LAPACK-3.2.0 from netlib.

Best wishes, 
Julien.

DSYEV
    w[0] = -1.3153e-10
    w[1] = +1.2980e-10
    w[2] = +7.3688e-10
    w[3] = +1.9914e-08
    w[4] = +1.1353e-06
    w[5] = +3.0201e-06
    w[6] = +1.9241e-04
    w[7] = +1.1083e-03
    w[8] = +3.1590e-01
    w[9] = +3.7261e+04
DSYEVR
    w[0] = -1.3123e-10
    w[1] = +1.2949e-10
    w[2] = +7.3688e-10
    w[3] = +1.9914e-08
    w[4] = +1.1353e-06
    w[5] = +3.0201e-06
    w[6] = +1.9241e-04
    w[7] = +1.1083e-03
    w[8] = +3.1590e-01
    w[9] = +3.7261e+04



On Sun, 25 Jan 2009, Stefano Zampini wrote:

Dear Sir/Madam,

I'm a PhD student from University of Milan.
For my research I have to compute eigenvalues and eigenvectors of
symmetric and positive definite matrices, and I need high relative
accuracy with eigenvalues.
When computing with EISPACK subroutine 'rs', I get them; whereas
LAPACK (version 3.1.1 provided by PETSc library version 3.0.0)
subroutine DSYEVR fails to compute the eigenvalues (and therefore the
eigenvectors).

Here is an example: the matrix

A =

 1.0e+003 *

  3.72392911775847   3.72417621094745   3.72441849434491
3.72465594898549   3.72488852275346   3.72511573898875
3.72533600901485   3.72554593008727   3.72574023321641
3.72591379061667
  3.72417621094745   3.72442889425013   3.72467621156890
3.72491825767305   3.72515506843033   3.72538627434410
3.72561047660650   3.72582459799911   3.72602380961074
3.72620331815688
  3.72441849434491   3.72467621156890   3.72492810287318
3.72517435105540   3.72541506239749   3.72564995818663
3.72587780689166   3.72609582239736   3.72629957153625
3.72648456246863
  3.72465594898549   3.72491825767305   3.72517435105540
3.72542448261057   3.72566881657343   3.72590715269008
3.72613840890325   3.72636006260045   3.72656804083824
3.72675812745475
  3.72488852275346   3.72515506843033   3.72541506239749
3.72566881657343   3.72591654406038   3.72615811370300
3.72639257893921   3.72661765868315   3.72682961178555
3.72702447710484
  3.72511573898875   3.72538627434410   3.72564995818663
3.72590715269008   3.72615811370300   3.72640277290883
3.72664031003418   3.72686867105446   3.72708442750212
3.72728385792687
  3.72533600901485   3.72561047660650   3.72587780689166
3.72613840890325   3.72639257893921   3.72664031003418
3.72688090743854   3.72711254372676   3.72733210136668
3.72753609194517
  3.72554593008727   3.72582459799911   3.72609582239736
3.72636006260045   3.72661765868315   3.72686867105446
3.72711254372676   3.72734769918653   3.72757135902932
3.72778027678920
  3.72574023321641   3.72602380961074   3.72629957153625
3.72656804083824   3.72682961178555   3.72708442750212
3.72733210136668   3.72757135902932   3.72779982861853
3.72801453929702
  3.72591379061667   3.72620331815688   3.72648456246863
3.72675812745475   3.72702447710484   3.72728385792687
3.72753609194517   3.72778027678920   3.72801453929702
3.72823624752337

EISPACK code return me

eigenvalues=

 1.329009434728462E-012
 1.236525483686476E-010
 7.374020242881216E-010
 1.991502375246769E-008
 1.135265945273430E-006
 3.020114028996641E-006
 1.924067270613097E-004
 1.108310608877240E-003
 0.315900687690585
  37260.9773916475

whereas LAPACK subroutine DSYEVR returns:

eigenvalues=

  18630.7659839204
  18630.7659839204
  18630.7659839204
  18630.7659839204
  18630.7659839204
  18630.7659839204
  18630.7659839204
  18630.7659839204
  18630.7659839204
  18630.7659839204

and eigenvectors all zeros.


The code I use with DSYEVR is the following (maybe I'm wrong here but
I don't think):

in my main:

     INTEGER N,LWORK,LIWORK
     DOUBLE PRECISION A(10,10),Z(10,10), W(10)
     N = 10
.....
     call query_lap(A,N,W,Z,LWORK,LIWORK)
     call lap_eig(A,N,W,Z,LWORK,LIWORK)
.....
     END

and the subroutines are defined by:

     subroutine query_lap(A,N,W,Z,LWORK,LIWORK)
!     routine to query dsyevr() for optimal workspace size
!     returned in LWORK and LIWORK

     implicit none

     external DSYEVR
     integer N,LWORK,LIWORK
     double precision A(N,N),W(N),Z(N,N)
!     working variables
     INTEGER LDA,IL,IU,M,INFO,LDZ,ISUPPZ(2*N),IWORK(1)
     DOUBLE PRECISION ABSTOL,VL,VU,WORK(1)
     CHARACTER JOBZ,RANGE,UPLO

     JOBZ = 'V'
     RANGE = 'A'
     UPLO = 'U'
     LDA = N
     VL = 0.
     VU = 0.
     ABSTOL = 0.
     IL = 0
     IU = 0
     LDZ = N
     LWORK = -1
     LIWORK = -1
     call DSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL,
    & M, W, Z, LDZ, ISUPPZ, WORK, LWORK,IWORK, LIWORK, INFO )
     LWORK = WORK(1)
     LIWORK = IWORK(1)

     end

     subroutine lap_eig(A,N,W,Z,LWORK,LIWORK)
!     routine to call dsyevr() RRR symmetric eigen routine
!     A is original matrix, N its dimension.
!     on exit W is eigenvalues, Z eigenvectors

     implicit none

     external DSYEVR
     integer N,LWORK,LIWORK
     double precision A(N,N), W(N),Z(N,N)
!     working variables
     INTEGER LDA,IL,IU,M,INFO,LDZ,ISUPPZ(2*N),IWORK(LIWORK)
     DOUBLE PRECISION ABSTOL,VL,VU,WORK(LWORK)
     CHARACTER JOBZ,RANGE,UPLO

     JOBZ = 'V'
     RANGE = 'A'
     UPLO = 'U'
     LDA = N
     VL = 0.
     VU = 0.
     ABSTOL = 0.
     IL = 0
     IU = 0
     LDZ = N
     call DSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL,
    & M, W, Z, LDZ, ISUPPZ, WORK, LWORK,IWORK, LIWORK, INFO )
     LWORK = WORK(1)
     LIWORK = IWORK(1)

     end



--
Stefano
_______________________________________________
Lapack mailing list
Lapack@Domain.Removed
http://lists.cs.utk.edu/listinfo/lapack

-------------- next part --------------
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/time.h>

static int max( int a, int b ){
   if (a>b) return(a); else return(b);
}

static int min( int a, int b ){
   if (a<b) return(a); else return(b);
}


extern void dsyev_( char *jobz, char *uplo, int *n, double *a,
   int *lda, double *w, double *work, int *lwork, int *info );

extern void dsyevr_( char *jobz, char *range, char *uplo,
   int *n, double *a, int *lda, double *vl, double *vu,
   int *il, int *iu, double *abstol, int *m, double *w,
   double *z, int *ldz, double *isuppz, double *work, int *lwork,
   int *iwork, int *liwork, int *info );

int main(int argc, char **argv) {

   int n;
   double *a, *work, *w; 
   int lwork, info;
   double vl, vu;
   int il, iu, m, liwork, *iwork;
   double abstol, *z, *isuppz;

   n = 10;
   a = (double *) malloc( n * n * sizeof( double ) );

a[0] = 3.72392911775847e+003; a[10] = 3.72417621094745e+003; a[20] = 
3.72441849434491e+003; a[30] = 3.72465594898549e+003; a[40] = 
3.72488852275346e+003; a[50] = 3.72511573898875e+003; a[60] = 
3.72533600901485e+003; a[70] = 3.72554593008727e+003; a[80] = 
3.72574023321641e+003; a[90] = 3.72591379061667e+003;
a[1] = 3.72417621094745e+003; a[11] = 3.72442889425013e+003; a[21] = 
3.72467621156890e+003; a[31] = 3.72491825767305e+003; a[41] = 
3.72515506843033e+003; a[51] = 3.72538627434410e+003; a[61] = 
3.72561047660650e+003; a[71] = 3.72582459799911e+003; a[81] = 
3.72602380961074e+003; a[91] = 3.72620331815688e+003;
a[2] = 3.72441849434491e+003; a[12] = 3.72467621156890e+003; a[22] = 
3.72492810287318e+003; a[32] = 3.72517435105540e+003; a[42] = 
3.72541506239749e+003; a[52] = 3.72564995818663e+003; a[62] = 
3.72587780689166e+003; a[72] = 3.72609582239736e+003; a[82] = 
3.72629957153625e+003; a[92] = 3.72648456246863e+003;
a[3] = 3.72465594898549e+003; a[13] = 3.72491825767305e+003; a[23] = 
3.72517435105540e+003; a[33] = 3.72542448261057e+003; a[43] = 
3.72566881657343e+003; a[53] = 3.72590715269008e+003; a[63] = 
3.72613840890325e+003; a[73] = 3.72636006260045e+003; a[83] = 
3.72656804083824e+003; a[93] = 3.72675812745475e+003;
a[4] = 3.72488852275346e+003; a[14] = 3.72515506843033e+003; a[24] = 
3.72541506239749e+003; a[34] = 3.72566881657343e+003; a[44] = 
3.72591654406038e+003; a[54] = 3.72615811370300e+003; a[64] = 
3.72639257893921e+003; a[74] = 3.72661765868315e+003; a[84] = 
3.72682961178555e+003; a[94] = 3.72702447710484e+003;
a[5] = 3.72511573898875e+003; a[15] = 3.72538627434410e+003; a[25] = 
3.72564995818663e+003; a[35] = 3.72590715269008e+003; a[45] = 
3.72615811370300e+003; a[55] = 3.72640277290883e+003; a[65] = 
3.72664031003418e+003; a[75] = 3.72686867105446e+003; a[85] = 
3.72708442750212e+003; a[95] = 3.72728385792687e+003;
a[6] = 3.72533600901485e+003; a[16] = 3.72561047660650e+003; a[26] = 
3.72587780689166e+003; a[36] = 3.72613840890325e+003; a[46] = 
3.72639257893921e+003; a[56] = 3.72664031003418e+003; a[66] = 
3.72688090743854e+003; a[76] = 3.72711254372676e+003; a[86] = 
3.72733210136668e+003; a[96] = 3.72753609194517e+003;
a[7] = 3.72554593008727e+003; a[17] = 3.72582459799911e+003; a[27] = 
3.72609582239736e+003; a[37] = 3.72636006260045e+003; a[47] = 
3.72661765868315e+003; a[57] = 3.72686867105446e+003; a[67] = 
3.72711254372676e+003; a[77] = 3.72734769918653e+003; a[87] = 
3.72757135902932e+003; a[97] = 3.72778027678920e+003;
a[8] = 3.72574023321641e+003; a[18] = 3.72602380961074e+003; a[28] = 
3.72629957153625e+003; a[38] = 3.72656804083824e+003; a[48] = 
3.72682961178555e+003; a[58] = 3.72708442750212e+003; a[68] = 
3.72733210136668e+003; a[78] = 3.72757135902932e+003; a[88] = 
3.72779982861853e+003; a[98] = 3.72801453929702e+003;
a[9] = 3.72591379061667e+003; a[19] = 3.72620331815688e+003; a[29] = 
3.72648456246863e+003; a[39] = 3.72675812745475e+003; a[49] = 
3.72702447710484e+003; a[59] = 3.72728385792687e+003; a[69] = 
3.72753609194517e+003; a[79] = 3.72778027678920; a[89] = 3.72801453929702e+003; 
a[99] = 3.72823624752337e+003;

   w = (double *) malloc( n * sizeof( double ) );
   work = (double *) malloc( 1 * sizeof( double ) );
   lwork = -1;
   dsyev_( "V", "U", &n, a, &n, w, work, &lwork, &info );
   lwork = (int) work[0];
   work = (double *) malloc( lwork * sizeof( double ) );
   dsyev_( "V", "U", &n, a, &n, w, work, &lwork, &info );
   free(work);
   printf("DSYEV\n");
   printf("   w[%d] = %+10.4e\n",0,w[0]);
   printf("   w[%d] = %+10.4e\n",1,w[1]);
   printf("   w[%d] = %+10.4e\n",2,w[2]);
   printf("   w[%d] = %+10.4e\n",3,w[3]);
   printf("   w[%d] = %+10.4e\n",4,w[4]);
   printf("   w[%d] = %+10.4e\n",5,w[5]);
   printf("   w[%d] = %+10.4e\n",6,w[6]);
   printf("   w[%d] = %+10.4e\n",7,w[7]);
   printf("   w[%d] = %+10.4e\n",8,w[8]);
   printf("   w[%d] = %+10.4e\n",9,w[9]);
   free(w);
a[0] = 3.72392911775847e+003; a[10] = 3.72417621094745e+003; a[20] = 
3.72441849434491e+003; a[30] = 3.72465594898549e+003; a[40] = 
3.72488852275346e+003; a[50] = 3.72511573898875e+003; a[60] = 
3.72533600901485e+003; a[70] = 3.72554593008727e+003; a[80] = 
3.72574023321641e+003; a[90] = 3.72591379061667e+003;
a[1] = 3.72417621094745e+003; a[11] = 3.72442889425013e+003; a[21] = 
3.72467621156890e+003; a[31] = 3.72491825767305e+003; a[41] = 
3.72515506843033e+003; a[51] = 3.72538627434410e+003; a[61] = 
3.72561047660650e+003; a[71] = 3.72582459799911e+003; a[81] = 
3.72602380961074e+003; a[91] = 3.72620331815688e+003;
a[2] = 3.72441849434491e+003; a[12] = 3.72467621156890e+003; a[22] = 
3.72492810287318e+003; a[32] = 3.72517435105540e+003; a[42] = 
3.72541506239749e+003; a[52] = 3.72564995818663e+003; a[62] = 
3.72587780689166e+003; a[72] = 3.72609582239736e+003; a[82] = 
3.72629957153625e+003; a[92] = 3.72648456246863e+003;
a[3] = 3.72465594898549e+003; a[13] = 3.72491825767305e+003; a[23] = 
3.72517435105540e+003; a[33] = 3.72542448261057e+003; a[43] = 
3.72566881657343e+003; a[53] = 3.72590715269008e+003; a[63] = 
3.72613840890325e+003; a[73] = 3.72636006260045e+003; a[83] = 
3.72656804083824e+003; a[93] = 3.72675812745475e+003;
a[4] = 3.72488852275346e+003; a[14] = 3.72515506843033e+003; a[24] = 
3.72541506239749e+003; a[34] = 3.72566881657343e+003; a[44] = 
3.72591654406038e+003; a[54] = 3.72615811370300e+003; a[64] = 
3.72639257893921e+003; a[74] = 3.72661765868315e+003; a[84] = 
3.72682961178555e+003; a[94] = 3.72702447710484e+003;
a[5] = 3.72511573898875e+003; a[15] = 3.72538627434410e+003; a[25] = 
3.72564995818663e+003; a[35] = 3.72590715269008e+003; a[45] = 
3.72615811370300e+003; a[55] = 3.72640277290883e+003; a[65] = 
3.72664031003418e+003; a[75] = 3.72686867105446e+003; a[85] = 
3.72708442750212e+003; a[95] = 3.72728385792687e+003;
a[6] = 3.72533600901485e+003; a[16] = 3.72561047660650e+003; a[26] = 
3.72587780689166e+003; a[36] = 3.72613840890325e+003; a[46] = 
3.72639257893921e+003; a[56] = 3.72664031003418e+003; a[66] = 
3.72688090743854e+003; a[76] = 3.72711254372676e+003; a[86] = 
3.72733210136668e+003; a[96] = 3.72753609194517e+003;
a[7] = 3.72554593008727e+003; a[17] = 3.72582459799911e+003; a[27] = 
3.72609582239736e+003; a[37] = 3.72636006260045e+003; a[47] = 
3.72661765868315e+003; a[57] = 3.72686867105446e+003; a[67] = 
3.72711254372676e+003; a[77] = 3.72734769918653e+003; a[87] = 
3.72757135902932e+003; a[97] = 3.72778027678920e+003;
a[8] = 3.72574023321641e+003; a[18] = 3.72602380961074e+003; a[28] = 
3.72629957153625e+003; a[38] = 3.72656804083824e+003; a[48] = 
3.72682961178555e+003; a[58] = 3.72708442750212e+003; a[68] = 
3.72733210136668e+003; a[78] = 3.72757135902932e+003; a[88] = 
3.72779982861853e+003; a[98] = 3.72801453929702e+003;
a[9] = 3.72591379061667e+003; a[19] = 3.72620331815688e+003; a[29] = 
3.72648456246863e+003; a[39] = 3.72675812745475e+003; a[49] = 
3.72702447710484e+003; a[59] = 3.72728385792687e+003; a[69] = 
3.72753609194517e+003; a[79] = 3.72778027678920; a[89] = 3.72801453929702e+003; 
a[99] = 3.72823624752337e+003;

   w = (double *) malloc( n * sizeof( double ) );
   z = (double *) malloc( n * n * sizeof( double ) );
   isuppz = (double *) malloc( 2 * n * sizeof( double ) );
   abstol =0e+00;
   lwork = -1;
   liwork = -1;
   work = (double *) malloc( 1 * sizeof( double ) );
   iwork = (int *) malloc( 1 * sizeof( int ) );
   dsyevr_( "V", "A", "U", &n, a, &n, &vl, &vu, &il, &iu, &abstol, &m, w,
      z, &n, isuppz, work, &lwork, iwork, &liwork, &info );
   lwork = (int) work[0];
   liwork = (int) iwork[0];
   free(work);
   free(iwork);
   work = (double *) malloc( lwork * sizeof( double ) );
   iwork = (int *) malloc( liwork * sizeof( int ) );
   dsyevr_( "V", "A", "U", &n, a, &n, &vl, &vu, &il, &iu, &abstol, &m, w,
      z, &n, isuppz, work, &lwork, iwork, &liwork, &info );
   free(isuppz);
   free(work);
   free(iwork);
   free(z);
   printf("DSYEVR\n");
   printf("   w[%d] = %+10.4e\n",0,w[0]);
   printf("   w[%d] = %+10.4e\n",1,w[1]);
   printf("   w[%d] = %+10.4e\n",2,w[2]);
   printf("   w[%d] = %+10.4e\n",3,w[3]);
   printf("   w[%d] = %+10.4e\n",4,w[4]);
   printf("   w[%d] = %+10.4e\n",5,w[5]);
   printf("   w[%d] = %+10.4e\n",6,w[6]);
   printf("   w[%d] = %+10.4e\n",7,w[7]);
   printf("   w[%d] = %+10.4e\n",8,w[8]);
   printf("   w[%d] = %+10.4e\n",9,w[9]);
   free(w);
   free(a);

   exit(0);

}

<Prev in Thread] Current Thread [Next in Thread>


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or