[SciPy-User] numpy svd

[Sturla Molden]

Den 23.03.2011 17:58, skrev Peter Spuhler:
> I've been porting some IDL code over to scipy and ran into a problem 
> with linalg.svd()
> The following would give me an error message running in a 32-bit 
> environment (epd 7.0.2).
> >>>numpy.linalg.svd(numpy.ones((10,10000)))
> Traceback (most recent call last):
>   File "C:\Program Files (x86)\Wing IDE 
> 4.0\src\debug\tserver\_sandbox.py", line 1, in <module>
>     # Used internally for debug sandbox under external interpreter
>   File "C:\Python27\Lib\site-packages\numpy\linalg\linalg.py", line 
> 1324, in svd
>     vt = vt.transpose().astype(result_t)
> MemoryError:
>
> The same function in 32-bit IDL seems to work fine (as well as in 
> Matlab and Mathematica)
> IDL>la_svd,dblarr(10,10000)+1,w,u,v
>
> They both use the gesdd lapack function on the backend.
> Why would the numpy routine have problems with this calculation when 
> the seemingly similar calculation works fine using IDL or Matlab?
>

Your largest returned array is 762 MB, so it should not run out of 
memory. But the transposition creates one or two temporary arrays (one 
for transpose, an possibly another for astype), which is above or close 
to the 2 GB limit for 32-bit systems. Try to pass a Fortran contiguous 
array to linalg.svd to suppress this. If that does not help, you can 
call LAPACK gesdd directly (see below). Remember to pass a Fortran 
contiguous array to avoid a temporary copy by f2py.

import scipy as sp
import scipy.linalg
from sp.linalg.flapack import dgesdd

64-bit will probably also solve the problem, if you have enough RAM.

Sturla