[Numpy-discussion] inverting and calculating eigenvalues for many small matrices

Daniel Wheeler daniel.wheeler2@gmail....
Tue Aug 16 09:28:20 CDT 2011

On Tue, Aug 16, 2011 at 5:53 AM, David Warde-Farley
<wardefar@iro.umontreal.ca> wrote:
> On 2011-08-15, at 4:11 PM, Daniel Wheeler wrote:
>> One thing that I know I'm doing wrong is
>> reassigning every sub-matrix to a new array. This may not be that
>> costly, but it seems fairly ugly. I wasn't sure how to pass the
>> address of the submatrix to the lapack routines so I'm assigning to a
>> new array and passing that instead.
> It looks like the arrays you're passing are C contiguous. Am I right about this? (I ask because I was under the impression that BLAS/LAPACK routines typically want Fortran-ordered input arrays).

Are you saying that fortran ordered arrays should be passed? The tests
pass when compared against doing numpy equivalents so I don't believe
that its currently broken (maybe suboptimal). There is a transpose and
copy here <http://matforge.org/fipy/browser/sandbox/smallMatrixTools/smt/smt.pyx#L65>
and here <http://matforge.org/fipy/browser/sandbox/smallMatrixTools/smt/smt.pyx#L80>.
I believe that reorders correctly. Maybe I should cast the arrays to
explicit fortran ordering rather than doing that (not sure how)?
However, the transpose and copy doesn't seem to be expensive compared
with the actual lapack routines.

> If your 3D array is also C-contiguous, you should be able to do pointer arithmetic with A.data and B.data. foo.strides[0] will tell you how many bytes you need to move to get to the next element along that axis.

Sounds complicated, but I'll try and figure it out. Thanks for the idea.

> If the 3D array is anything but C contiguous, then I believe the copy is necessary. You should check for that in your Python-visible "solve" wrapper, and make a copy of it that is C contiguous if necessary (check foo.flags.c_contiguous), as this will be likely faster than copying to the same buffer each time in the loop.

The copy is required after the transpose (which is required for
fortran ordering). I'll look into the pointer arithmetic stuff and see
if that helps any. Thanks.

Daniel Wheeler

More information about the NumPy-Discussion mailing list