[Numpy-discussion] Array views

[Pauli Virtanen]

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