eigen(symmetric=TRUE) for complex matrices

On Jun 18, 2013, at 21:49 , peter dalgaard wrote:

On Jun 18, 2013, at 21:23 , Berend Hasselman wrote:

So it seems that the blocked algorithm is the cause of the error and that using the (possibly slow) unblocked algorithm gives the correct result.

Thanks Berend,

The finger seems to be pointing to the internal libRlapack/libRblas. The apparent pattern is that things work when R is linked against external libraries and not when the internal ones are used. So it could be time to start looking for differences between R's lapack module and the original LAPACK code.

Now done and no (relevant) differences found. Hmmm....

There could be compiler issues. It could also be that the internal LAPACK uses a newer version than the external libs and a bug was introduced in between them.

One option could be to bypass R by coding Robin's example in pure Fortran and see whether issues persist when linked against libRlapack, vecLib, and the relevant subset of current LAPACK sources. (The bogus 1.0000 eigenvalue in the 33x33 variant of Robin's example should make it rather easy to spot whether things work or not.)

I have made that pure Fortran program.
I have run it on OS X 10.8.4 with

- libRblas libRlapack
- my private refblas and Lapack
- framework Accelerate
- downloaded zheev.f + dependencies + libRblas
- downloaded zheev.f + dependencies + my refblas

The Fortran program and the shell script to do the compiling and running and the output file follow at the end of this mail.
The shell scripts runs otool -L on each executable to make sure it is linked to the intended libraries.
The fortran program runs zheev with the minimal lwork and the optimal lwork.

Summary of the output:

- libRblas libRlapack  ==> Bad
- my private refblas and Lapack  ==> OK
- framework Accelerate ==> OK
- downloaded zheev.f + dependencies + libRblas ==> Bad
- downloaded zheev.f + dependencies + my refblas ==> OK

It looks like libRblas is the cause of the problem.
I haven't done any further investigation of the matter.

Berend

Output is:

<output>
Running hb with Rlapack and Rblas
./hbRlapack:
	/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRblas.dylib (compatibility version 0.0.0, current version 0.0.0)
	/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRlapack.dylib (compatibility version 3.0.0, current version 3.0.1)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
 lwork=         199
 info=           0
 min eig value= -4.433595556845234E-008
 lwork=        3300
 info=           0
 min eig value= -0.359347003948635     

Running hb with Haslapack and Hasblas (Lapack 3.4.2)
./hbHaslapack:
	/Users/berendhasselman/Library/lib/x86_64/librefblas.dylib (compatibility version 0.0.0, current version 0.0.0)
	/Users/berendhasselman/Library/lib/x86_64/liblapack.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
 lwork=         199
 info=           0
 min eig value= -4.433595556845234E-008
 lwork=        3300
 info=           0
 min eig value= -4.433595559424842E-008

Running hb with -framework Accelerate
./hbveclib:
	/System/Library/Frameworks/Accelerate.framework/Versions/A/Accelerate (compatibility version 1.0.0, current version 4.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
 lwork=         199
 info=           0
 min eig value= -4.433595568797253E-008
 lwork=        3300
 info=           0
 min eig value= -4.433595550683581E-008

Compile downloaded zheev + dependencies
/Users/berendhasselman/Documents/Programming/R/problems/bug-error/eigbug
Running hb with downloaded zheev and Rblas
./hbzheev:
	/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRblas.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
 lwork=         199
 info=           0
 min eig value= -4.433595556845234E-008
 lwork=        3300
 info=           0
 min eig value= -0.359347003948635     

Running hb with downloaded zheev and Hasrefblas
./hbzheev:
	/Users/berendhasselman/Library/lib/x86_64/librefblas.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
 lwork=         199
 info=           0
 min eig value= -4.433595556845234E-008
 lwork=        3300
 info=           0
 min eig value= -4.433595559424842E-008
</output>

Fortran program:
<fortran>
      program test

      integer n
      parameter(n=100)
      complex*16 A(n,n)
      double precision w(n), rwork(3*n-2)

      integer lwork, info

      complex*16 work, wtmp(1)
      allocatable work(:)

      call genmat(A,n)
      lwork = 2*n-1
      allocate(work(lwork))
      print *, "lwork=",lwork
      call zheev("N","L",n,A,n,w,work,lwork,rwork,info)
      print *, "info=",info
      print *, "min eig value=",minval(w)
      deallocate(work)

      call genmat(A,n)
      lwork = -1
      call zheev("N","L",n,A,n,w,wtmp,lwork,rwork,info)
      lwork = real(wtmp(1))
      allocate(work(lwork))
      print *, "lwork=",lwork
      call zheev("N","L",n,A,n,w,work,lwork,rwork,info)
      print *, "info=",info
      print *, "min eig value=",minval(w)

      stop
      end

      subroutine genmat(A,n)
      integer n
      complex*16 A(n,*)
      integer i, j
      double precision tmp

      do i=1,n
          do j=1,n
              tmp = exp(-0.1D0 * (i-j)**2)
              A(i,j) = cmplx(tmp,0.0D0)
          enddo
      enddo

      return
      end
</fortran>


Shell script to do the compiling and running:
<script>
gfortran -c hankinbug.f -o hankinbug.o
gfortran hankinbug.o -L/Library/Frameworks/R.framework/Libraries -lRblas -lRlapack -o hbRlapack
echo "Running hb with Rlapack and Rblas"
otool -L ./hbRlapack
./hbRlapack
echo ""

gfortran hankinbug.o -L${HOME}/Library/lib/x86_64 -lrefblas -llapack -o hbHaslapack
echo "Running hb with Haslapack and Hasblas (Lapack 3.4.2)"
otool -L ./hbHaslapack
./hbHaslapack
echo ""

gfortran hankinbug.o -framework Accelerate -o hbveclib
echo "Running hb with -framework Accelerate"
otool -L ./hbveclib
./hbveclib
echo ""

echo "Compile downloaded zheev + dependencies"
cd zheev
gfortran -c *.f
cd -
gfortran hankinbug.o zheev/*.o -L/Library/Frameworks/R.framework/Libraries -lRblas -o hbzheev
echo "Running hb with downloaded zheev and Rblas"
otool -L ./hbzheev
./hbzheev
echo ""

gfortran hankinbug.o zheev/*.o -L${HOME}/Library/lib/x86_64 -lrefblas -o hbzheev
echo "Running hb with downloaded zheev and Hasrefblas"
otool -L ./hbzheev
./hbzheev
echo ""
</script>