[SciPy-User] scipy+lapack array copy
Bruce Southey
bsouthey@gmail....
Tue Mar 9 08:38:10 CST 2010
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
More information about the SciPy-User
mailing list