Dear Simon,
Thank you for reporting this problem with dlarrb. In fact, the problem
has been detected before and has been fixed. A new routine with the
corresponding fix and more advanced features will soon be available in
an upcoming release of LAPACK.
With best wishes,
Osni
Hi,
A user of statistical software R (cran.rproject.org) reported a problem
when
trying to get an eigen decomposition of a particular symmetric matrix. The
problem can be tracked down to lapack routine dsyevr failing to return.
Specifically, routine DLARRB gets stuck in an infinite loop: the reason for
this is that while refining eigenvalues this routine compares an error bound
on an eigenvalue to the difference between two other eigenvalues 
unfortunately this difference has become negative, so the convergence
condition is never met. (The matrix does have high proportion of repeated and
zero eigenvalues, but it's not a problem to decompose it using other lapack
routines.)
Several R developers have replicated the problem (using R, which calls
lapack)
on a variety of platforms (on one machine we got convergence), but we are at
a loss about how to fix it, and are unable to find any reference to such a
problem by internet searching. I have therefore replicated the problem using
a small fortran wrapper program, which loads the offending matrix, and calls
lapack directly.
I used lapack 3 from netlib, with the source patches from
www.netlib.org/lapack/patch
and the blas supplied with lapack. The test suites indicate no problems. I
am
running suse linux 9.3 on a pentium xeon machine, and g77 (gcc 3.3.5). I've
tried turning off all optimization when building the blas, lapack and the
test program, but the problem persists.
I've attached the wrapper program, lap.f (please excuse my fortran) and the
problematic matrix dumped as binary, thematrix.bin, to be read by lap.f.
Any help would be much appreciated!
best,
Simon

c wrapper to call lapack::dsyevr with troublesome matrix
c Program hangs while refining eigenvalues in DLARRB
c compiling lapack, blas and this without optimization
c doesn't help: neither does reducing ABSTOL to the safe minumum
c see also R/eigen.bug.r for further details, including
c code to dumb matrix appropriately.
INTEGER N,LWORK,LIWORK
DOUBLE PRECISION A(72,72),Z(72,72), W(72)
N = 72
call readmat(A,N)
call query_lap(A,N,W,Z,LWORK,LIWORK)
call lap_eig(A,N,W,Z,LWORK,LIWORK)
call hello
STOP
END
SUBROUTINE hello
Print *, 'hello'
END subroutine hello
subroutine readmat(A,n)
c Reads matrix (written with R `writeBin' command)
integer n,ios
double precision A(n,n)
open (2,iostat=ios,file='thematrix.bin',form='unformatted',
$status='old')
read(2) A
close(2)
c Print *, A
end subroutine readmat
subroutine query_lap(A,N,W,Z,LWORK,LIWORK)
c routine to query dsyevr() for optimal workspace size
c returned in LWORK and LIWORK
external DSYEVR
integer N,LWORK,LIWORK
double precision A(N,N),W(N),Z(N,N)
c 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 query_lap
subroutine lap_eig(A,N,W,Z,LWORK,LIWORK)
c routine to call dsyevr() RRR symmetric eigen routine
c A is original matrix, N its dimension.
c on exit W is eigenvalues, Z eigenvectors
external DSYEVR
external DLAMCH
integer N,LWORK,LIWORK
double precision A(N,N), W(N),Z(N,N)
c 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.
c DLAMCH('Safe minimum')
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 subroutine lap_eig

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