[SciPy-User] Using GEMM with specified output-array

[Sturla Molden]

Thomas Unterthiner <thomas_unterthiner at web.de> wrote:

> However, this gives me "failed in converting 2nd keyword `c' of 
> _fblas.sgemm to C/Fortran array". I assume this is because 
> Z.flags['OWNDATA'] is False, due to the transpose. Is that correct?

I am not sure. Only Pearu knows all the details about f2py.

In general f2py is written to produce correct results, not generate the
most efficient interface code.
 
> The problem is that other parts of my program expect all my matrices in 
> row-major order, so  sing `order="F"`  is not really an option for me. 
> Is there any way around this?

Call cblas_sgemm from Cython or use numpy.dot would be the easiest.

If you need to grab the BLAS SGEMM function from SciPy, sgemm._cpointer is
a PyCObject that encapsulates a function pointer to the Fortran SGEMM
subroutine. You can e.g. use this in Cython to call SGEMM without the f2py
layer. This requires that you know the ABI of the Fortran compiler used to
build BLAS.

I would also like to add that if you use the Fortran 2003 ISO C bindings,
it is almost impossible to avoid that the compiler generates temporary
copies of your arrays. If you cast C pointers to Fortran pointers with
f_c_pointer, the aliasing rules in Fortran mandates that copy-in copy-out
is used when you call SGEMM with the Fortran pointers as dummy variables.

Sturla