Re: [Numpy-discussion] Fwd: Re: Calling BLAS functions from Python

On 8/27/19 8:04 PM, Ilhan Polat wrote:
> The inplace overwriting is done if f2py can forward the original array 
> down to the low level.
> 
> When it is not contiguous then it has to somehow marshall the view into 
> a compatible array and that is when the inbetween array is formed. And 
> also that array can also be overwritten but that would not be the 
> original view you started with. Hence it is kind of a convenience cost 
> you pay.
> 
> Cython might be a better option for you such that you can pass things 
> around via memory views and cython wrappers of BLAS.

Thanks - I'll look into that.

Jens Jørgen

> On Tue, Aug 27, 2019, 16:40 Jens Jørgen Mortensen <[email protected] 
> <mailto:[email protected]>> wrote:
> 
>     Sorry!  Stupid me, asking scipy questions on numpy-discussion.  Now
>     continuing on scipy-user.  Any help is much appreciated.  See short
>     numpy-discussion thread here:
>     https://mail.python.org/pipermail/numpy-discussion/2019-August/079945.html
> 
> 
>     Hi!
> 
>     I'm trying to use dgemm, zgemm and friends from scipy.linalg.blas to
>     multiply matrices efficiently.  As an example, I'd like to do:
> 
>           c += a.dot(b)
> 
>     using whatever BLAS scipy is linked to and I want to avoid copies of
>     large matrices.  This works the way I want it:
> 
>       >>> import numpy as np
>       >>> from scipy.linalg.blas import dgemm
>       >>> a = np.ones((2, 3), order='F')
>       >>> b = np.ones((3, 4), order='F')
>       >>> c = np.zeros((2, 4), order='F')
>       >>> dgemm(1.0, a, b, 1.0, c, 0, 0, 1)
>     array([[3., 3., 3., 3.],
>              [3., 3., 3., 3.]])
>       >>> print(c)
>     [[3. 3. 3. 3.]
>        [3. 3. 3. 3.]]
> 
>     but if c is not contiguous, then c is not overwritten:
> 
>       >>> c = np.zeros((7, 4), order='F')[:2, :]
>       >>> dgemm(1.0, a, b, 1.0, c, 0, 0, 1)
>     array([[3., 3., 3., 3.],
>              [3., 3., 3., 3.]])
>       >>> print(c)
>     [[0. 0. 0. 0.]
>        [0. 0. 0. 0.]]
> 
>     Which is also what the docs say, but I think the raw BLAS function
>     dgemm
>     could do the update of c in-place by setting LDC=7.  See here:
> 
>     http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f.html
> 
>     Is there a way to call the raw BLAS function from Python?
> 
>     I found this capsule thing, but I don't know if there is a way to call
>     that (maybe using ctypes):
> 
>       >>> from scipy.linalg import cython_blas
>       >>> cython_blas.__pyx_capi__['dgemm']
>     <capsule object "void (char *, char *, int *, int *, int *,
>     __pyx_t_5scipy_6linalg_11cython_blas_d *,
>     __pyx_t_5scipy_6linalg_11cython_blas_d *, int *,
>     __pyx_t_5scipy_6linalg_11cython_blas_d *, int *,
>     __pyx_t_5scipy_6linalg_11cython_blas_d *,
>     __pyx_t_5scipy_6linalg_11cython_blas_d *, int *)" at 0x7f06fe1d2ba0>
> 
>     Best,
>     Jens Jørgen
>     _______________________________________________
>     NumPy-Discussion mailing list
>     [email protected] <mailto:[email protected]>
>     https://mail.python.org/mailman/listinfo/numpy-discussion
> 
> 
> _______________________________________________
> SciPy-User mailing list
> [email protected]
> https://mail.python.org/mailman/listinfo/scipy-user
> 

_______________________________________________
SciPy-User mailing list
[email protected]
https://mail.python.org/mailman/listinfo/scipy-user