Ah! very nice. I did not know that numpy-1.6.1 supports in place 'dot', and neither the fact that you could access the underlying BLAS functions like so. This is pretty neat. Thanks. Now I at least have an idea how the sparse version might work.
If I get time I will probably give numpy-1.6.1 a shot. I already have the MKL libraries thanks to free version of epd for students. On Sat, Mar 26, 2011 at 2:34 PM, Pauli Virtanen <[email protected]> wrote: > > 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! >
_______________________________________________ NumPy-Discussion mailing list [email protected] http://mail.scipy.org/mailman/listinfo/numpy-discussion
