[Numpy-discussion] Efficient ?axpy operation without copy (B += a*A)
Charles R Harris
charlesr.harris at gmail.com
Tue Jun 23 15:31:03 EDT 2009
On Tue, Jun 23, 2009 at 12:46 PM, Pauli Virtanen<pav at iki.fi> wrote:
> 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
> 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.
>> 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.
Now and then I've thought about adding it as a ufunc. It wouldn't be hard.
More information about the NumPy-Discussion