LAPACK Archives

[Lapack] dsyevr fails to terminate (while calling DLARRB)

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.r-project.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
  


<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