Re: [Numpy-discussion] problem, bug?

[Jonathan B. <bo...@ac...>]

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/