[Rd] eigen(symmetric=TRUE) for complex matrices
Berend Hasselman
bhh at xs4all.nl
Wed Jun 19 12:43:45 CEST 2013
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>
More information about the R-devel
mailing list