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

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Ilhan Polat
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. 

On Tue, Aug 27, 2019, 16:40 Jens Jørgen Mortensen <[hidden email]> 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
[hidden email]
https://mail.python.org/mailman/listinfo/numpy-discussion

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

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

Jens Jørgen Mortensen
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 <[hidden email]
> <mailto:[hidden email]>> 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
>     [hidden email] <mailto:[hidden email]>
>     https://mail.python.org/mailman/listinfo/numpy-discussion
>
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> https://mail.python.org/mailman/listinfo/scipy-user
>

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