C implementations for LAPACK functions

[brandonwillard]

While we're still using the C backend as the default, we should be aware that there are no real C implementations for LAPACK function-based Ops.

We should see if there's a quick way to add some C implementations. There's already support for Netlib-providing libraries like MKL and ATLAS, so this could be straightforward in some cases, but we might need more explicit checking in order to support other LAPACK functions.

See #1020 for a discussion mentioning specific, useful LAPACK Ops.

[twiecki]

I don't quite understand what the issue is here, could you add a bit more explanation?

there are no real C implementations for LAPACK function-based Ops

Why is that the case? Because it's not implemented or because LAPACK doesn't provide the right C hooks to develop these Ops?

[brandonwillard]

This issue is that we don't have COp implementations for LAPACK Ops like SVD, Cholesky, QRFull, etc.

Since those Ops are backed by NumPy/SciPy implementations that likely call their respective BLAS/LAPACK subset shared libraries, I'm guessing that we can access those LAPACK routines at the C level like aesara.tensor.blas.GemmRelated subclasses do for BLAS routines.

[jessegrabowski]

Would the numba C implementation of the LAPACK functions (see here) serve as a starting point? I guess they are something like "half-wrappers" around the cython functions in scipy.linalg.cython_lapack in C, I guess to simplify the function signatures and returns, and to cache the functions in memory. The actual functions are implemented using ctypes in numba.np.linalg

I've already been working on LAPACK interfaces of linear algebra functions for numba-scipy, so maybe some of this work could be improved and re-used here? The implementations in aesara.tensor.blas.GemmRelated and aesara.tensor.blas_c , however, look significantly more sophisticated than what is being done in the numba libraries I linked. @brandonwillard I'm curious if you think there is any bridge between the two approaches.

[brandonwillard]

@brandonwillard I'm curious if you think there is any bridge between the two approaches.

The challenge here has more to do with working inside Aesara's old C backend. We could always create Ops that don't have C implementations and only Numba ones. We're moving in that direction anyway, so it's a viable alternative.