[Numpy-discussion] inplace dot products

[David Warde-Farley]

Hi Olivier,

There was this idea posted on the Scipy-user list a while back:

	http://projects.scipy.org/pipermail/scipy-user/2008-August/017954.html

but it doesn't look like he got anywhere with it, or even got a  
response.

I just tried it and I observe the same behaviour. A quick look at the  
SciPy sources tells me there is something fishy.

subroutine  
< 
tchar=s,d,c,z>gemm(m,n,k,alpha,a,b,beta,c,trans_a,trans_b,lda,ka,ldb,kb)
   ! c = gemm(alpha,a,b,beta=0,c=0,trans_a=0,trans_b=0,overwrite_c=0)
   ! Calculate C <- alpha * op(A) * op(B) + beta * C

I don't read Fortran very well, but it seems to me as though the  
Fortran prototype doesn't match the python prototype.

I'll poke around a little more, but in summary: there's no numpy- 
sanctioned way to specify an output array for a dot(), AFAIK. This is  
a bit of an annoyance, I agree, though I seem to remember Robert Kern  
offering a fairly compelling argument why it's hard. I just don't know  
what that argument is :)

David


On 18-Feb-09, at 4:06 AM, Olivier Grisel wrote:

> Hi numpist people,
>
> I discovered the ufunc and there ability to compute the results on
> preallocated arrays:
>
>>>> a = arange(10, dtype=float32)
>>>> b = arange(10, dtype=float32) + 1
>>>> c = add(a, b, a)
>>>> c is a
>    True
>>>> a
>    array([  1.,   3.,   5.,   7.,   9.,  11.,  13.,  15.,  17.,
> 19.], dtype=float32)
>
> My questions is : is there a way to have an equivalent for the dot
> product operation:
> I want atlas to build my dot products without allocating a temporary
> array and reuse
> a preallocated array of results. Suppose I have:
>
>>>> results = array((10, 3), dtype=float32)
>>>> W = arange(6, dtype=float32).reshape((2, 3))
>>>> x = arange(20, dtype=float32).reshape((10, 2))
>
> What I want is the equivalent of the following without the
> intermediate call to malloc:
>
>>>> results[:] = dot(x, W)
>
> Any idea? I tried to introspect the various docstring of numpy core
> modules but I
> could not get any lead.
>
> -- 
> Olivier
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion