On Tue, Aug 16, 2011 at 5:53 AM, David Warde-Farley <[email protected]> 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 _______________________________________________ NumPy-Discussion mailing list [email protected] http://mail.scipy.org/mailman/listinfo/numpy-discussion
