[SciPy-User] scipy+lapack array copy

[Bruce Southey]

On 8. mars 2010, at 10.10, Skipper Seabold wrote:

 > On Mon, Mar 8, 2010 at 12:47 PM, Paul Anton Letnes
 > <paul.anton.letnes <at> gmail.com> wrote:
 >> Hi everyone.
 >>
 >> I have a numerical project where I link a fortran module to python 
using f2py. The fortran module
basically sets up the left hand side coefficient matrix of a linear 
equation system. The right hand side is
small, so I set that up in python. So far, everything works great.
 >>
 >> Now, when I call scipy-lapack, I am using:
 >> from scipy.linalg.lapack import flapack as lapack
 >> lu, piv, x, info = lapack.zgesv(A, b)
 >> A is created in the fortran part of the program. b is created from 
numpy.zeros(). Both report as
F_CONTIGUOUS by printing A.flags and b.flags.
 >>
 >> The problem is that I'd like A to be large (>1GB), but the lapack 
wrapper (or lapack itself?) copies the
array, doubling the memory use. Aside from being unnecessary, this also 
reduces the practical problem
size I can solve by a factor 2, which of course is very annoying.
 >>
 >> How can I get around this problem?
 >>

If memory usage is a concern then you can directly form the X.T*X matrix 
rather than doing the transpose and multiplication. The advantage is 
that X.T*X has a fixed dimensions regardless of the number of 
observations. Further, you can exploit symmetry and use half-storage 
(there are lapack routines). The disadvantage is that is more complex 
and probably not as fast as doing the multiplication.  If numerical 
stability is a concern then you probably need to use a better solving 
than LU with partial pivoting anyhow.

Bruce