eigen(symmetric=TRUE) for complex matrices
On Jun 19, 2013, at 12:43 , Berend Hasselman wrote:
On 19-06-2013, at 10:24, peter dalgaard <pdalgd at gmail.com> wrote:
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>
Thanks. I think I have it nailed down now. The culprit was indeed in our reference BLAS (I had only checked the LAPACK code), cmplxblas.f to be specific. Revision 53001 had a number of IF statements being commented out, but two of the changes looked like this:
@@ -1561,7 +1561,7 @@
C( J, J ) = DBLE( C( J, J ) )
END IF
DO 170 L = 1, K
- IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
+c IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
$ THEN
TEMP1 = ALPHA*DCONJG( B( J, L ) )
TEMP2 = DCONJG( ALPHA*A( J, L ) )
Notice that the continuation line was NOT commented out. So FORTRAN, being what it is, continues the line before the comment and parses it as
DO 170 L = 1, KTHEN
with KTHEN uninitialized! and things go downhill from there.
(The uninitialized variable was actually hinted at in PR14964 and the fact that I could get one of my builds to segfault also helped.)
Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com