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
|