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

[Pauli Virtanen]

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