diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl
index 373cf5b5..10ae2fc4 100644
--- a/src/LinearMaps.jl
+++ b/src/LinearMaps.jl
@@ -72,33 +72,15 @@ function Base.:(*)(A::LinearMap, x::AbstractVector)
     size(A, 2) == length(x) || throw(DimensionMismatch("mul!"))
     return @inbounds mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
 end
-function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=true, β::Number=false)
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector,
+                    α::Number=true, β::Number=false)
     @boundscheck check_dim_mul(y, A, x)
-    if isone(α)
-        iszero(β) && (A_mul_B!(y, A, x); return y)
-        isone(β) && (y .+= A * x; return y)
-        # β != 0, 1
-        rmul!(y, β)
-        y .+= A * x
-        return y
-    elseif iszero(α)
-        iszero(β) && (fill!(y, zero(eltype(y))); return y)
-        isone(β) && return y
-        # β != 0, 1
-        rmul!(y, β)
-        return y
-    else # α != 0, 1
-        iszero(β) && (A_mul_B!(y, A, x); rmul!(y, α); return y)
-        isone(β) && (y .+= rmul!(A * x, α); return y)
-        # β != 0, 1
-        rmul!(y, β)
-        y .+= rmul!(A * x, α)
-        return y
-    end
+    _muladd!(MulStyle(A), y, A, x, α, β)
 end
 # the following is of interest in, e.g., subspace-iteration methods
-Base.@propagate_inbounds function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix, α::Number=true, β::Number=false)
-    @boundscheck check_dim_mul(Y, A, X)
+function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix, α::Number=true, β::Number=false)
+    (size(Y, 1) == size(A, 1) && size(X, 1) == size(A, 2) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!"))
+    # TODO: reuse intermediate z if necessary
     @inbounds @views for i = 1:size(X, 2)
         mul!(Y[:, i], A, X[:, i], α, β)
     end
@@ -109,13 +91,55 @@ Base.@propagate_inbounds function LinearAlgebra.mul!(Y::AbstractMatrix, A::Linea
     return Y
 end
 
-A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector)  = mul!(y, A, x)
-At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
-Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)
+# multiplication helper functions
+@inline _muladd!(::FiveArg, y, A::MapOrMatrix, x, α, β) = mul!(y, A, x, α, β)
+@inline _muladd!(::ThreeArg, y, A::MapOrMatrix, x, α, β) =
+    iszero(β) ? muladd!(y, A, x, α, β, nothing) : muladd!(y, A, x, α, β, similar(y))
+@inline _muladd!(::FiveArg, y, A::MapOrMatrix, x, α, β, _) = mul!(y, A, x, α, β)
+@inline _muladd!(::ThreeArg, y, A::MapOrMatrix, x, α, β, z) = muladd!(y, A, x, α, β, z)
 
+"""
+    muladd!(y, A, x, α, β, z)
+
+Compute 5-arg multiplication for `LinearMap`s that do not have a 5-arg `mul!`,
+but only a 3-arg `mul!`. For storing the intermediate `A*x*α`, a vector (or
+matrix, as appropriate) `z` is (already) provided.
+"""
+@inline function muladd!(y, A::MapOrMatrix, x, α, β, z)
+    if iszero(α)
+        iszero(β) && (fill!(y, zero(eltype(y))); return y)
+        isone(β) && return y
+        # β != 0, 1
+        rmul!(y, β)
+        return y
+    else
+        if iszero(β)
+            mul!(y, A, x)
+            !isone(α) && rmul!(y, α)
+            return y
+        else
+            mul!(z, A, x)
+            if isone(α)
+                if isone(β)
+                    y .+= z
+                else
+                    y .= y .* β .+ z
+                end
+            else
+                if isone(β)
+                    y .+= z .* α
+                else
+                    y .= y .* β .+ z .* α
+                end
+            end
+            return y
+        end
+    end
+end
+
+include("transpose.jl") # transposing linear maps
 include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
 include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I
-include("transpose.jl") # transposing linear maps
 include("linearcombination.jl") # defining linear combinations of linear maps
 include("composition.jl") # composition of linear maps
 include("functionmap.jl") # using a function as linear map
diff --git a/src/blockmap.jl b/src/blockmap.jl
index f5caad7a..f6246839 100644
--- a/src/blockmap.jl
+++ b/src/blockmap.jl
@@ -344,24 +344,9 @@ end
 # multiplication with vectors & matrices
 ############
 
-Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) =
-    mul!(y, A, x)
-
-Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::TransposeMap{<:Any,<:BlockMap}, x::AbstractVector) =
-    mul!(y, A, x)
-
-Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) =
-    mul!(y, transpose(A), x)
-
-Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::AdjointMap{<:Any,<:BlockMap}, x::AbstractVector) =
-    mul!(y, A, x)
-
-Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) =
-    mul!(y, adjoint(A), x)
-
 for Atype in (AbstractVector, AbstractMatrix)
     @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, A::BlockMap, x::$Atype,
-                        α::Number=true, β::Number=false)
+                                α::Number=true, β::Number=false)
         require_one_based_indexing(y, x)
         @boundscheck check_dim_mul(y, A, x)
         return _blockmul!(y, A, x, α, β)
@@ -369,7 +354,7 @@ for Atype in (AbstractVector, AbstractMatrix)
 
     for (maptype, transform) in ((:(TransposeMap{<:Any,<:BlockMap}), :transpose), (:(AdjointMap{<:Any,<:BlockMap}), :adjoint))
         @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, wrapA::$maptype, x::$Atype,
-                        α::Number=true, β::Number=false)
+                                α::Number=true, β::Number=false)
             require_one_based_indexing(y, x)
             @boundscheck check_dim_mul(y, wrapA, x)
             return _transblockmul!(y, wrapA.lmap, x, α, β, $transform)
@@ -457,15 +442,6 @@ LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}(
 
 Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) = (eltype(A) == eltype(B) && A.maps == B.maps)
 
-Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) =
-    mul!(y, A, x, true, false)
-
-Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) =
-    mul!(y, transpose(A), x, true, false)
-
-Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) =
-    mul!(y, adjoint(A), x, true, false)
-
 for Atype in (AbstractVector, AbstractMatrix)
     @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, A::BlockDiagonalMap, x::$Atype,
                         α::Number=true, β::Number=false)
@@ -483,3 +459,14 @@ end
     end
     return y
 end
+
+Base.@propagate_inbounds function muladd!(y, A::BlockDiagonalMap, x, α, β, z)
+    require_one_based_indexing(y, x, z)
+    @boundscheck check_dim_mul(y, A, x)
+    @boundscheck size(y) == size(z)
+    maps, yinds, xinds = A.maps, A.rowranges, A.colranges
+    @views @inbounds for i in 1:length(maps)
+        _muladd!(MulStyle(maps[i]), selectdim(y, 1, yinds[i]), maps[i], selectdim(x, 1, xinds[i]), α, β, selectdim(z, 1, yinds[i]))
+    end
+    return y
+end
diff --git a/src/composition.jl b/src/composition.jl
index de644ca1..d7b4ab3c 100644
--- a/src/composition.jl
+++ b/src/composition.jl
@@ -16,6 +16,9 @@ struct CompositeMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
 end
 CompositeMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = CompositeMap{T, As}(maps)
 
+# treat CompositeMaps as FiveArg's because the 5-arg mul! is optimized
+MulStyle(::CompositeMap) = FiveArg()
+
 # basic methods
 Base.size(A::CompositeMap) = (size(A.maps[end], 1), size(A.maps[1], 2))
 Base.isreal(A::CompositeMap) = all(isreal, A.maps) # sufficient but not necessary
@@ -120,17 +123,21 @@ LinearAlgebra.adjoint(A::CompositeMap{T}) where {T}   = CompositeMap{T}(map(adjo
 Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps)
 
 # multiplication with vectors
-function A_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector)
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, A::CompositeMap, x::AbstractVector,
+                    α::Number=true, β::Number=false)
     # no size checking, will be done by individual maps
     N = length(A.maps)
+    A1 = first(A.maps)
     if N==1
-        A_mul_B!(y, A.maps[1], x)
+        mul!(y, A1, x, α, β)
     else
         dest = similar(y, size(A.maps[1], 1))
-        A_mul_B!(dest, A.maps[1], x)
+        _muladd!(MulStyle(A1), dest, A1, x, α, false)
         source = dest
-        if N>2
-            dest = similar(y, size(A.maps[2], 1))
+        if N>2 || (MulStyle(A.maps[2]) === ThreeArg() && !iszero(β))
+            # we need a second intermediate vector if there are more than two maps
+            # or if the second out of two is a ThreeArg and we need to do muladd
+            dest = Array{T}(undef, size(A.maps[2], 1))
         end
         for n in 2:N-1
             try
@@ -142,14 +149,22 @@ function A_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector)
                     rethrow(err)
                 end
             end
-            A_mul_B!(dest, A.maps[n], source)
+            mul!(dest, A.maps[n], source)
             dest, source = source, dest # alternate dest and source
         end
-        A_mul_B!(y, A.maps[N], source)
+        Alast = last(A.maps)
+        if (MulStyle(Alast) === ThreeArg() && !iszero(β))
+            try
+                resize!(dest, size(Alast, 1))
+            catch err
+                if err == ErrorException("cannot resize array with shared data")
+                    dest = Array{T}(undef, size(Alast, 1))
+                else
+                    rethrow(err)
+                end
+            end
+        end
+        _muladd!(MulStyle(Alast), y, Alast, source, true, β, dest)
     end
     return y
 end
-
-At_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
-
-Ac_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)
diff --git a/src/functionmap.jl b/src/functionmap.jl
index 4e1b086f..2421144f 100644
--- a/src/functionmap.jl
+++ b/src/functionmap.jl
@@ -107,15 +107,16 @@ function Base.:(*)(A::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector)
     end
 end
 
-function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
-    (length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch("A_mul_B!"))
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
+    @boundscheck check_dim_mul(y, A, x)
     ismutating(A) ? A.f(y, x) : copyto!(y, A.f(x))
     return y
 end
 
-function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
-    (issymmetric(A) || (isreal(A) && ishermitian(A))) && return A_mul_B!(y, A, x)
-    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("At_mul_B!"))
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, transA::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector)
+    A = transA.lmap
+    issymmetric(A) && return mul!(y, A, x)
+    @boundscheck check_dim_mul(y, transA, x)
     if A.fc !== nothing
         if !isreal(A)
             x = conj(x)
@@ -135,9 +136,10 @@ function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
     end
 end
 
-function Ac_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
-    ishermitian(A) && return A_mul_B!(y, A, x)
-    (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch("Ac_mul_B!"))
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, adjA::AdjointMap{<:Any,<:FunctionMap}, x::AbstractVector)
+    A = adjA.lmap
+    ishermitian(A) && return mul!(y, A, x)
+    @boundscheck check_dim_mul(y, adjA, x)
     if A.fc !== nothing
         ismutating(A) ? A.fc(y, x) : copyto!(y, A.fc(x))
         return y
diff --git a/src/kronecker.jl b/src/kronecker.jl
index 80fd502a..8db9b252 100644
--- a/src/kronecker.jl
+++ b/src/kronecker.jl
@@ -109,9 +109,9 @@ Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps
     temp2 = similar(y, nb)
     @views @inbounds for i in 1:ma
         v[i] = one(T)
-        A_mul_B!(temp1, At, v)
-        A_mul_B!(temp2, X, temp1)
-        A_mul_B!(y[((i-1)*mb+1):i*mb], B, temp2)
+        mul!(temp1, At, v)
+        mul!(temp2, X, temp1)
+        mul!(y[((i-1)*mb+1):i*mb], B, temp2)
         v[i] = zero(T)
     end
     return y
@@ -131,7 +131,7 @@ end
 # multiplication with vectors
 #################
 
-Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T}
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T}
     require_one_based_indexing(y)
     @boundscheck check_dim_mul(y, L, x)
     A, B = L.maps
@@ -139,7 +139,7 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,
     _kronmul!(y, B, X, transpose(A), T)
     return y
 end
-Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T}
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T}
     require_one_based_indexing(y)
     @boundscheck check_dim_mul(y, L, x)
     A = first(L.maps)
@@ -150,35 +150,29 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}
 end
 # mixed-product rule, prefer the right if possible
 # (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ)
-Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector)
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector)
     B, A = L.maps
     if length(A.maps) == length(B.maps) && all(M -> check_dim_mul(M[1], M[2]), zip(A.maps, B.maps))
-        A_mul_B!(y, kron(map(*, A.maps, B.maps)...), x)
+        mul!(y, kron(map(*, A.maps, B.maps)...), x)
     else
-        A_mul_B!(y, LinearMap(A)*B, x)
+        mul!(y, LinearMap(A)*B, x)
     end
 end
 # mixed-product rule, prefer the right if possible
 # (A₁ ⊗ B₁)*(A₂⊗B₂)*...*(Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ)
-Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T}
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T}
     As = map(AB -> AB.maps[1], L.maps)
     Bs = map(AB -> AB.maps[2], L.maps)
     As1, As2 = Base.front(As), Base.tail(As)
     Bs1, Bs2 = Base.front(Bs), Base.tail(Bs)
     apply = all(A -> check_dim_mul(A...), zip(As1, As2)) && all(A -> check_dim_mul(A...), zip(Bs1, Bs2))
     if apply
-        A_mul_B!(y, kron(prod(As), prod(Bs)), x)
+        mul!(y, kron(prod(As), prod(Bs)), x)
     else
-        A_mul_B!(y, CompositeMap{T}(map(LinearMap, L.maps)), x)
+        mul!(y, CompositeMap{T}(map(LinearMap, L.maps)), x)
     end
 end
 
-Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) =
-    A_mul_B!(y, transpose(A), x)
-
-Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) =
-    A_mul_B!(y, adjoint(A), x)
-
 ###############
 # KroneckerSumMap
 ###############
@@ -249,6 +243,8 @@ where `A` can be a square `AbstractMatrix` or a `LinearMap`.
 
 Base.:(^)(A::MapOrMatrix, ::KronSumPower{p}) where {p} = kronsum(ntuple(n -> convert_to_lmaps_(A), Val(p))...)
 
+MulStyle(A::KronSumPower) = MulStyle(Matrix{Float64}(undef, 0, 0), A.maps...)
+
 Base.size(A::KroneckerSumMap, i) = prod(size.(A.maps, i))
 Base.size(A::KroneckerSumMap) = (size(A, 1), size(A, 2))
 
@@ -261,20 +257,30 @@ LinearAlgebra.transpose(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(ma
 
 Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = (eltype(A) == eltype(B) && A.maps == B.maps)
 
-Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector)
+Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector,
+                    α::Number=true, β::Number=false)
     @boundscheck check_dim_mul(y, L, x)
     A, B = L.maps
     ma, na = size(A)
     mb, nb = size(B)
     X = reshape(x, (nb, na))
     Y = reshape(y, (nb, na))
-    mul!(Y, X, convert(AbstractMatrix, transpose(A)))
-    mul!(Y, B, X, true, true)
+    _muladd!(MulStyle(X), Y, X, convert(AbstractMatrix, transpose(A)), true, β)
+    _muladd!(MulStyle(B), Y, B, X, α, true)
     return y
 end
 
-Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) =
-    A_mul_B!(y, transpose(A), x)
-
-Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) =
-    A_mul_B!(y, adjoint(A), x)
+Base.@propagate_inbounds function muladd!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector,
+                    α, β, z::AbstractVector)
+    @boundscheck check_dim_mul(y, L, x)
+    @boundscheck size(y) == size(z)
+    A, B = L.maps
+    ma, na = size(A)
+    mb, nb = size(B)
+    X = reshape(x, (nb, na))
+    Y = reshape(y, (nb, na))
+    Z = reshape(z, (nb, na))
+    _muladd!(MulStyle(X), Y, X, convert(AbstractMatrix, transpose(A)), true, β, Z)
+    _muladd!(MulStyle(B), Y, B, X, α, true, Z)
+    return y
+end
diff --git a/src/linearcombination.jl b/src/linearcombination.jl
index 7f556918..1732c4a4 100644
--- a/src/linearcombination.jl
+++ b/src/linearcombination.jl
@@ -75,9 +75,21 @@ for Atype in (AbstractVector, AbstractMatrix)
             rmul!(y, β)
             return y
         else
-            mul!(y, first(A.maps), x, α, β)
-            return _mul!(MulStyle(A), y, A, x, α)
+            A1 = first(A.maps)
+            if MulStyle(A1) === ThreeArg() && !iszero(β)
+                # if we need an intermediate vector, allocate here and reuse in
+                # LinearCombination multiplication
+                z = similar(y)
+                muladd!(y, A1, x, α, β, z)
+                __mul!(y, Base.tail(A.maps), x, α, z)
+            else # MulStyle(A1) === FiveArg()
+                # this is allocation-free
+                mul!(y, A1, x, α, β)
+                # let _mul! decide whether an intermediate vector needs to be allocated
+                _mul!(MulStyle(A), y, A, x, α)
+            end
         end
+        return y
     end
 end
 
@@ -91,21 +103,4 @@ end
 @inline __mul!(y, As::Tuple{Vararg{LinearMap}}, x, α, z) =
     __mul!(__mul!(y, first(As), x, α, z), Base.tail(As), x, α, z)
 @inline __mul!(y, A::Tuple{LinearMap}, x, α, z) = __mul!(y, first(A), x, α, z)
-@inline __mul!(y, A::LinearMap, x, α, z) = muladd!(MulStyle(A), y, A, x, α, z)
-
-@inline muladd!(::FiveArg, y, A, x, α, _) = mul!(y, A, x, α, true)
-@inline function muladd!(::ThreeArg, y, A, x, α, z)
-    mul!(z, A, x)
-    if isone(α)
-        y .+= z
-    else
-        y .+= z .* α
-    end
-    return y
-end
-
-A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, A, x)
-
-At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, transpose(A), x)
-
-Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = mul!(y, adjoint(A), x)
+@inline __mul!(y, A::LinearMap, x, α, z) = _muladd!(MulStyle(A), y, A, x, α, true, z)
diff --git a/src/transpose.jl b/src/transpose.jl
index 60caf3c2..664f3c3b 100644
--- a/src/transpose.jl
+++ b/src/transpose.jl
@@ -32,18 +32,18 @@ Base.:(==)(A::LinearMap, B::TransposeMap)    = issymmetric(A) && B.lmap == A
 Base.:(==)(A::LinearMap, B::AdjointMap)      = ishermitian(A) && B.lmap == A
 
 # multiplication with vector
-A_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) =
-    issymmetric(A.lmap) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x)
+Base.@propagate_inbounds LinearAlgebra.mul!(y::AbstractVector, A::TransposeMap, x::AbstractVector,
+                α::Number=true, β::Number=false) =
+    issymmetric(A.lmap) ? _muladd!(MulStyle(A), y, A.lmap, x, α, β) : _muladd!(MulStyle(A), y, A, x, α, β)
 
-At_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
+Base.@propagate_inbounds LinearAlgebra.mul!(y::AbstractVector, A::AdjointMap{<:Any,<:TransposeMap}, x::AbstractVector,
+                α::Number=true, β::Number=false) =
+    isreal(A.lmap) ? _muladd!(MulStyle(A), y, A.lmap.lmap, x, α, β) : (_muladd!(MulStyle(A), y, A.lmap.lmap, conj(x), conj(α), conj(β)); conj!(y))
 
-Ac_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) =
-    isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y))
+Base.@propagate_inbounds LinearAlgebra.mul!(y::AbstractVector, A::AdjointMap, x::AbstractVector,
+                α::Number=true, β::Number=false) =
+    ishermitian(A.lmap) ? _muladd!(MulStyle(A), y, A.lmap, x, α, β) : _muladd!(MulStyle(A), y, A, x, α, β)
 
-A_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) =
-    ishermitian(A.lmap) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x)
-
-At_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) =
-    isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y))
-
-Ac_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
+Base.@propagate_inbounds LinearAlgebra.mul!(y::AbstractVector, A::TransposeMap{<:Any,<:AdjointMap}, x::AbstractVector,
+                α::Number=true, β::Number=false) =
+    isreal(A.lmap) ? _muladd!(MulStyle(A), y, A.lmap.lmap, x, α, β) : (_muladd!(MulStyle(A), y, A.lmap.lmap, conj(x), conj(α), conj(β)); conj!(y))
diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl
index 537467c2..bb3548a9 100644
--- a/src/uniformscalingmap.jl
+++ b/src/uniformscalingmap.jl
@@ -76,18 +76,12 @@ function _scaling!(y, λ::Number, x, α::Number, β::Number)
             y .+= λ .* x .* α
             return y
         else # β != 0, 1
-            isone(λ) && (y .= y .* β .+ x .* α; return y)
             y .= y .* β .+ λ .* x .* α
             return y
         end # β-cases
     end # α-cases
 end # function _scaling!
 
-A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) = mul!(y, A, x)
-At_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
-Ac_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)
-
-
 # combine LinearMap and UniformScaling objects in linear combinations
 Base.:(+)(A₁::LinearMap, A₂::UniformScaling) = A₁ + UniformScalingMap(A₂.λ, size(A₁, 1))
 Base.:(+)(A₁::UniformScaling, A₂::LinearMap) = UniformScalingMap(A₁.λ, size(A₂, 1)) + A₂
diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl
index 49ee410d..2f6e089d 100644
--- a/src/wrappedmap.jl
+++ b/src/wrappedmap.jl
@@ -36,33 +36,6 @@ LinearAlgebra.issymmetric(A::WrappedMap) = A._issymmetric
 LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian
 LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef
 
-# multiplication with vectors & matrices
-A_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)
-Base.:(*)(A::WrappedMap, x::AbstractVector) = *(A.lmap, x)
-
-At_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) =
-    (issymmetric(A) || (isreal(A) && ishermitian(A))) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x)
-
-Ac_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) =
-    ishermitian(A) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x)
-
-if VERSION ≥ v"1.3.0-alpha.115"
-    for Atype in (AbstractVector, AbstractMatrix)
-        @eval Base.@propagate_inbounds LinearAlgebra.mul!(y::$Atype, A::WrappedMap, x::$Atype,
-                        α::Number=true, β::Number=false) =
-            mul!(y, A.lmap, x, α, β)
-    end
-else
-# This is somewhat suboptimal, because the absence of 5-arg mul! for MatrixMaps
-# doesn't allow to define a 5-arg mul! for WrappedMaps which do have a 5-arg mul!
-# I'd assume, however, that 5-arg mul! becomes standard in Julia v≥1.3 anyway
-# the idea is to let the fallback handle 5-arg calls
-    for Atype in (AbstractVector, AbstractMatrix)
-        @eval Base.@propagate_inbounds LinearAlgebra.mul!(Y::$Atype, A::WrappedMap, X::$Atype) =
-            mul!(Y, A.lmap, X)
-    end
-end # VERSION
-
 # combine LinearMap and Matrix objects: linear combinations and map composition
 Base.:(+)(A₁::LinearMap, A₂::AbstractMatrix) = +(A₁, WrappedMap(A₂))
 Base.:(+)(A₁::AbstractMatrix, A₂::LinearMap) = +(WrappedMap(A₁), A₂)
@@ -71,3 +44,47 @@ Base.:(-)(A₁::AbstractMatrix, A₂::LinearMap) = -(WrappedMap(A₁), A₂)
 
 Base.:(*)(A₁::LinearMap, A₂::AbstractMatrix) = *(A₁, WrappedMap(A₂))
 Base.:(*)(A₁::AbstractMatrix, A₂::LinearMap) = *(WrappedMap(A₁), A₂)
+
+# multiplication with vector/matrix
+for Atype in (AbstractVector, AbstractMatrix)
+    @eval Base.@propagate_inbounds LinearAlgebra.mul!(y::$Atype, A::WrappedMap, x::$Atype,
+                            α::Number, β::Number=false) =
+        _muladd!(MulStyle(A), y, A.lmap, x, α, β)
+    @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, A::TransposeMap{<:Any,<:WrappedMap}, x::$Atype,
+                            α::Number, β::Number=false)
+        if (issymmetric(A) || (isreal(A) && ishermitian(A)))
+            _muladd!(MulStyle(A), y, A.lmap, x, α, β)
+        else
+            _muladd!(MulStyle(A), y, transpose(A.lmap.lmap), x, α, β)
+        end
+        return y
+    end
+    @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, adjA::AdjointMap{<:Any,<:WrappedMap}, x::$Atype,
+                            α::Number, β::Number=false)
+        if ishermitian(adjA)
+            _muladd!(MulStyle(adjA), y, adjA.lmap, x, α, β)
+        else
+            _muladd!(MulStyle(adjA), y, adjoint(adjA.lmap.lmap), x, α, β)
+        end
+        return y
+    end
+end
+
+Base.@propagate_inbounds muladd!(y, A::WrappedMap, x, α, β, z) =
+    _muladd!(MulStyle(A), y, A.lmap, x, α, β, z)
+Base.@propagate_inbounds function muladd!(y, A::TransposeMap{<:Any,<:WrappedMap}, x, α, β, z)
+    if (issymmetric(A) || (isreal(A) && ishermitian(A)))
+        _muladd!(MulStyle(A), y, A.lmap, x, α, β, z)
+    else
+        _muladd!(MulStyle(A), y, transpose(A.lmap.lmap), x, α, β, z)
+    end
+    return y
+end
+Base.@propagate_inbounds function muladd!(y, adjA::AdjointMap{<:Any,<:WrappedMap}, x, α, β, z)
+    if ishermitian(adjA)
+        _muladd!(MulStyle(adjA), y, adjA.lmap, x, α, β, z)
+    else
+        _muladd!(MulStyle(adjA), y, adjoint(adjA.lmap.lmap), x, α, β, z)
+    end
+    return y
+end
diff --git a/test/blockmap.jl b/test/blockmap.jl
index 49103056..33c09b01 100644
--- a/test/blockmap.jl
+++ b/test/blockmap.jl
@@ -1,4 +1,4 @@
-using Test, LinearMaps, LinearAlgebra, SparseArrays
+using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools
 
 @testset "block maps" begin
     @testset "hcat" begin
@@ -203,4 +203,28 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
             end
         end
     end
+
+    @testset "function block map" begin
+        CS! = LinearMap{ComplexF64}(cumsum!,
+                                    (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y); reverse!(y)), 10;
+                                    ismutating=true)
+        A = rand(ComplexF64, 10, 10)
+        B = rand(10, 10)
+        L = [CS! LinearMap(A); LinearMap(B) CS!]
+        M = [Matrix(CS!) A; B Matrix(CS!)]
+        u = rand(ComplexF64, 20)
+        v = rand(ComplexF64, 20)
+        for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64))
+            for transform in (identity, adjoint)
+                @test mul!(copy(v), transform(L), u, α, β) ≈ transform(M)*u*α + v*β
+                @test mul!(copy(v), transform(LinearMap(L)), u, α, β) ≈ transform(M)*u*α + v*β
+                @test mul!(copy(v), LinearMap(transform(L)), u, α, β) ≈ transform(M)*u*α + v*β
+                # bmap = @benchmarkable mul!($(copy(v)), $(transform(L)), $u, $α, $β)
+                # bmat = @benchmarkable mul!($(copy(v)), $(transform(Lmat)), $u, $α, $β)
+                # @show run(bmap, samples=3).allocs
+                # @show run(bmat, samples=3).allocs
+                # @test run(blm, samples=3).allocs <= 1
+            end
+        end
+    end
 end
diff --git a/test/composition.jl b/test/composition.jl
index 4ca162fc..d20c7e8d 100644
--- a/test/composition.jl
+++ b/test/composition.jl
@@ -53,6 +53,24 @@ using Test, LinearMaps, LinearAlgebra
     mul!(w, L, v)
     @test w ≈ LF * v
 
+    # test new type
+    F = SimpleFunctionMap(cumsum, 10)
+    FC = SimpleComplexFunctionMap(cumsum, 10)
+    @test @inferred ndims(F) == 2
+    @test @inferred size(F, 1) == 10
+    @test @inferred length(F) == 100
+    @test @inferred !issymmetric(F)
+    @test @inferred !ishermitian(F)
+    @test @inferred !ishermitian(FC)
+    @test @inferred !isposdef(F)
+    w = similar(v)
+    mul!(w, F, v)
+    @test w == F * v
+    # @test_throws StackOverflowError F' * v
+    # @test_throws StackOverflowError transpose(F) * v
+    # @test_throws StackOverflowError mul!(w, adjoint(F), v)
+    # @test_throws StackOverflowError mul!(w, transpose(F), v)
+
     # test composition of several maps with shared data #31
     global sizes = ( (5, 2), (3, 3), (3, 2), (2, 2), (9, 2), (7, 1) )
     N = length(sizes) - 1
diff --git a/test/functionmap.jl b/test/functionmap.jl
index 0eec87c5..d2a5255e 100644
--- a/test/functionmap.jl
+++ b/test/functionmap.jl
@@ -56,9 +56,10 @@ using Test, LinearMaps, LinearAlgebra
     @test @inferred mul!(similar(v), CS!, v) == cv
     @test_throws ErrorException CS!'v
     @test_throws ErrorException adjoint(CS!) * v
-    @test_throws ErrorException mul!(similar(v), CS!', v)
-    @test_throws ErrorException mul!(similar(v), transpose(CS!), v)
-    CS! = LinearMap{ComplexF64}(cumsum!, (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y)), 10; ismutating=true)
+    CS! = LinearMap{ComplexF64}(cumsum!,
+                                (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y); reverse!(y)), 10;
+                                ismutating=true)
+    M = Matrix(CS!)
     @inferred adjoint(CS!)
     @test @inferred LinearMaps.ismutating(CS!)
     @test @inferred CS! * v == cv
@@ -67,12 +68,27 @@ using Test, LinearMaps, LinearAlgebra
     @test @inferred CS' * v == reverse!(cumsum(reverse(v)))
     @test @inferred mul!(similar(v), transpose(CS), v) == reverse!(cumsum(reverse(v)))
     @test @inferred mul!(similar(v), adjoint(CS), v) == reverse!(cumsum(reverse(v)))
+    u = rand(ComplexF64, 10)
+    v = rand(ComplexF64, 10)
+    for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64))
+        for transform in (identity, adjoint, transpose)
+            @test mul!(copy(v), transform(CS!), u, α, β) ≈ transform(M)*u*α + v*β
+            @test mul!(copy(v), transform(LinearMap(CS!)), u, α, β) ≈ transform(M)*u*α + v*β
+            @test mul!(copy(v), LinearMap(transform(CS!)), u, α, β) ≈ transform(M)*u*α + v*β
+            if transform != transpose
+                bm = @benchmarkable mul!($(copy(v)), $(transform(CS!)), $u, $α, $β)
+                @test run(bm, samples=3).allocs <= 1
+            end
+        end
+    end
 
     # Test fallback methods:
     L = @inferred LinearMap(x -> x, x -> x, 10)
     v = randn(10)
     @test @inferred (2 * L)' * v ≈ 2 * v
     @test @inferred transpose(2 * L) * v ≈ 2 * v
+    @test @inferred (L*2)' * v ≈ 2 * v
+    @test @inferred transpose(L*2) * v ≈ 2 * v
     L = @inferred LinearMap{ComplexF64}(x -> x, x -> x, 10)
     v = rand(ComplexF64, 10)
     w = similar(v)
diff --git a/test/linearcombination.jl b/test/linearcombination.jl
index 31b40b55..37588b9d 100644
--- a/test/linearcombination.jl
+++ b/test/linearcombination.jl
@@ -2,7 +2,7 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
 
 @testset "linear combinations" begin
     CS! = LinearMap{ComplexF64}(cumsum!,
-                                (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y)), 10;
+                                (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y); reverse!(y)), 10;
                                 ismutating=true)
     v = rand(ComplexF64, 10)
     u = similar(v)
@@ -14,7 +14,9 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
     b = @benchmarkable mul!($u, $L, $v)
     @test run(b, samples=5).allocs <= 1
     for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64))
-        @test mul!(copy(u), L, v, α, β) ≈ Matrix(L)*v*α + u*β
+        for transform in (identity, adjoint, transpose)
+            @test mul!(copy(u), transform(L), v, α, β) ≈ transform(Matrix(L))*v*α + u*β
+        end
     end
     V = rand(ComplexF64, 10, 3)
     U = similar(V)
@@ -44,17 +46,35 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
         b = @benchmarkable mul!($w, $LC, $v)
         @test run(b, samples=3).allocs == 0
         for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64))
-            if testallocs
-                b = @benchmarkable mul!($w, $LC, $v, $α, $β)
-                @test run(b, samples=3).allocs == 0
-                b = @benchmarkable mul!($w, $(I + LC), $v, $α, $β)
-                @test run(b, samples=3).allocs == 0
-                b = @benchmarkable mul!($w, $(LC + I), $v, $α, $β)
-                @test run(b, samples=3).allocs == 0
+            for λ in (true, rand(ComplexF64)), sz in (10, (10, 10))
+                v = rand(ComplexF64, sz)
+                w = similar(v)
+                if testallocs
+                    blmap = @benchmarkable mul!($w, $LC, $v, $α, $β)
+                    bmat = @benchmarkable begin
+                        mul!($w, $A, $v, $α, $β)
+                        mul!($w, $B, $v, $α, true)
+                    end
+                    lallocs = run(blmap, samples=3).allocs
+                    mallocs = run(bmat, samples=3).allocs
+                    @test lallocs <= mallocs
+                    @test lallocs == 0
+
+                    blmap = @benchmarkable mul!($w, $(LC + λ*I), $v, $α, $β)
+                    bmat = @benchmarkable begin
+                        mul!($w, $A, $v, $α, $β)
+                        mul!($w, $B, $v, $α, true)
+                        mul!($w, $(λ*I), $v, $α, true)
+                    end
+                    lallocs = run(blmap, samples=3).allocs
+                    mallocs = run(bmat, samples=3).allocs
+                    @test lallocs <= mallocs
+                    @test lallocs == 0
+                end
+                y = rand(ComplexF64, sz)
+                @test mul!(copy(y), LC, v, α, β) ≈ (A + B)*v*α + y*β
+                @test mul!(copy(y), LC+λ*I, v, α, β) ≈ (A + B + λ*I)*v*α + y*β
             end
-            y = rand(ComplexF64, size(v))
-            @test mul!(copy(y), LC, v, α, β) ≈ Matrix(LC)*v*α + y*β
-            @test mul!(copy(y), LC+I, v, α, β) ≈ Matrix(LC + I)*v*α + y*β
         end
     end
     # @test_throws ErrorException LinearMaps.LinearCombination{ComplexF64}((M, N), (1, 2, 3))
diff --git a/test/runtests.jl b/test/runtests.jl
index b5444f0c..42e9bbf0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,7 +3,7 @@ import LinearMaps: FiveArg, ThreeArg
 
 const matrixstyle = VERSION ≥ v"1.3.0-alpha.115" ? FiveArg() : ThreeArg()
 
-const testallocs = true
+const testallocs = VERSION ≥ v"1.4.0-"
 
 include("linearmaps.jl")
 
diff --git a/test/transpose.jl b/test/transpose.jl
index 705af876..59c65dbb 100644
--- a/test/transpose.jl
+++ b/test/transpose.jl
@@ -66,6 +66,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
     x = rand(ComplexF64, 10)
     for transform1 in (adjoint, transpose), transform2 in (adjoint, transpose)
         @test transform2(LinearMap(transform1(CS))) * x ≈ transform2(transform1(M))*x
+        @test transform2(transform1(CS)) * x ≈ transform2(transform1(M))*x
     end
 
     id = @inferred LinearMap(identity, identity, 10; issymmetric=true, ishermitian=true, isposdef=true)
diff --git a/test/uniformscalingmap.jl b/test/uniformscalingmap.jl
index f83715d4..f2229836 100644
--- a/test/uniformscalingmap.jl
+++ b/test/uniformscalingmap.jl
@@ -25,10 +25,10 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
     @test (3 * I - 2 * M') * v == -2 * A'v + 3v
     @test transpose(LinearMap(2 * M' + 3 * I)) * v ≈ transpose(2 * A' + 3 * I) * v
     @test LinearMap(2 * M' + 0I)' * v ≈ (2 * A')' * v
-    for λ in (0, 1, rand()), α in (0, 1, rand()), β in (0, 1, rand())
+    for λ in (0, 1, rand()), α in (0, 1, rand()), β in (0, 1, rand()), sz in (10, (10,5))
         Λ = @inferred LinearMap(λ*I, 10)
-        x = rand(10)
-        y = rand(10)
+        x = rand(Float64, sz)
+        y = rand(Float64, sz)
         if testallocs
             b = @benchmarkable mul!($y, $Λ, $x, $α, $β)
             @test run(b, samples=3).allocs == 0
diff --git a/test/wrappedmap.jl b/test/wrappedmap.jl
index 57171b25..7151652b 100644
--- a/test/wrappedmap.jl
+++ b/test/wrappedmap.jl
@@ -1,4 +1,4 @@
-using Test, LinearMaps, LinearAlgebra
+using Test, LinearMaps, LinearAlgebra, BenchmarkTools
 
 @testset "wrapped maps" begin
     A = rand(10, 20)
@@ -14,4 +14,35 @@ using Test, LinearMaps, LinearAlgebra
     @test @inferred !issymmetric(MB)
     @test @inferred isposdef(MA)
     @test @inferred isposdef(MB)
+    B = rand(ComplexF64, 20, 20)
+    v = rand(ComplexF64, 20)
+    u = rand(ComplexF64, 20)
+    LB = LinearMap(B)
+    for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64))
+        for transform in (identity, adjoint, transpose)
+            @test mul!(copy(u), transform(LB), v, α, β) ≈ transform(B)*v*α + u*β
+            @test mul!(copy(u), transform(LinearMap(LB)), v, α, β) ≈ transform(B)*v*α + u*β
+            @test mul!(copy(u), transform(LinearMap(LinearMap(LB))), v, α, β) ≈ transform(B)*v*α + u*β
+            @test mul!(copy(u), LinearMap(transform(LinearMap(LB))), v, α, β) ≈ transform(B)*v*α + u*β
+            if testallocs
+                blmap = @benchmarkable mul!($(copy(u)), $(transform(LB)), $v, $α, $β)
+                @test run(blmap, samples=3).allocs == 0
+            end
+        end
+    end
+
+    CS! = LinearMap{ComplexF64}(cumsum!,
+                                (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y); reverse!(y)), 10;
+                                ismutating=true)
+    v = rand(ComplexF64, 10)
+    u = rand(ComplexF64, 10)
+    M = Matrix(CS!)
+    for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64))
+        for transform in (identity, adjoint, transpose)
+            @test mul!(copy(u), transform(CS!), v, α, β) ≈ transform(M)*v*α + u*β
+            @test mul!(copy(u), transform(LinearMap(CS!)), v, α, β) ≈ transform(M)*v*α + u*β
+            @test mul!(copy(u), transform(LinearMap(LinearMap(CS!))), v, α, β) ≈ transform(M)*v*α + u*β
+            @test mul!(copy(u), LinearMap(transform(LinearMap(CS!))), v, α, β) ≈ transform(M)*v*α + u*β
+        end
+    end
 end