On Sat, 26 Mar 2011 19:13:43 +0000, Pauli Virtanen wrote:
[clip]
> If you want to have control over temporaries, you can make use of the
> out= argument of ufuncs (`numpy.dot` will gain it in 1.6.1 --- you can
> call LAPACK routines from scipy.lib in the meantime, if your data is in
> Fortran order).
Like so:
# Fortran-order for efficient DGEMM -- each column must be contiguous
A = np.random.randn(4,4).copy('F')
b = np.random.randn(4,10).copy('F')
def updater(b, col_idx):
# This will work in Numpy 1.6.1
dot(A, b[:,col_idx].copy(), out=b[:,col_idx])
In the meantime you can do
A = np.random.randn(4,4).copy('F')
b = np.random.randn(4,10).copy('F')
from scipy.lib.blas import get_blas_funcs
gemm, = get_blas_funcs(['gemm'], [A, b]) # get correct type func
def updater(b, col_idx):
bcol = b[:,col_idx]
c = gemm(1.0, A, bcol.copy(), 0.0, bcol, overwrite_c=True)
assert c is bcol # check that it didn't make copies!
Note that DGEMM and `dot` cannot do in-place multiplication -- at least
the BLAS library I have fails when the B and C arguments point to the
same memory, so you'll anyway end up with one temporary. (This has
nothing to do with Scipy -- same occurs in Fortran).
--
Pauli Virtanen
_______________________________________________
NumPy-Discussion mailing list
[email protected]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
