[Numpy-discussion] Efficient ?axpy operation without copy (B += a*A)

[Charles R Harris]

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

Now and then I've thought about adding it as a ufunc. It wouldn't be hard.

Chuck