In-place matrix multiplication compatible with 0.5 and 0.6


#1

I want to perform in-place matrix multiplications where the matrices can be sparse or dense. As a simple case, consider evaluating residuals from a regression model

resid .= y - X*β

where resid is a pre-allocated dense vector, y and β are both dense vectors, and X is a StridedArray or a SparseMatrixCSC.

For the sparse X I can call

A_mul_B!(-1, X, β, 1, copy!(resid, y))

but that is not defined for dense X. I can drop down to BLAS.gemv! or BLAS.gemm! for the dense case but I would prefer a method call that applies to both sparse and dense and doesn’t cause allocation or introduce inefficiencies. I would also prefer syntax that applies to both v0.5 and v0.6


#2

Are you OK with adding a dependency ? Then I believe https://github.com/andreasnoack/LinearAlgebra.jl has the A_mul_B! definitions you’re llooking for.


#3

Similar methods should be defined in the dense case. I have most of the definitions in https://github.com/andreasnoack/LinearAlgebra.jl/blob/master/src/juliaBLAS.jl but it requires a bit of work to migrate them to Base and I haven’t had time for that yet. A solution could be to register LinearAlgebra.jl and depend on it but the package doesn’t have documentation yet. For now, the easiest solution might be to copy the methods to your package.


#4

Thanks @andreasnoack, will do.