[Numpy-discussion] inplace dot products

[Olivier Grisel]

2009/2/20 David Warde-Farley <dwf at cs.toronto.edu>:
> 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 :)
>

Alright, thanks for the reply.

Is there a canonical way /sample code to gain low level access to blas / lapack
atlas routines using ctypes from numpy / scipy code?

I don't mind fixing the dimensions and the ndtype of my array if it can
decrease the memory overhead.

BTW, Robert if your insight on this topic would be very much appreciated.

-- 
Olivier