Are there any 5-argument versions of `mul!` or similar in LinearAlgebra now?

dmbates · September 5, 2018, 3:42pm

The BLAS routines for multiplication and solving linear systems allow for scalar multipliers, usually written alpha and beta, so that the general mutating version of a matrix product is equivalent to

C .= α .* (A * B) + β .* C

I thought I saw 5-argument versions of mul! with a signature like

mul!(C::Matrix{T}, A::Matrix{T}, B::Matrix{T}, alpha::T, beta::T)

but it appears I was hallucinating.

Are there methods for mul!, etc. that would end up calling LinearAlgebra.BLAS.gemm!with alpha and beta arguments?

The reason I want this is because a calculation in a tight inner loop needs to replace C by C - A*B and I would much prefer to do this in-place. If there are alternative ways to accomplish this I would welcome hearing of them.

kristoffer.carlsson · September 5, 2018, 4:37pm

AFAIK, the gemm_wrapper! which is the entry point the the BLAS gemm still hardcodes alpha and beta in its call:

github.com

JuliaLang/julia/blob/8d993569b86608f9588458930a203fa4bc29ad50/stdlib/LinearAlgebra/src/matmul.jl#L458


      
              end
          
          
    if mA == 2 && nA == 2 && nB == 2
                  return matmul2x2!(C,tA,tB,A,B)
              end
              if mA == 3 && nA == 3 && nB == 3
                  return matmul3x3!(C,tA,tB,A,B)
              end
          
          
    if stride(A, 1) == stride(B, 1) == stride(C, 1) == 1 && stride(A, 2) >= size(A, 1) && stride(B, 2) >= size(B, 1) && stride(C, 2) >= size(C, 1)
                  return BLAS.gemm!(tA, tB, one(T), A, B, zero(T), C)
              end
              generic_matmatmul!(C, tA, tB, A, B)
          end
          
          
# blas.jl defines matmul for floats; other integer and mixed precision
          # cases are handled here
          
          
lapack_size(t::AbstractChar, M::AbstractVecOrMat) = (size(M, t=='N' ? 1 : 2), size(M, t=='N' ? 2 : 1))
          
          
function copyto!(B::AbstractVecOrMat, ir_dest::UnitRange{Int}, jr_dest::UnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::UnitRange{Int}, jr_src::UnitRange{Int})

Sparse arrays do have this 5 arg version though, heh:

github.com

JuliaLang/julia/blob/8d993569b86608f9588458930a203fa4bc29ad50/stdlib/SparseArrays/src/linalg.jl#L32


      
          function sppromote(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}
              Tv = promote_type(TvA, TvB)
              Ti = promote_type(TiA, TiB)
              A  = convert(SparseMatrixCSC{Tv,Ti}, A)
              B  = convert(SparseMatrixCSC{Tv,Ti}, B)
              A, B
          end
          
          
# In matrix-vector multiplication, the correct orientation of the vector is assumed.
          
          
function mul!(C::StridedVecOrMat, A::SparseMatrixCSC, B::StridedVecOrMat, α::Number, β::Number)
              A.n == size(B, 1) || throw(DimensionMismatch())
              A.m == size(C, 1) || throw(DimensionMismatch())
              size(B, 2) == size(C, 2) || throw(DimensionMismatch())
              nzv = A.nzval
              rv = A.rowval
              if β != 1
                  β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
              end
              for k = 1:size(C, 2)
                  @inbounds for col = 1:A.n