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);
}
|