fix lincomb mul with matrices

Author: dkarrasch

src/linearcombination.jl

@@ -95,8 +95,7 @@ end
 
 @inline muladd!(::FiveArg, y, A, x, α, _) = mul!(y, A, x, α, true)
 @inline function muladd!(::ThreeArg, y, A, x, α, z)
-    # TODO: replace by mul!(z, A, x)
-    A_mul_B!(z, A, x)
+    mul!(z, A, x)
     if isone(α)
         y .+= z
     else

test/linearcombination.jl

@@ -16,6 +16,10 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools
     for α in (false, true, rand(ComplexF64)), β in (false, true, rand(ComplexF64))
         @test mul!(copy(u), L, v, α, β) ≈ Matrix(L)*v*α + u*β
     end
+    V = rand(ComplexF64, 10, 3)
+    U = similar(V)
+    @test mul!(U, L, V) ≈ n*cumsum(V, dims=1)
+    @test mul!(U, LinearMap(L), V) ≈ n*cumsum(V, dims=1)
 
     A = 2 * rand(ComplexF64, (10, 10)) .- 1
     B = rand(ComplexF64, size(A)...)