[SciPy-User] RZ Factorization

[Thomas Unterthiner]

On 2014-09-03 20:02, Sturla Molden wrote:
> Del Citto Francesco <francesco.delcitto at sauber-motorsport.com> wrote:
>
>> Is there any other way of accessing these functions from a clean SciPy installations?
>> Are there plans to include the RZ factorization with a high-level wrapper
>> as done for the QR?
> I often find that I need BLAS and LAPACK functions not exposed by SciPy.
> And so I end up doing exactly what you have done here, except I do not
> modify SciPy but make my own extension module. (I prefer to use Cython
> instead of f2py to better control the behavior of wrapper, though.) Anyhow,
> it raises the question if scipy.linalg.blas and scipy.linalg.lapack perhaps
> should be more complete? Today it mostly exposes the BLAS and LAPACK that
> SciPy needs internally. But perhaps it is tine sto acknowledge that some
> users of SciPy also need other parts of these libraries. If we do what you
> have done here, and extended SciPy's f2py wrappers with missing LAPACK and
> BLAS subroutines, will PRs be accepted? Or there a policy that SciPy should
> only expose the parts of LAPACK and BLAS which it needs internally?
>
> Sturla
>

I also often to run into the situation of needing BLAS/LAPACK 
functionality that is not exposed by SciPy. It would be nice if more/all 
the calls would be wrapped. However, I find it very annoying that the 
f2py wrappers often silently copy and/or transpose data. I understand 
that this is a F-order vs C-order issue, but sometimes e.g. I know that 
my matrices are PSD, yet f2py will copy/transpose anyway. Such behavior 
is annoying and often makes the SciPy wrappers useless for my purposes.

So if talk of extending the LAPACK/BLAS support for SciPy is on the 
table, then I'd like to propose to get rid of that behavior: low-level 
wrappers should not secretly make copies of the data. I'm using 
low-level functions to save speed, yet the copying takes away much of 
the speed gains.


With all that said, I use CFFI for my wrappers, and that typically ends 
up like this:

     # NOTE: both MKL and OpenBLAS (from 0.2.10) know saxpby
     #       but this will fail on older OpenBLAS versions
     __blas_ffi = FFI()
     __blas_ffi.cdef("""
     void cblas_saxpby(const int N, const float alpha,
                       const float *X, const int incX,
                       const float beta, float *Y, const int incY);
         """)
     import os
     c = np.__config__.get_info('blas_opt_info')
     soname = os.path.join(c['library_dirs'][0], "lib"  + 
c['libraries'][0] + '.so')
     __blas_cffi = __blas_ffi.dlopen(soname)


     ### USAGE ###
     px = __blas_ffi.cast("float*", x.ctypes.data)
     py = __blas_ffi.cast("const float*", y.ctypes.data)
     __blas_cffi.cblas_saxpby(x.size, alpha, py, 1, beta, px, 1)


Which I think is fairly straight forward. I think it would be easily 
doable to wrap all of LAPACKE/CBLAS this way (thus getting rid of C- vs 
F-order, since that can be specified by the user) and replace the 
current f2py wrappers with that. If time permits and there is enough 
interested, I could see if I can hack together a prototype/PR for this.

Cheers

Thomas