[Numpy-discussion] SVD errors
Pauli Virtanen
pav@iki...
Tue Feb 3 02:46:12 CST 2009
Mon, 02 Feb 2009 18:27:05 -0600, Robert Kern wrote:
> On Mon, Feb 2, 2009 at 18:21, <mtrumpis@berkeley.edu> wrote:
>> Hello list.. I've run into two SVD errors over the last few days. Both
>> errors are identical in numpy/scipy.
>>
>> I've submitted a ticket for the 1st problem (numpy ticket #990).
>> Summary is: some builds of the lapack_lite module linking against
>> system LAPACK (not the bundled dlapack_lite.o, etc) give a
>> "LinAlgError: SVD did not converge" exception on my matrix. This error
>> does occur using Mac's Accelerate framework LAPACK, and a coworker's
>> Ubuntu LAPACK version. It does not seem to happen using ATLAS LAPACK
>> (nor using Octave/Matlab on said Ubuntu)
>
> These are almost certainly issues with the particular implementations of
> LAPACK that you are using. I don't think there is anything we can do
> from numpy or scipy to change this.
Yes, this is almost certainly a LAPACK problem. If in doubt, you can test
it with the following F90 program (be sure to link it against the same
LAPACK as Numpy). Save the matrices with 'np.savetxt("foo.txt", x.ravel
())' and run './test < foo.txt'.
----
program test
implicit none
integer, parameter :: N = 128
double precision :: A(N*N), S(N), U(N,N), Vh(N,N)
double precision, allocatable :: WORK(:)
double precision :: tmp
integer :: IWORK(8*N)
integer :: INFO = 0, LWORK
read(*,*) A
A = reshape(transpose(reshape(A, (/N, N/))), (/ N*N /))
call dgesdd('A', N, N, A, N, S, U, N, Vh, N, tmp, -1, IWORK, INFO)
LWORK = tmp
if (info .ne. 0) stop 'lwork query failed'
write(*,*) 'lwork:', lwork
allocate(WORK(LWORK))
call dgesdd('A', N, N, A, N, S, U, N, Vh, N, WORK, LWORK, IWORK, INFO)
write(*,*) 'info:', INFO
write(*,*) 'min(S):', minval(S)
if (INFO .ne. 0) then
write(*,*) ' -> SVD failed to converge'
end if
if (minval(S) .lt. 0) then
write(*,*) ' -> negative singular value'
end if
end program
----
More information about the Numpy-discussion
mailing list