Simplify multiplication code (#51)

* simplify multiplication code, more consistent code style

* remove eachcol again

Author: dkarrasch

src/LinearMaps.jl

@@ -50,6 +50,10 @@ function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap{T}, X::AbstractMatri
     @inbounds @views for i = 1:size(X, 2)
         mul!(Y[:, i], A, X[:, i], α, β)
     end
+    # starting from Julia v1.1, we could use the `eachcol` iterator
+    # for (Xi, Yi) in zip(eachcol(X), eachcol(Y))
+    #     mul!(Yi, A, Xi, α, β)
+    # end
     return Y
 end
 

src/composition.jl

@@ -131,63 +131,7 @@ function A_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector)
     end
     return y
 end
-function At_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector)
-    # no size checking, will be done by individual maps
-    N = length(A.maps)
-    if N == 1
-        At_mul_B!(y, A.maps[1], x)
-    else
-        T = eltype(y)
-        dest = Array{T}(undef, size(A.maps[N], 2))
-        At_mul_B!(dest, A.maps[N], x)
-        source = dest
-        if N>2
-            dest = Array{T}(undef, size(A.maps[N-1], 2))
-        end
-        for n = N-1:-1:2
-            try
-                resize!(dest, size(A.maps[n], 2))
-            catch err
-                if err == ErrorException("cannot resize array with shared data")
-                    dest = Array{T}(undef, size(A.maps[n], 2))
-                else
-                    rethrow(err)
-                end
-            end
-            At_mul_B!(dest, A.maps[n], source)
-            dest, source = source, dest # alternate dest and source
-        end
-        At_mul_B!(y, A.maps[1], source)
-    end
-    return y
-end
-function Ac_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector)
-    # no size checking, will be done by individual maps
-    N = length(A.maps)
-    if N == 1
-        Ac_mul_B!(y, A.maps[1], x)
-    else
-        T = eltype(y)
-        dest = Array{T}(undef, size(A.maps[N], 2))
-        Ac_mul_B!(dest, A.maps[N], x)
-        source = dest
-        if N>2
-            dest = Array{T}(undef, size(A.maps[N-1], 2))
-        end
-        for n = N-1:-1:2
-            try
-                resize!(dest, size(A.maps[n], 2))
-            catch err
-                if err == ErrorException("cannot resize array with shared data")
-                    dest = Array{T}(undef, size(A.maps[n], 2))
-                else
-                    rethrow(err)
-                end
-            end
-            Ac_mul_B!(dest, A.maps[n], source)
-            dest, source = source, dest # alternate dest and source
-        end
-        Ac_mul_B!(y, A.maps[1], source)
-    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)

src/functionmap.jl

@@ -38,11 +38,6 @@ ismutating(A::FunctionMap) = A._ismutating
 _ismutating(f) = first(methods(f)).nargs == 3
 
 # multiplication with vector
-function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
-    (length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch())
-    ismutating(A) ? A.f(y, x) : copyto!(y, A.f(x))
-    return y
-end
 function Base.:(*)(A::FunctionMap, x::AbstractVector)
     length(x) == A.N || throw(DimensionMismatch())
     if ismutating(A)
@@ -54,6 +49,12 @@ function Base.:(*)(A::FunctionMap, x::AbstractVector)
     return y
 end
 
+function A_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
+    (length(x) == A.N && length(y) == A.M) || throw(DimensionMismatch())
+    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) && return A_mul_B!(y, A, x)
     (length(x) == A.M && length(y) == A.N) || throw(DimensionMismatch())

src/linearcombination.jl

@@ -59,29 +59,7 @@ function A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
     end
     return y
 end
-function At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
-    # no size checking, will be done by individual maps
-    At_mul_B!(y, A.maps[1], x)
-    l = length(A.maps)
-    if l>1
-        z = similar(y)
-        for n = 2:l
-            At_mul_B!(z, A.maps[n], x)
-            y .+= z
-        end
-    end
-    return y
-end
-function Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
-    # no size checking, will be done by individual maps
-    Ac_mul_B!(y, A.maps[1], x)
-    l = length(A.maps)
-    if l>1
-        z = similar(y)
-        for n=2:l
-            Ac_mul_B!(z, A.maps[n], x)
-            y .+= z
-        end
-    end
-    return y
-end
+
+At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
+
+Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)

src/wrappedmap.jl

@@ -40,4 +40,4 @@ Base.:(-)(A1::LinearMap, A2::AbstractMatrix) = -(A1, WrappedMap(A2))
 Base.:(-)(A1::AbstractMatrix, A2::LinearMap) = -(WrappedMap(A1), A2)
 
 Base.:(*)(A1::LinearMap, A2::AbstractMatrix) = *(A1, WrappedMap(A2))
-Base.:(*)(A1::AbstractMatrix, A2::LinearMap) = *(WrappedMap(A1) ,A2)
+Base.:(*)(A1::AbstractMatrix, A2::LinearMap) = *(WrappedMap(A1), A2)