<div dir="auto">The inplace overwriting is done if f2py can forward the original array down to the low level.<div dir="auto"><br></div><div dir="auto">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. </div><div dir="auto"><br></div><div dir="auto">Cython might be a better option for you such that you can pass things around via memory views and cython wrappers of BLAS. </div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Aug 27, 2019, 16:40 Jens Jørgen Mortensen <<a href="mailto:jjmo@dtu.dk">jjmo@dtu.dk</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Sorry!  Stupid me, asking scipy questions on numpy-discussion.  Now <br>
continuing on scipy-user.  Any help is much appreciated.  See short <br>
numpy-discussion thread here: <br>
<a href="https://mail.python.org/pipermail/numpy-discussion/2019-August/079945.html" rel="noreferrer noreferrer" target="_blank">https://mail.python.org/pipermail/numpy-discussion/2019-August/079945.html</a><br>
<br>
<br>
Hi!<br>
<br>
I'm trying to use dgemm, zgemm and friends from scipy.linalg.blas to <br>
multiply matrices efficiently.  As an example, I'd like to do:<br>
<br>
     c += a.dot(b)<br>
<br>
using whatever BLAS scipy is linked to and I want to avoid copies of <br>
large matrices.  This works the way I want it:<br>
<br>
 >>> import numpy as np<br>
 >>> from scipy.linalg.blas import dgemm<br>
 >>> a = np.ones((2, 3), order='F')<br>
 >>> b = np.ones((3, 4), order='F')<br>
 >>> c = np.zeros((2, 4), order='F')<br>
 >>> dgemm(1.0, a, b, 1.0, c, 0, 0, 1)<br>
array([[3., 3., 3., 3.],<br>
        [3., 3., 3., 3.]])<br>
 >>> print(c)<br>
[[3. 3. 3. 3.]<br>
  [3. 3. 3. 3.]]<br>
<br>
but if c is not contiguous, then c is not overwritten:<br>
<br>
 >>> c = np.zeros((7, 4), order='F')[:2, :]<br>
 >>> dgemm(1.0, a, b, 1.0, c, 0, 0, 1)<br>
array([[3., 3., 3., 3.],<br>
        [3., 3., 3., 3.]])<br>
 >>> print(c)<br>
[[0. 0. 0. 0.]<br>
  [0. 0. 0. 0.]]<br>
<br>
Which is also what the docs say, but I think the raw BLAS function dgemm <br>
could do the update of c in-place by setting LDC=7.  See here:<br>
<br>
     <a href="http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f.html" rel="noreferrer noreferrer" target="_blank">http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f.html</a><br>
<br>
Is there a way to call the raw BLAS function from Python?<br>
<br>
I found this capsule thing, but I don't know if there is a way to call <br>
that (maybe using ctypes):<br>
<br>
 >>> from scipy.linalg import cython_blas<br>
 >>> cython_blas.__pyx_capi__['dgemm']<br>
<capsule object "void (char *, char *, int *, int *, int *, <br>
__pyx_t_5scipy_6linalg_11cython_blas_d *, <br>
__pyx_t_5scipy_6linalg_11cython_blas_d *, int *, <br>
__pyx_t_5scipy_6linalg_11cython_blas_d *, int *, <br>
__pyx_t_5scipy_6linalg_11cython_blas_d *, <br>
__pyx_t_5scipy_6linalg_11cython_blas_d *, int *)" at 0x7f06fe1d2ba0><br>
<br>
Best,<br>
Jens Jørgen<br>
_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@python.org" target="_blank" rel="noreferrer">NumPy-Discussion@python.org</a><br>
<a href="https://mail.python.org/mailman/listinfo/numpy-discussion" rel="noreferrer noreferrer" target="_blank">https://mail.python.org/mailman/listinfo/numpy-discussion</a><br>
</blockquote></div>