On 19-06-2013, at 10:24, peter dalgaard <pda...@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>
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
