minor fixes

Author: dkarrasch

src/LinearMaps.jl

@@ -47,20 +47,17 @@ Base.ndims(::LinearMap) = 2
 Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
 Base.length(A::LinearMap) = size(A)[1] * size(A)[2]
 
-# check dimension consistency for right multiply: y = A*x and Y = A*X
-function check_dim_mul(Y::AbstractVecOrMat, A::LinearMap, X::AbstractVecOrMat)
+# check dimension consistency for multiplication C = A*B
+function check_dim_mul(C, A, B)
     # @info "checked vector dimensions" # uncomment for testing
-    m, n = size(A)
-    (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) ||
-        throw(DimensionMismatch("mul! $(size(X)) $(size(Y)) $(size(A))"))
-    return nothing
-end
-
-# check dimension consistency for left multiply: X = Y*A (e.g., Y=y')
-function check_dim_mul(X::AbstractVecOrMat, Y::AbstractVecOrMat, A::LinearMap)
-    m, n = size(A)
-    (n == size(X, 2) && m == size(Y, 2) && size(Y, 1) == size(X, 1)) ||
-        throw(DimensionMismatch("left mul!"))
+    mA, nA = size(A) # A always has two dimensions
+    mB, nB = size(B, 1), size(B, 2)
+    if mB != nA
+        throw(DimensionMismatch("left factor has dimensions ($mA,$nA), right factor has dimensions ($mB,$nB)"))
+    end
+    if size(C, 1) != mA || size(C, 2) != nB
+        throw(DimensionMismatch("result has dimensions $(size(C)), needs ($mA,$nB)"))
+    end
     return nothing
 end
 

src/left.jl

@@ -15,24 +15,28 @@ Base.:(*)(y::AdjointAbsVec, A::LinearMap) = adjoint(*(A', y'))
 Base.:(*)(y::TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y))
 
 # mul!(x, y', A)
-LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec, A::LinearMap) =
-    mul!(x, y, A, true, false)
-
-# todo: not sure if we need bounds checks and propagate inbounds stuff here
+Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::AdjointAbsVec, A::LinearMap)
+    @boundscheck check_dim_mul(x, y, A)
+    @inbounds mul!(adjoint(x), A', y')
+    return adjoint(x)
+end
 
-# mul!(x, y', A, α, β)
-# the key here is that "adjoint" and "transpose" are lazy "views"
-function LinearAlgebra.mul!(x::AdjointAbsVec, y::AdjointAbsVec,
+Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::AdjointAbsVec,
         A::LinearMap, α::Number, β::Number)
-    check_dim_mul(x, y, A)
-    mul!(adjoint(x), A', y', α, β)
+    @boundscheck check_dim_mul(x, y, A)
+    @inbounds mul!(adjoint(x), A', y', α, β)
     return adjoint(x)
 end
 
-# mul!(x, transpose(y), A, α, β)
-function LinearAlgebra.mul!(x::TransposeAbsVec, y::TransposeAbsVec,
+Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::TransposeAbsVec, A::LinearMap)
+    @boundscheck check_dim_mul(x, y, A)
+    @inbounds mul!(transpose(x), transpose(A), transpose(y))
+    return transpose(x)
+end
+
+Base.@propagate_inbounds function LinearAlgebra.mul!(x::AbstractMatrix, y::TransposeAbsVec,
          A::LinearMap, α::Number, β::Number)
-    check_dim_mul(x, y, A)
-    mul!(transpose(x), transpose(A), transpose(y), α, β)
+    @boundscheck check_dim_mul(x, y, A)
+    @inbounds mul!(transpose(x), transpose(A), transpose(y), α, β)
     return transpose(x)
 end

test/left.jl

@@ -5,7 +5,7 @@ using LinearAlgebra: mul!
 
 
 function left_tester(L::LinearMap{T}) where {T}
-    M,N = size(L)
+    M, N = size(L)
     A = Matrix(L)
 
     x = rand(T, N)
@@ -26,12 +26,10 @@ function left_tester(L::LinearMap{T}) where {T}
 
     b1 = transpose(y) * A
     b2 = similar(b1)
-@static if VERSION ≥ v"1.4.0"
-    mul!(b2, transpose(y), L) # 3-arg # todo: fails in Julia 1.0
+    mul!(b2, transpose(y), L) # 3-arg
     @test b1 ≈ b2
     mul!(b2, transpose(y), L, true, false) # 5-arg
     @test b1 ≈ b2
-end # version
 
     true
 end