WIP: review before version 3

src/LinearMaps.jl

@@ -16,8 +16,8 @@ end
 
 abstract type LinearMap{T} end
 
-const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}}
-const RealOrComplex = Union{Real,Complex}
+const MapOrMatrix{T} = Union{LinearMap{T}, AbstractMatrix{T}}
+const RealOrComplex = Union{Real, Complex}
 
 Base.eltype(::LinearMap{T}) where {T} = T
 
@@ -45,7 +45,8 @@ LinearAlgebra.ishermitian(::LinearMap) = false # default assumptions
 LinearAlgebra.isposdef(::LinearMap) = false # default assumptions
 
 Base.ndims(::LinearMap) = 2
-Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
+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]
 
 """
@@ -115,8 +116,9 @@ end
     mul!(Y::AbstractVecOrMat, A::LinearMap, B::AbstractVector) -> Y
     mul!(Y::AbstractMatrix, A::LinearMap, B::AbstractMatrix) -> Y
 
-Calculates the action of the linear map `A` on the vector or matrix `B` and stores the result in `Y`,
-overwriting the existing value of `Y`. Note that `Y` must not be aliased with either `A` or `B`.
+Calculates the action of the linear map `A` on the vector or matrix `B` and stores the
+result in `Y`, overwriting the existing value of `Y`. Note that `Y` must not be aliased
+with either `A` or `B`.
 
 ## Examples
 ```jldoctest; setup=(using LinearAlgebra, LinearMaps)
@@ -144,26 +146,26 @@ end
     mul!(C::AbstractVecOrMat, A::LinearMap, B::AbstractVector, α, β) -> C
     mul!(C::AbstractMatrix, A::LinearMap, B::AbstractMatrix, α, β) -> C
 
-Combined inplace multiply-add ``A B α + C β``. The result is stored in `C` by overwriting it.
-Note that `C` must not be aliased with either `A` or `B`.
+Combined inplace multiply-add ``A B α + C β``. The result is stored in `C` by overwriting
+it. Note that `C` must not be aliased with either `A` or `B`.
 
 ## Examples
 ```jldoctest; setup=(using LinearAlgebra, LinearMaps)
 julia> A=LinearMap([1.0 2.0; 3.0 4.0]); B=[1.0, 1.0]; C=[1.0, 3.0];
-  
+
 julia> mul!(C, A, B, 100.0, 10.0) === C
 true
-  
+
 julia> C
 2-element Array{Float64,1}:
  310.0
  730.0
 
 julia> A=LinearMap([1.0 2.0; 3.0 4.0]); B=[1.0 1.0; 1.0 1.0]; C=[1.0 2.0; 3.0 4.0];
-  
+
 julia> mul!(C, A, B, 100.0, 10.0) === C
 true
-  
+
 julia> C
 2×2 Array{Float64,2}:
  310.0  320.0
@@ -257,15 +259,17 @@ with the purpose of redefining its properties via the keyword arguments `kwargs`
 a `UniformScaling` object `J` with specified (square) dimension `M`; or
 from a function or callable object `f`. In the latter case, one also needs to specify
 the size of the equivalent matrix representation `(M, N)`, i.e., for functions `f` acting
-on length `N` vectors and producing length `M` vectors (with default value `N=M`). Preferably,
-also the `eltype` `T` of the corresponding matrix representation needs to be specified, i.e.
-whether the action of `f` on a vector will be similar to, e.g., multiplying by numbers of type `T`.
-If not specified, the devault value `T=Float64` will be assumed. Optionally, a corresponding
-function `fc` can be specified that implements the adjoint (=transpose in the real case) of `f`.
+on length `N` vectors and producing length `M` vectors (with default value `N=M`).
+Preferably, also the `eltype` `T` of the corresponding matrix representation needs to be
+specified, i.e. whether the action of `f` on a vector will be similar to, e.g., multiplying
+by numbers of type `T`. If not specified, the devault value `T=Float64` will be assumed.
+Optionally, a corresponding function `fc` can be specified that implements the adjoint
+(=transpose in the real case) of `f`.
 
 The keyword arguments and their default values for the function-based constructor are:
 *   `issymmetric::Bool = false` : whether `A` or `f` acts as a symmetric matrix
-*   `ishermitian::Bool = issymmetric & T<:Real` : whether `A` or `f` acts as a Hermitian matrix
+*   `ishermitian::Bool = issymmetric & T<:Real` : whether `A` or `f` acts as a Hermitian
+    matrix
 *   `isposdef::Bool = false` : whether `A` or `f` acts as a positive definite matrix.
 For existing linear maps or matrices `A`, the default values will be taken by calling
 `issymmetric`, `ishermitian` and `isposdef` on the existing object `A`.
@@ -274,8 +278,8 @@ For the function-based constructor, there is one more keyword argument:
 *   `ismutating::Bool` : flags whether the function acts as a mutating matrix multiplication
     `f(y,x)` where the result vector `y` is the first argument (in case of `true`),
     or as a normal matrix multiplication that is called as `y=f(x)` (in case of `false`).
-    The default value is guessed by looking at the number of arguments of the first occurrence
-    of `f` in the method table.
+    The default value is guessed by looking at the number of arguments of the first
+    occurrence of `f` in the method table.
 """
 LinearMap(A::MapOrMatrix; kwargs...) = WrappedMap(A; kwargs...)
 LinearMap(J::UniformScaling, M::Int) = UniformScalingMap(J.λ, M)

src/left.jl

@@ -50,6 +50,7 @@ end
 function mul!(x::AbstractMatrix, y::AdjointAbsVecOrMat, A::LinearMap, α::Number, β::Number)
     check_dim_mul(x, y, A)
     _unsafe_mul!(x', conj(α)*A', y', true, conj(β))
+    # why not?: _unsafe_mul!(x', A', y', conj(α), conj(β))
     return x
 end
 
@@ -62,5 +63,6 @@ end
 function mul!(x::AbstractMatrix, y::TransposeAbsVecOrMat, A::LinearMap, α::Number, β::Number)
     check_dim_mul(x, y, A)
     _unsafe_mul!(transpose(x), α*transpose(A), transpose(y), true, β)
+    # why not?: _unsafe_mul!(transpose(x), transpose(A), transpose(y), α, β)
     return x
 end

src/scaledmap.jl

@@ -24,7 +24,7 @@ Base.adjoint(A::ScaledMap) = conj(A.λ) * adjoint(A.lmap)
 
 # comparison (sufficient, not necessary)
 Base.:(==)(A::ScaledMap, B::ScaledMap) =
-    (eltype(A) == eltype(B) && A.lmap == B.lmap) && A.λ == B.λ
+    eltype(A) == eltype(B) && A.lmap == B.lmap && A.λ == B.λ
 
 # scalar multiplication and division
 Base.:(*)(α::RealOrComplex, A::LinearMap) = ScaledMap(α, A)
@@ -43,12 +43,12 @@ Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B)
 Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ
 
 # multiplication with vectors/matrices
-for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
+for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
     @eval begin
-        function _unsafe_mul!(y::$outtype, A::ScaledMap, x::$intype)
+        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In)
             return _unsafe_mul!(y, A.lmap, x, A.λ, false)
         end
-        function _unsafe_mul!(y::$outtype, A::ScaledMap, x::$intype, α::Number, β::Number)
+        function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In, α::Number, β::Number)
             return _unsafe_mul!(y, A.lmap, x, A.λ * α, β)
         end
     end

src/transpose.jl

@@ -52,54 +52,52 @@ Base.:(==)(A::LinearMap, B::AdjointMap)      = ishermitian(A) && B.lmap == A
 
 # multiplication with vector/matrices
 # # TransposeMap
-function _unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector)
-    issymmetric(A.lmap) ? _unsafe_mul!(y, A.lmap, x) : error("transpose not implemented for $(A.lmap)")
-end
-function _unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix)
-    issymmetric(A.lmap) ? _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x)
-end
-function _unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector, α::Number, β::Number)
-    issymmetric(A.lmap) ? _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β)
-end
-function _unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix, α::Number, β::Number)
-    issymmetric(A.lmap) ? _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β)
-end
+_unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector) =
+    issymmetric(A.lmap) ?
+        _unsafe_mul!(y, A.lmap, x) : error("transpose not implemented for $(A.lmap)")
+_unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix) =
+    issymmetric(A.lmap) ?
+        _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x)
+_unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector, α::Number, β::Number)=
+    issymmetric(A.lmap) ?
+        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β)
+_unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix, α::Number, β::Number) =
+    issymmetric(A.lmap) ?
+        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β)
 # # AdjointMap
-function _unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector)
-    ishermitian(A.lmap) ? _unsafe_mul!(y, A.lmap, x) : error("adjoint not implemented for $(A.lmap)")
-end
-function _unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix)
-    ishermitian(A.lmap) ? _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x)
-end
-function _unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector, α::Number, β::Number)
-    ishermitian(A.lmap) ? _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β)
-end
-function _unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix, α::Number, β::Number)
-    ishermitian(A.lmap) ? _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β)
-end
+_unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector) =
+    ishermitian(A.lmap) ?
+        _unsafe_mul!(y, A.lmap, x) : error("adjoint not implemented for $(A.lmap)")
+_unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix) =
+    ishermitian(A.lmap) ?
+        _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x)
+_unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector, α::Number, β::Number) =
+    ishermitian(A.lmap) ?
+        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β)
+_unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix, α::Number, β::Number) =
+    ishermitian(A.lmap) ?
+        _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β)
+
 # # ConjugateMap
-for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
+const ConjugateMap = AdjointMap{<:Any, <:TransposeMap}
+# canonical order of adjoint followed by transpose
+LinearAlgebra.transpose(A::AdjointMap) = adjoint(transpose(A.lmap))
+for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
     @eval begin
-        function _unsafe_mul!(y::$outtype, Ac::AdjointMap{<:Any,<:TransposeMap}, x::$intype)
+        function _unsafe_mul!(y::$Out, Ac::ConjugateMap, x::$In)
             return _conjmul!(y, Ac.lmap.lmap, x)
         end
-        function _unsafe_mul!(y::$outtype, Ac::AdjointMap{<:Any,<:TransposeMap}, x::$intype, α::Number, β::Number)
+        function _unsafe_mul!(y::$Out, Ac::ConjugateMap, x::$In, α::Number, β::Number)
             return _conjmul!(y, Ac.lmap.lmap, x, α, β)
         end
-        function _unsafe_mul!(y::$outtype, At::TransposeMap{<:Any,<:AdjointMap}, x::$intype)
-            return _conjmul!(y, At.lmap.lmap, x)
-        end
-        function _unsafe_mul!(y::$outtype, At::TransposeMap{<:Any,<:AdjointMap}, x::$intype, α::Number, β::Number)
-            return _conjmul!(y, At.lmap.lmap, x, α, β)
-        end
     end
 end
 
 # multiplication helper function
 _conjmul!(y, A, x) = conj!(mul!(y, A, conj(x)))
 function _conjmul!(y, A, x::AbstractVector, α, β)
     xca = conj!(x * α)
-    z = A*xca
+    z = A * xca
     y .= y .* β + conj.(z)
     return y
 end

src/uniformscalingmap.jl

@@ -7,9 +7,11 @@ struct UniformScalingMap{T} <: LinearMap{T}
     end
 end
 UniformScalingMap(λ::Number, M::Int, N::Int) =
-    (M == N ? UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square"))
+    (M == N ?
+        UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square"))
 UniformScalingMap(λ::T, sz::Dims{2}) where {T} =
-    (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square"))
+    (sz[1] == sz[2] ?
+        UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square"))
 
 MulStyle(::UniformScalingMap) = FiveArg()
 
@@ -42,13 +44,13 @@ Base.:(*)(A::UniformScalingMap, x::AbstractVector) =
     length(x) == A.M ? A.λ * x : throw(DimensionMismatch("*"))
 
 # multiplication with vector/matrix
-for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
+for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
     @eval begin
-        function _unsafe_mul!(y::$outtype, J::UniformScalingMap, x::$intype)
+        function _unsafe_mul!(y::$Out, J::UniformScalingMap, x::$In)
             _scaling!(y, J.λ, x, true, false)
             return y
         end
-        function _unsafe_mul!(y::$outtype, J::UniformScalingMap, x::$intype,
+        function _unsafe_mul!(y::$Out, J::UniformScalingMap, x::$In,
                     α::Number, β::Number)
             _scaling!(y, J.λ, x, α, β)
             return y
@@ -58,41 +60,27 @@ end
 
 function _scaling!(y, λ, x, α, β)
     if (iszero(α) || iszero(λ))
-        iszero(β) && (fill!(y, zero(eltype(y))); return y)
+        iszero(β) && return fill!(y, zero(eltype(y)))
         isone(β) && return y
-        # β != 0, 1
-        rmul!(y, β)
-        return y
+        return rmul!(y, β)
+    elseif isone(α) && isone(λ)
+        iszero(β) && return copyto!(y, x)
+        isone(β) && return y .+= x
+        return y .= y .* β .+ x
     elseif isone(α)
-        if iszero(β)
-            isone(λ) && return copyto!(y, x)
-            y .= λ .* x
-            return y
-        elseif isone(β)
-            isone(λ) && return y .+= x
-            y .+= λ .* x
-            return y
-        else # β != 0, 1
-            isone(λ) && (axpby!(one(eltype(x)), x, β, y); return y)
-            y .= y .* β .+ λ .* x
-            return y
-        end
-    else # α != 0, 1
-        if iszero(β)
-            isone(λ) && (y .= x .* α; return y)
-            y .= λ .* x .* α
-            return y
-        elseif isone(β)
-            isone(λ) && (axpby!(α, x, β, y); return y)
-            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!
+        iszero(β) && return y .= λ .* x
+        isone(β) && return y .+= λ .* x
+        return y .= y .* β .+ λ .* x
+    elseif isone(λ)
+        iszero(β) && return y .= x .* α
+        isone(β) && return y .+= x .* α
+        return y .= y .* β .+ x .* α
+    else
+        iszero(β) && return y .= λ .* x .* α
+        isone(β) && return y .+= λ .* x .* α
+        return y .= y .* β .+ λ .* x .* α
+    end
+end
 
 # combine LinearMap and UniformScaling objects in linear combinations
 Base.:(+)(A₁::LinearMap, A₂::UniformScaling) = A₁ + UniformScalingMap(A₂.λ, size(A₁, 1))