[Numpy-discussion] Efficient ?axpy operation without copy (B += a*A)
Pauli Virtanen
pav at iki.fi
Tue Jun 23 14:46:23 EDT 2009
On 2009-06-23, Michael McNeil Forbes <mforbes at physics.ubc.ca> wrote:
> Is there a way of performing vectorized ?axpy (daxpy) operations
> without making copies or dropping into C?
>
> i.e: I want to do
>
> big = (10000,5000)
> A = np.ones(big,dtype=float)
> B = np.ones(big,dtype=float)
> a = 1.5
> B += a*A
I think the only available choice is to use the BLAS routines
from scipy.lib:
>>> from scipy.lib.blas import get_blas_funcs
>>> axpy, = get_blas_funcs(['axpy'], [A, B])
>>> res = axpy(A.ravel(), B.ravel(), A.size, a)
>>> res.base is B
True
>>> B[0,0]
2.5
Works provided A and B are initially in C-order so that ravel()
doesn't create copies. If unsure, it's best to make use of the
return value of axpy and not assume B is modified in-place.
[clip]
> There are exposed blas daxpy operations in scipy, but in the version I
> have (EPD), these also seem to make copies (though recent version seem
> to be fixed by looking at the source.)
I don't see many relevant changes in scipy.lib recently, so I'm
not sure what change you mean by the above.
--
Pauli Virtanen
More information about the NumPy-Discussion
mailing list