[Numpy-discussion] problem, bug?
Jonathan Bober
bober at acm.cs.nyu.edu
Thu Mar 25 10:44:06 EST 2004
I've looked a little more at this and find that it seems to be a bug the
the lapack_lite code, but I really can't decipher it at the moment. I
believe I have a fix in the python code that works, though.
An offlist email from Arnd Baecker pointed me to:
"LWORK (input) INTEGER
The dimension of the array WORK. LWORK >= 1.
[...]
If JOBZ = 'S' or 'A' LWORK >=
3*min(M,N)*min(M,N) +
max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)). For good performance,
LWORK should generally be larger.
If LWORK = -1 but other input arguments are legal, WORK(1) returns
the optimal LWORK."
The singular_value_decomposition function first calls
lapack_lite2.dgesdd with lwork=-1 in order to get the optimal lwork
value, and then calls it with the optimal lwork value returned. However,
an invalid lwork value is being returned sometimes perhaps when
4*min(M,N)*min(M,N)+4*min(M,N) > max(M,N).
Because singular_value_decomposition always uses JOBZ = 'S' or 'A', I
changed
lwork = 3*min(m,n)*min(m,n) +
max(max(m,n),4*min(m,n)*min(m,n)+4*min(m,n)) + 500
lwork = 1
work = num.zeros((lwork,), t)
results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, work, -1,
iwork, 0)
lwork = int(work[0])
print 'lwork = ' , lwork
work = num.zeros((lwork,), t)
results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, work,
lwork, iwork, 0)
to
lwork = 3*min(m,n)*min(m,n) +
max(max(m,n),4*min(m,n)*min(m,n)+4*min(m,n)) + 500
work = num.zeros((lwork,), t)
results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, work,
lwork, iwork, 0)
and things seem to be working now. Before, doing
A = ones((60, 100), type='Float64')
x,y,z = svd(A)
resulted in the error, while
A = one((100,100), type='Float64')
x,y,z = svd(A)
did not.
Now both seem to work, and I am currently running
for x in range(1,1000):
for y in range(1,1000):
A = ones((x,y), type='Float64')
a,b,c = svd(A)
which hasn't failed yet, but will of course take quite a long time. If
this all works, I'll submit a patch at some point, perhaps, though the
real fix would be a fix in the lapack_lite module.
--
Jon Bober
http://acm.cs.nyu.edu/~bober/
More information about the NumPy-Discussion
mailing list