From 12c13fbdd2c658a5e3444346312685711585f3db Mon Sep 17 00:00:00 2001 From: krcools Date: Thu, 2 Jun 2022 13:38:13 +0100 Subject: [PATCH] Specialisation of mul! for left multiplication with scalars (#173) Co-authored-by: Daniel Karrasch --- docs/Project.toml | 2 +- docs/src/custom.jl | 30 ++++++-- docs/src/history.md | 8 ++ docs/src/types.md | 2 + src/LinearMaps.jl | 118 +++++++++++++++++++++++------ src/blockmap.jl | 152 ++++++++++++++++++++++++++++++++------ src/conversion.jl | 58 ++++++--------- src/embeddedmap.jl | 11 +++ src/fillmap.jl | 7 ++ src/linearcombination.jl | 14 ++++ src/scaledmap.jl | 5 ++ src/transpose.jl | 12 +++ src/uniformscalingmap.jl | 10 +++ src/wrappedmap.jl | 5 ++ test/blockmap.jl | 22 ++++-- test/embeddedmap.jl | 3 +- test/fillmap.jl | 3 +- test/linearmaps.jl | 1 + test/transpose.jl | 6 +- test/uniformscalingmap.jl | 2 + 20 files changed, 374 insertions(+), 97 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 4c6dd40b..1a251c99 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,4 +8,4 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" [compat] BenchmarkTools = "1" Documenter = "0.25, 0.26, 0.27" -Literate = "2" \ No newline at end of file +Literate = "2" diff --git a/docs/src/custom.jl b/docs/src/custom.jl index 92a09588..1a52902f 100644 --- a/docs/src/custom.jl +++ b/docs/src/custom.jl @@ -29,10 +29,13 @@ end Base.size(A::MyFillMap) = A.size # By a couple of defaults provided for all subtypes of `LinearMap`, we only need to define -# a `LinearAlgebra.mul!` method to have minimal, operational type. +# a `LinearMaps._unsafe_mul!` method to have a minimal, operational type. The (internal) +# function `_unsafe_mul!` is called by `LinearAlgebra.mul!`, constructors, and conversions +# and only needs to be concerned with the bare computing kernel. Dimension checking is done +# on the level of `mul!` etc. Factoring out dimension checking is done to minimise overhead +# caused by repetitive checking. -function LinearAlgebra.mul!(y::AbstractVecOrMat, A::MyFillMap, x::AbstractVector) - LinearMaps.check_dim_mul(y, A, x) +function LinearMaps._unsafe_mul!(y::AbstractVecOrMat, A::MyFillMap, x::AbstractVector) return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x)) end @@ -84,7 +87,7 @@ using BenchmarkTools LinearMaps.MulStyle(A::MyFillMap) = FiveArg() -function LinearAlgebra.mul!( +function LinearMaps._unsafe_mul!( y::AbstractVecOrMat, A::MyFillMap, x::AbstractVector, @@ -126,7 +129,7 @@ typeof(A') try A'x catch e println(e) end # If the operator is symmetric or Hermitian, the transpose and the adjoint, respectively, -# of the linear map `A` is given by `A` itself. So let's define corresponding checks. +# of the linear map `A` is given by `A` itself. So let us define corresponding checks. LinearAlgebra.issymmetric(A::MyFillMap) = A.size[1] == A.size[2] LinearAlgebra.ishermitian(A::MyFillMap) = isreal(A.λ) && A.size[1] == A.size[2] @@ -154,22 +157,35 @@ try MyFillMap(5.0, (3, 4))' * ones(3) catch e println(e) end # The first option is to write `LinearAlgebra.mul!` methods for the corresponding wrapped # map types; for instance, -function LinearAlgebra.mul!( +function LinearMaps._unsafe_mul!( y::AbstractVecOrMat, transA::LinearMaps.TransposeMap{<:Any,<:MyFillMap}, x::AbstractVector ) - LinearMaps.check_dim_mul(y, transA, x) λ = transA.lmap.λ return fill!(y, iszero(λ) ? zero(eltype(y)) : transpose(λ)*sum(x)) end +# Now, the adjoint multiplication works. + +MyFillMap(5.0, (3, 4))' * ones(3) + # If you have set the `MulStyle` trait to `FiveArg()`, you should provide a corresponding # 5-arg `mul!` method for `LinearMaps.TransposeMap{<:Any,<:MyFillMap}` and # `LinearMaps.AdjointMap{<:Any,<:MyFillMap}`. # ### Path 2: Invariant `LinearMap` subtypes +# Before we start, let us delete the previously defined method to make sure we use the +# following definitions. + +Base.delete_method( + first(methods( + LinearMaps._unsafe_mul!, + (AbstractVecOrMat, LinearMaps.TransposeMap{<:Any,<:MyFillMap}, AbstractVector)) + ) +) + # The seconnd option is when your class of linear maps that are modelled by your custom # `LinearMap` subtype are invariant under taking adjoints and transposes. diff --git a/docs/src/history.md b/docs/src/history.md index a83d02e9..27909a3e 100644 --- a/docs/src/history.md +++ b/docs/src/history.md @@ -1,6 +1,14 @@ # Version history ## What's new in v3.7 +* `mul!(M::AbstractMatrix, A::LinearMap, s::Number, a, b)` methods are provided, mimicking + similar methods in `Base.LinearAlgebra`. This version allows for the memory efficient + implementation of in-place addition and conversion of a `LinearMap` to `Matrix`. + Efficient specialisations for `WrappedMap`, `ScaledMap`, and `LinearCombination` are + provided. If users supply the corresponding `_unsafe_mul!` method for their custom maps, + conversion, construction, and inplace addition will benefit from this supplied efficient + implementation. If no specialisation is supplied, a generic fallback is used that is based + on feeding the canonical basis of unit vectors to the linear map. * A new map type called `EmbeddedMap` is introduced. It is a wrapper of a "small" `LinearMap` (or a suitably converted `AbstractVecOrMat`) embedded into a "larger" zero map. Hence, diff --git a/docs/src/types.md b/docs/src/types.md index 7191fe34..e8330316 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -103,6 +103,8 @@ Base.:*(::AbstractMatrix,::LinearMap) LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::AbstractVector) LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::AbstractVector,::Number,::Number) LinearAlgebra.mul!(::AbstractMatrix,::AbstractMatrix,::LinearMap) +LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::Number) +LinearAlgebra.mul!(::AbstractVecOrMat,::LinearMap,::Number,::Number,::Number) *(::LinearAlgebra.AdjointAbsVec,::LinearMap) *(::LinearAlgebra.TransposeAbsVec,::LinearMap) ``` diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index c17ef46f..52504975 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -143,14 +143,14 @@ 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]; Y = similar(B); mul!(Y, A, B); +julia> A=LinearMap([1.0 2.0; 3.0 4.0]); B=ones(2); Y = similar(B); mul!(Y, A, B); julia> Y 2-element Array{Float64,1}: 3.0 7.0 -julia> A=LinearMap([1.0 2.0; 3.0 4.0]); B=[1.0 1.0; 1.0 1.0]; Y = similar(B); mul!(Y, A, B); +julia> A=LinearMap([1.0 2.0; 3.0 4.0]); B=ones(4,4); Y = similar(B); mul!(Y, A, B); julia> Y 2×2 Array{Float64,2}: @@ -162,6 +162,35 @@ function mul!(y::AbstractVecOrMat, A::LinearMap, x::AbstractVector) check_dim_mul(y, A, x) return _unsafe_mul!(y, A, x) end +# the following is of interest in, e.g., subspace-iteration methods +function mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix) + check_dim_mul(Y, A, X) + return _unsafe_mul!(Y, A, X) +end + +""" + mul!(Y::AbstractMatrix, A::LinearMap, b::Number) -> Y + +Scales the matrix representation of the linear map `A` by `b` and stores the result in `Y`, +overwriting the existing value of `Y`. + +## Examples +```jldoctest; setup=(using LinearAlgebra, LinearMaps) +julia> A = LinearMap{Int}(cumsum, 3); b = 2; Y = Matrix{Int}(undef, (3,3)); + +julia> mul!(Y, A, b) +3×3 Matrix{Int64}: + 2 0 0 + 2 2 0 + 2 2 2 +``` +""" +function mul!(y::AbstractVecOrMat, A::LinearMap, s::Number) + size(y) == size(A) || + throw( + DimensionMismatch("y has size $(size(y)), A has size $(size(A)).")) + return _unsafe_mul!(y, A, s) +end """ mul!(C::AbstractVecOrMat, A::LinearMap, B::AbstractVector, α, β) -> C @@ -197,6 +226,46 @@ function mul!(y::AbstractVecOrMat, A::LinearMap, x::AbstractVector, α::Number, check_dim_mul(y, A, x) return _unsafe_mul!(y, A, x, α, β) end +function mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix, α::Number, β::Number) + check_dim_mul(Y, A, X) + return _unsafe_mul!(Y, A, X, α, β) +end + +""" + mul!(Y::AbstractMatrix, A::LinearMap, b::Number, α::Number, β::Number) -> Y + +Scales the matrix representation of the linear map `A` by `b*α`, adds the result to `Y*β` +and stores the final result in `Y`, overwriting the existing value of `Y`. + +## Examples +```jldoctest; setup=(using LinearAlgebra, LinearMaps) +julia> A = LinearMap{Int}(cumsum, 3); b = 2; Y = ones(Int, (3,3)); + +julia> mul!(Y, A, b, 2, 1) +3×3 Matrix{Int64}: + 5 1 1 + 5 5 1 + 5 5 5 +``` +""" +function mul!(y::AbstractMatrix, A::LinearMap, s::Number, α::Number, β::Number) + size(y) == size(A) || + throw( + DimensionMismatch("y has size $(size(y)), A has size $(size(A)).")) + return _unsafe_mul!(y, A, s, α, β) +end + +_unsafe_mul!(y, A::MapOrVecOrMat, x) = mul!(y, A, x) +_unsafe_mul!(y, A::AbstractVecOrMat, x, α, β) = mul!(y, A, x, α, β) +_unsafe_mul!(y::AbstractVecOrMat, A::LinearMap, x::AbstractVector, α, β) = + _generic_mapvec_mul!(y, A, x, α, β) +_unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix) = + _generic_mapmat_mul!(y, A, x) +_unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix, α::Number, β::Number) = + _generic_mapmat_mul!(y, A, x, α, β) +_unsafe_mul!(Y::AbstractMatrix, A::LinearMap, s::Number) = _generic_mapnum_mul!(Y, A, s) +_unsafe_mul!(Y::AbstractMatrix, A::LinearMap, s::Number, α::Number, β::Number) = + _generic_mapnum_mul!(Y, A, s, α, β) function _generic_mapvec_mul!(y, A, x, α, β) # this function needs to call mul! for, e.g., AdjointMap{...,<:CustomMap} @@ -226,33 +295,40 @@ function _generic_mapvec_mul!(y, A, x, α, β) end end -# the following is of interest in, e.g., subspace-iteration methods -function mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix) - check_dim_mul(Y, A, X) - return _unsafe_mul!(Y, A, X) -end -function mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix, α::Number, β::Number) - check_dim_mul(Y, A, X) - return _unsafe_mul!(Y, A, X, α, β) +function _generic_mapmat_mul!(Y, A, X) + for (Xi, Yi) in zip(eachcol(X), eachcol(Y)) + mul!(Yi, A, Xi) + end + return Y end - -function _generic_mapmat_mul!(Y, A, X, α=true, β=false) +function _generic_mapmat_mul!(Y, A, X, α, β) for (Xi, Yi) in zip(eachcol(X), eachcol(Y)) mul!(Yi, A, Xi, α, β) end return Y end -_unsafe_mul!(y, A::MapOrVecOrMat, x) = mul!(y, A, x) -_unsafe_mul!(y, A::AbstractVecOrMat, x, α, β) = mul!(y, A, x, α, β) -function _unsafe_mul!(y::AbstractVecOrMat, A::LinearMap, x::AbstractVector, α, β) - return _generic_mapvec_mul!(y, A, x, α, β) -end -function _unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix) - return _generic_mapmat_mul!(y, A, x) +function _generic_mapnum_mul!(Y, A, s) + T = promote_type(eltype(A), typeof(s)) + ax2 = axes(A)[2] + xi = zeros(T, ax2) + @inbounds for (i, Yi) in zip(ax2, eachcol(Y)) + xi[i] = s + mul!(Yi, A, xi) + xi[i] = zero(T) + end + return Y end -function _unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix, α, β) - return _generic_mapmat_mul!(y, A, x, α, β) +function _generic_mapnum_mul!(Y, A, s, α, β) + T = promote_type(eltype(A), typeof(s)) + ax2 = axes(A)[2] + xi = zeros(T, ax2) + @inbounds for (i, Yi) in zip(ax2, eachcol(Y)) + xi[i] = s + mul!(Yi, A, xi, α, β) + xi[i] = zero(T) + end + return Y end include("left.jl") # left multiplication by a transpose or adjoint vector diff --git a/src/blockmap.jl b/src/blockmap.jl index afc7d8fd..67de2145 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -283,7 +283,7 @@ Base.:(==)(A::BlockMap, B::BlockMap) = # multiplication helper functions ############ -function _blockmul!(y, A::BlockMap, x, α, β) +function _blockmul!(y, A, x, α, β) if iszero(α) iszero(β) && return fill!(y, zero(eltype(y))) isone(β) && return y @@ -292,11 +292,37 @@ function _blockmul!(y, A::BlockMap, x, α, β) return __blockmul!(MulStyle(A), y, A, x, α, β) end -# provide one global intermediate storage vector if necessary -__blockmul!(::FiveArg, y, A, x, α, β) = ___blockmul!(y, A, x, α, β, nothing) -__blockmul!(::ThreeArg, y, A, x, α, β) = ___blockmul!(y, A, x, α, β, similar(y)) +__blockmul!(::MulStyle, y, A, x::Number, α, β) = ___blockmul!(y, A, x, α, β) +function ___blockmul!(y, A, x::Number, α, β) + maps, rows, yinds, xinds = A.maps, A.rows, A.rowranges, A.colranges + mapind = 0 + if iszero(β) + s = x*α + for (row, yi) in zip(rows, yinds) + mapind += 1 + _unsafe_mul!(view(y, yi, xinds[mapind]), maps[mapind], s) + for _ in 2:row + mapind += 1 + _unsafe_mul!(view(y, yi, xinds[mapind]), maps[mapind], s) + end + end + else + for (row, yi) in zip(rows, yinds) + mapind += 1 + _unsafe_mul!(view(y, yi, xinds[mapind]), maps[mapind], x, α, β) + for _ in 2:row + mapind += 1 + _unsafe_mul!(view(y, yi, xinds[mapind]), maps[mapind], x, α, β) + end + end + end + return y +end -function ___blockmul!(y, A, x, α, β, ::Nothing) +# provide one global intermediate storage vector if necessary +__blockmul!(::FiveArg, y, A, x::AbstractVecOrMat, α, β) = ___blockmul!(y, A, x, α, β, nothing) +__blockmul!(::ThreeArg, y, A, x::AbstractVecOrMat, α, β) = ___blockmul!(y, A, x, α, β, similar(y)) +function ___blockmul!(y, A, x::AbstractVecOrMat, α, β, ::Nothing) maps, rows, yinds, xinds = A.maps, A.rows, A.rowranges, A.colranges mapind = 0 for (row, yi) in zip(rows, yinds) @@ -310,7 +336,7 @@ function ___blockmul!(y, A, x, α, β, ::Nothing) end return y end -function ___blockmul!(y, A, x, α, β, z) +function ___blockmul!(y, A, x::AbstractVecOrMat, α, β, z) maps, rows, yinds, xinds = A.maps, A.rows, A.rowranges, A.colranges mapind = 0 for (row, yi) in zip(rows, yinds) @@ -333,30 +359,65 @@ function ___blockmul!(y, A, x, α, β, z) return y end -function _transblockmul!(y, A::BlockMap, x, α, β, transform) - maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges +function _transblockmul!(y, A, x, α, β, transform) if iszero(α) iszero(β) && return fill!(y, zero(eltype(y))) isone(β) && return y return rmul!(y, β) else + return __transblockmul!(y, A, x, α, β, transform) + end +end +function __transblockmul!(y, A, x::Number, α, β, transform) + maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges + if iszero(β) + s = x*α # first block row (rowind = 1) of A, meaning first block column of A', fill all of y - xrow = selectdim(x, 1, first(xinds)) for rowind in 1:first(rows) - yrow = selectdim(y, 1, yinds[rowind]) - _unsafe_mul!(yrow, transform(maps[rowind]), xrow, α, β) + _unsafe_mul!(view(y, yinds[rowind], first(xinds)), transform(maps[rowind]), s) + end + mapind = first(rows) + # subsequent block rows of A (block columns of A') + @inbounds for i in 2:length(rows), _ in 1:rows[i] + mapind +=1 + _unsafe_mul!(view(y, yinds[mapind], xinds[i]), transform(maps[mapind]), s) + end + else + # first block row (rowind = 1) of A, meaning first block column of A', fill all of y + for rowind in 1:first(rows) + ytile = view(y, yinds[rowind], first(xinds)) + _unsafe_mul!(ytile, transform(maps[rowind]), x, α, β) end mapind = first(rows) # subsequent block rows of A (block columns of A'), # add results to corresponding parts of y # TODO: think about multithreading - @inbounds for i in 2:length(rows) - xrow = selectdim(x, 1, xinds[i]) - for _ in 1:rows[i] - mapind +=1 - yrow = selectdim(y, 1, yinds[mapind]) - _unsafe_mul!(yrow, transform(maps[mapind]), xrow, α, true) - end + @inbounds for i in 2:length(rows), _ in 1:rows[i] + mapind +=1 + ytile = view(y, yinds[mapind], xinds[i]) + _unsafe_mul!(ytile, transform(maps[mapind]), x, α, β) + end + end + return y +end +function __transblockmul!(y, A, x, α, β, transform) + maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges + # first block row (rowind = 1) of A, meaning first block column of A', fill all of y + xrow = selectdim(x, 1, first(xinds)) + for rowind in 1:first(rows) + yrow = selectdim(y, 1, yinds[rowind]) + _unsafe_mul!(yrow, transform(maps[rowind]), xrow, α, β) + end + mapind = first(rows) + # subsequent block rows of A (block columns of A'), + # add results to corresponding parts of y + # TODO: think about multithreading + @inbounds for i in 2:length(rows) + xrow = selectdim(x, 1, xinds[i]) + for _ in 1:rows[i] + mapind +=1 + yrow = selectdim(y, 1, yinds[mapind]) + _unsafe_mul!(yrow, transform(maps[mapind]), xrow, α, true) end end return y @@ -393,6 +454,24 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM end end +############ +# multiplication with a scalar +############ + +function _unsafe_mul!(Y::AbstractMatrix, A::BlockMap, s::Number, α::Number=true, β::Number=false) + require_one_based_indexing(Y, s) + return _blockmul!(Y, A, s, α, β) +end +for (MT, transform) in ((:TransposeMap, :transpose), (:AdjointMap, :adjoint)) + @eval begin + function _unsafe_mul!(Y::AbstractMatrix, wrapA::$MT{<:Any, <:BlockMap}, s::Number, + α::Number=true, β::Number=false) + require_one_based_indexing(Y) + return _transblockmul!(Y, wrapA.lmap, s, α, β, $transform) + end + end +end + ############ # BlockDiagonalMap ############ @@ -478,11 +557,11 @@ LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} = Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) = (eltype(A) == eltype(B) && all(A.maps .== B.maps)) -for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) +for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix), (Number, AbstractMatrix)) @eval begin function _unsafe_mul!(y::$Out, A::BlockDiagonalMap, x::$In) require_one_based_indexing(y, x) - return _blockscaling!(y, A, x, true, false) + return _blockscaling!(y, A, x) end function _unsafe_mul!(y::$Out, A::BlockDiagonalMap, x::$In, α::Number, β::Number) require_one_based_indexing(y, x) @@ -491,11 +570,38 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM end end -function _blockscaling!(y, A::BlockDiagonalMap, x, α, β) +function _blockscaling!(y, A, x::Number) + maps, yinds, xinds = A.maps, A.rowranges, A.colranges + fill!(y, zero(eltype(y))) + # TODO: think about multi-threading here + @inbounds for (yind, map, xind) in zip(yinds, maps, xinds) + _unsafe_mul!(view(y, yind, xind), map, x) + end + return y +end +function _blockscaling!(y, A, x::Number, α, β) + maps, yinds, xinds = A.maps, A.rowranges, A.colranges + LinearAlgebra._rmul_or_fill!(y, β) + # TODO: think about multi-threading here + @inbounds for (yind, map, xind) in zip(yinds, maps, xinds) + _unsafe_mul!(view(y, yind, xind), map, x, α, true) + end + return y +end + +function _blockscaling!(y, A, x) + maps, yinds, xinds = A.maps, A.rowranges, A.colranges + # TODO: think about multi-threading here + @inbounds for (yind, map, xind) in zip(yinds, maps, xinds) + _unsafe_mul!(selectdim(y, 1, yind), map, selectdim(x, 1, xind)) + end + return y +end +function _blockscaling!(y, A, x, α, β) maps, yinds, xinds = A.maps, A.rowranges, A.colranges # TODO: think about multi-threading here - @inbounds for i in 1:length(maps) - _unsafe_mul!(selectdim(y, 1, yinds[i]), maps[i], selectdim(x, 1, xinds[i]), α, β) + @inbounds for (yind, map, xind) in zip(yinds, maps, xinds) + _unsafe_mul!(selectdim(y, 1, yind), map, selectdim(x, 1, xind), α, β) end return y end diff --git a/src/conversion.jl b/src/conversion.jl index 8262fdf7..810dc976 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -1,21 +1,22 @@ # Matrix: create matrix representation of LinearMap function Base.Matrix{T}(A::LinearMap) where {T} - M, N = size(A) - mat = Matrix{T}(undef, (M, N)) - v = fill(zero(T), N) - @inbounds for i in 1:N - v[i] = one(T) - # need mul!, e.g., for TransposeMap{<:CustomMap} - mul!(view(mat, :, i), A, v) - v[i] = zero(T) - end - return mat + mat = Matrix{T}(undef, size(A)) + mul!(mat, A, one(T)) end + +function Base.AbstractMatrix{T}(A::LinearMap) where {T} + mat = similar(Array{T}, axes(A)) + mul!(mat, A, one(T)) +end + Base.Matrix(A::LinearMap{T}) where {T} = Matrix{T}(A) +Base.AbstractMatrix(A::LinearMap{T}) where {T} = AbstractMatrix{T}(A) + Base.Array(A::LinearMap) = Matrix(A) Base.convert(::Type{T}, A::LinearMap) where {T<:Matrix} = T(A) Base.convert(::Type{Array}, A::LinearMap) = convert(Matrix, A) -Base.convert(::Type{AbstractMatrix}, A::LinearMap) = convert(Matrix, A) +# Base.convert(::Type{AbstractMatrix}, A::LinearMap) = convert(Matrix, A) +Base.convert(::Type{AbstractMatrix}, A::LinearMap) = AbstractMatrix(A) Base.convert(::Type{AbstractArray}, A::LinearMap) = convert(AbstractMatrix, A) # sparse: create sparse matrix representation of LinearMap @@ -48,13 +49,10 @@ SparseArrays.SparseMatrixCSC(A::LinearMap) = sparse(A) # special cases # ScaledMap -Base.Matrix{T}(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) where {T} = - convert(Matrix{T}, A.λ * A.lmap.lmap) SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) = A.λ * sparse(A.lmap.lmap) # UniformScalingMap -Base.Matrix{T}(J::UniformScalingMap) where {T} = Matrix{T}(J.λ*I, size(J)) Base.convert(::Type{AbstractMatrix}, J::UniformScalingMap) = Diagonal(fill(J.λ, J.M)) # WrappedMap @@ -73,24 +71,18 @@ for (T, t) in ((AdjointMap, adjoint), (TransposeMap, transpose)) end # LinearCombination -function Base.Matrix{T}(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}}) where {T} - maps = ΣA.maps - mats = map(A->getfield(A, :lmap), maps) - return Matrix{T}(sum(mats)) -end function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}}) - maps = ΣA.maps - mats = map(A->getfield(A, :lmap), maps) - return convert(SparseMatrixCSC, sum(mats)) + mats = map(A->getfield(A, :lmap), ΣA.maps) + return sum(sparse, mats) end # CompositeMap function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, LinearMap}}) where {T} B, A = AB.maps require_one_based_indexing(B) - Y = Matrix{eltype(AB)}(undef, size(AB)) - @views for i in 1:size(Y, 2) - _unsafe_mul!(Y[:, i], A, B.lmap[:, i]) + Y = Matrix{T}(undef, size(AB)) + for (yi, bi) in zip(eachcol(Y), eachcol(B.lmap)) + _unsafe_mul!(yi, A, bi) end return Y end @@ -105,33 +97,30 @@ for ((TA, fieldA), (TB, fieldB)) in (((VecOrMatMap, :lmap), (VecOrMatMap, :lmap) end function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap}}) where {T} B, A = AB.maps - return convert(Matrix{T}, A.lmap*B.lmap) + return mul!(Matrix{T}(undef, size(AB)), A.lmap, B.lmap) end function SparseArrays.sparse(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap}}) B, A = AB.maps - return convert(SparseMatrixCSC, A.lmap*B.lmap) + return sparse(A.lmap)*sparse(B.lmap) end function Base.Matrix{T}(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}}) where {T} A, J = λA.maps - return convert(Matrix{T}, J.λ*A.lmap) + return mul!(Matrix{T}(undef, size(λA)), J.λ, A.lmap) end function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}}) A, J = λA.maps - return convert(SparseMatrixCSC, J.λ*A.lmap) + return J.λ*sparse(A.lmap) end function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}}) where {T} J, A = Aλ.maps - return convert(Matrix{T}, A.lmap*J.λ) + return mul!(Matrix{T}(undef, size(Aλ)), A.lmap, J.λ) end function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}}) J, A = Aλ.maps - return convert(SparseMatrixCSC, A.lmap*J.λ) + return sparse(A.lmap)*J.λ end # BlockMap & BlockDiagonalMap -Base.Matrix{T}(A::BlockMap) where {T} = hvcat(A.rows, convert.(Matrix{T}, A.maps)...) -Base.convert(::Type{AbstractMatrix}, A::BlockMap) = - hvcat(A.rows, convert.(AbstractMatrix, A.maps)...) function SparseArrays.sparse(A::BlockMap) return hvcat( A.rows, @@ -139,7 +128,6 @@ function SparseArrays.sparse(A::BlockMap) convert.(AbstractArray, _tail(A.maps))... ) end -Base.Matrix{T}(A::BlockDiagonalMap) where {T} = Base._cat((1,2), convert.(Matrix{T}, A.maps)...) Base.convert(::Type{AbstractMatrix}, A::BlockDiagonalMap) = sparse(A) function SparseArrays.sparse(A::BlockDiagonalMap) return blockdiag(convert.(SparseMatrixCSC, A.maps)...) diff --git a/src/embeddedmap.jl b/src/embeddedmap.jl index fb48e437..1e659a38 100644 --- a/src/embeddedmap.jl +++ b/src/embeddedmap.jl @@ -81,3 +81,14 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM return y end end + +function _unsafe_mul!(Y::AbstractMatrix, A::EmbeddedMap, x::Number) + fill!(Y, zero(eltype(Y))) + _unsafe_mul!(view(Y, A.rows, A.cols), A.lmap, x) + return Y +end +function _unsafe_mul!(Y::AbstractMatrix, A::EmbeddedMap, x::Number, α::Number, β::Number) + LinearAlgebra._rmul_or_fill!(Y, β) + _unsafe_mul!(view(Y, A.rows, A.cols), A.lmap, x, α, !iszero(β)) + return Y +end diff --git a/src/fillmap.jl b/src/fillmap.jl index c415740f..ea017c83 100644 --- a/src/fillmap.jl +++ b/src/fillmap.jl @@ -36,6 +36,13 @@ function _unsafe_mul!(y::AbstractVecOrMat, A::FillMap, x::AbstractVector) return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x)) end +_unsafe_mul!(Y::AbstractMatrix, A::FillMap, x::Number) = fill!(Y, A.λ*x) +function _unsafe_mul!(Y::AbstractMatrix, A::FillMap, x::Number, α::Number, β::Number) + LinearAlgebra._rmul_or_fill!(Y, β) + Y .+= A.λ*x*α + return Y +end + function _unsafe_mul!(y::AbstractVecOrMat, A::FillMap, x::AbstractVector, α::Number, β::Number) if iszero(α) !isone(β) && rmul!(y, β) diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 9b0078e0..b0f8d21f 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -112,6 +112,20 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM end end +function _unsafe_mul!(M::AbstractMatrix, L::LinearCombination, s::Number) + _unsafe_mul!(M, first(L.maps), s) + _mul!(MulStyle(L), M, L, s, true) + return M +end + +function _unsafe_mul!(M::AbstractMatrix, L::LinearCombination, s::Number, α::Number, β::Number) + LinearAlgebra._rmul_or_fill!(M, β) + for map in L.maps + _unsafe_mul!(M, map, s, α, true) + end + return M +end + _mul!(::FiveArg, y, A::LinearCombination, x, α) = __mul!(y, _tail(A.maps), x, α, nothing) _mul!(::ThreeArg, y, A::LinearCombination, x, α) = __mul!(y, _tail(A.maps), x, α, similar(y)) diff --git a/src/scaledmap.jl b/src/scaledmap.jl index e62aa014..a2efdef9 100644 --- a/src/scaledmap.jl +++ b/src/scaledmap.jl @@ -66,3 +66,8 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), end end end + +_unsafe_mul!(Y::AbstractMatrix, X::ScaledMap, c::Number) = + _unsafe_mul!(Y, X.lmap, X.λ*c) +_unsafe_mul!(Y::AbstractMatrix, X::ScaledMap, c::Number, α::Number, β::Number) = + _unsafe_mul!(Y, X.lmap, X.λ*c, α, β) diff --git a/src/transpose.jl b/src/transpose.jl index 97f61d04..89a8b45c 100644 --- a/src/transpose.jl +++ b/src/transpose.jl @@ -57,12 +57,18 @@ _unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector) = _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::AbstractMatrix, A::TransposeMap, x::Number) = + issymmetric(A.lmap) ? + _unsafe_mul!(y, A.lmap, x) : _generic_mapnum_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, α, β) +_unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::Number, α::Number, β::Number) = + issymmetric(A.lmap) ? + _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapnum_mul!(y, A, x, α, β) # # AdjointMap _unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector) = ishermitian(A.lmap) ? @@ -70,12 +76,18 @@ _unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector) = _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::AbstractMatrix, A::AdjointMap, x::Number) = + ishermitian(A.lmap) ? + _unsafe_mul!(y, A.lmap, x) : _generic_mapnum_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, α, β) +_unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::Number, α::Number, β::Number) = + ishermitian(A.lmap) ? + _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapnum_mul!(y, A, x, α, β) # # ConjugateMap const ConjugateMap = AdjointMap{<:Any, <:TransposeMap} diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 5530f896..fae06a35 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -70,6 +70,16 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM end end +function _unsafe_mul!(M::AbstractMatrix, L::UniformScalingMap, s::Number, + α::Number=true, β::Number=false) + LinearAlgebra._rmul_or_fill!(M, β) + c = L.λ * s * α + for (i,j) in zip(axes(L)...) + M[i,j] += c + end + return M +end + function _scaling!(y, λ, x, α, β) if (iszero(α) || iszero(λ)) iszero(β) && return fill!(y, zero(eltype(y))) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 356d1383..72c5c476 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -110,6 +110,11 @@ end mul!(X::AbstractMatrix, Y::AbstractMatrix, A::VecOrMatMap, α::Number, β::Number) = mul!(X, Y, A.lmap, α, β) + +_unsafe_mul!(Y::AbstractMatrix, A::VecOrMatMap, s::Number) = _unsafe_mul!(Y, A.lmap, s) +_unsafe_mul!(Y::AbstractMatrix, A::VecOrMatMap, s::Number, α::Number, β::Number) = + _unsafe_mul!(Y, A.lmap, s, α, β) + # the following 2 methods are needed for disambiguation with left-multiplication function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::AbstractMatrix{<:RealOrComplex}, A::VecOrMatMap{<:RealOrComplex}, α::RealOrComplex, β::RealOrComplex) diff --git a/test/blockmap.jl b/test/blockmap.jl index e8f8b5e8..36c016f6 100644 --- a/test/blockmap.jl +++ b/test/blockmap.jl @@ -23,11 +23,11 @@ using LinearMaps: FiveArg A = [A11 A12] x = rand(n+n2) @test size(L) == size(A) == size(Lv) - @test Matrix(L) == A == Matrix(Lv) + @test Matrix(L) == A == Matrix(Lv) == mul!(copy(A), L, 1, true, false) @test L * x ≈ A * x ≈ Lv * x L = @inferred hcat(LinearMap(A11), LinearMap(A12), LinearMap(A11)) A = [A11 A12 A11] - @test Matrix(L) == A + @test Matrix(L) == A == mul!(zero(A), L, 1, true, false) A = [I I I A11 A11 A11 a] @test (@which [A11 A11 A11]).module != LinearMaps @test (@which [I I I A11 A11 A11]).module != LinearMaps @@ -72,7 +72,8 @@ using LinearMaps: FiveArg A = [A11; A21] x = rand(elty, n) @test size(L) == size(A) - @test Matrix(L) == Matrix(Lv) == A + @test Matrix(L) == A == mul!(copy(A), L, 1, true, false) + @test Matrix(Lv) == A == mul!(copy(A), Lv, 1, true, false) @test L * x ≈ Lv * x ≈ A * x A = [I; I; I; A11; A11; A11; reduce(hcat, fill(v, n))] @test (@which [I; I; I; A11; A11; A11; v v v v v v v v v v]).module != LinearMaps @@ -109,7 +110,10 @@ using LinearMaps: FiveArg @test L isa LinearMaps.BlockMap{elty} @test size(L) == size(A) @test L * x ≈ Lv * x ≈ A * x - @test Matrix(L) == Matrix(Lv) == A + @test Matrix(L) == Matrix(Lv) == A == mul!(zero(A), L, 1) + for α in (false, 1, rand()), β in (false, 1, rand()) + @test mul!(copy(A), L, 2, α, β) ≈ A*(2α + β) + end @test convert(AbstractMatrix, L) == A A = [I A12; A21 I] @test (@which [I A12; A21 I]).module != LinearMaps @@ -184,6 +188,9 @@ using LinearMaps: FiveArg @test Lt * x ≈ transform(A) * x @test convert(AbstractMatrix, Lt) == transform(A) @test sparse(transform(L)) == transform(A) + for α in (false, 1, rand()), β in (false, 1, rand()) + @test mul!(copy(transform(A)), Lt, 2, α, β) ≈ transform(A)*(2α + β) + end Lt = @inferred transform(LinearMap(L)) @test Lt * x ≈ transform(A) * x @test Matrix(Lt) == Matrix(transform(A)) @@ -210,10 +217,10 @@ using LinearMaps: FiveArg M2 = randn(elty, m, n+2); L2 = LinearMap(M2) M3 = randn(elty, m, n+3); L3 = LinearMap(M3) - # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse: if elty <: Complex @test_throws ErrorException LinearMaps.BlockDiagonalMap{Float64}((L1, L2, L3, L2, L1)) end + # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse: Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...)) @test (@which blockdiag(sparse.((M1, M2, M3, M2, M1))...)).module != LinearMaps @test (@which cat(M1, M2, M3, M2, M1; dims=(1,2))).module != LinearMaps @@ -224,7 +231,10 @@ using LinearMaps: FiveArg @test Bdv.maps isa Vector @test @inferred(LinearMaps.MulStyle(Bd)) === FiveArg() @test occursin("$(5m)×$(5n+9) LinearMaps.BlockDiagonalMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), Bd)) - @test Matrix(Bd) == Md + @test Matrix(Bd) == Md == mul!(zero(Md), Bd, 1) + for α in (false, rand()), β in (false, rand()) + @test mul!(copy(Md), Bd, 2, α, β) ≈ Md*(2α + β) + end @test convert(AbstractMatrix, Bd) isa SparseMatrixCSC @test sparse(Bd) == Md @test Matrix(@inferred blockdiag(L1)) == M1 diff --git a/test/embeddedmap.jl b/test/embeddedmap.jl index 754451ff..8398fb72 100644 --- a/test/embeddedmap.jl +++ b/test/embeddedmap.jl @@ -9,7 +9,8 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays zeros(size(M,1), offset[2]) M] BL = @inferred LinearMap(L, size(BM); offset=offset) s1, s2 = size(BM) - @test (@inferred Matrix(BL)) == BM + @test (@inferred Matrix(BL)) == BM == mul!(zero(BM), BL, 1) + @test mul!(copy(BM), BL, 2, true, true) == 3BM @test (@inferred Matrix(BL')) == BM' @test (@inferred Matrix(transpose(BL))) == transpose(BM) diff --git a/test/fillmap.jl b/test/fillmap.jl index 03ce8d8e..7cc0de58 100644 --- a/test/fillmap.jl +++ b/test/fillmap.jl @@ -16,7 +16,8 @@ using LinearMaps, LinearAlgebra, Test @test size(L) == (M, N) @test adjoint(L) == FillMap(adjoint(λ), (3,2)) @test transpose(L) == FillMap(λ, (3,2)) - @test Matrix(L) == A + @test Matrix(L) == A == mul!(copy(A), L, 1) == mul!(copy(A), L, 1, true, false) + @test mul!(copy(1A), L, 2, true, true) == 3A @test L * x ≈ A * x @test mul!(w, L, x) ≈ A * x @test mul!(W, L, X) ≈ A * X diff --git a/test/linearmaps.jl b/test/linearmaps.jl index 713c88d5..2783cfe6 100644 --- a/test/linearmaps.jl +++ b/test/linearmaps.jl @@ -77,6 +77,7 @@ LinearAlgebra.mul!(y::AbstractVector, A::Union{SimpleFunctionMap,SimpleComplexFu @test mul!(copy(W), F, V) ≈ L*V @test mul!(copy(W), F, V, α, β) ≈ L*V*α + W*β @test mul!(copy(W), F, V, 0, β) ≈ W*β + @test mul!(copy(FM), F, 1, true, true) == 2FM Fs = sparse(F) @test SparseMatrixCSC(F) == Fs == L diff --git a/test/transpose.jl b/test/transpose.jl index ef240505..17c81604 100644 --- a/test/transpose.jl +++ b/test/transpose.jl @@ -13,7 +13,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays @test !(transpose(CS) == adjoint(CS)) @test !(adjoint(CS) == transpose(CS)) M = Matrix(CS) - @test M == LowerTriangular(ones(10,10)) + @test M == LowerTriangular(ones(10,10)) == mul!(copy(M), CS, 1, true, false) x = rand(ComplexF64, 10); w = rand(ComplexF64, 10) X = rand(ComplexF64, 10, 3); W = rand(ComplexF64, 10, 3) α, β = rand(ComplexF64, 2) @@ -26,14 +26,16 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays @test mul!(copy(w), transform(CS), x, α, β) ≈ transform(M)*x*α + w*β @test mul!(copy(W), transform(CS), X, α, β) ≈ transform(M)*X*α + W*β end - for transform1 in (adjoint, transpose), transform2 in (adjoint, transpose) + for transform1 in (adjoint, transpose), transform2 in (identity, adjoint, transpose) @test transform1(transform1(CS)) === CS @test transform2(transform1(transform2(CS))) === transform1(CS) @test transform2(transform1(CS)) * x ≈ transform2(transform1(M))*x @test mul!(copy(w), transform2(transform1(CS)), x) ≈ transform2(transform1(M))*x @test mul!(copy(W), transform2(transform1(CS)), X) ≈ transform2(transform1(M))*X + @test mul!(copy(M), transform2(transform1(CS)), 2) ≈ transform2(transform1(M))*2 @test mul!(copy(w), transform2(transform1(CS)), x, α, β) ≈ transform2(transform1(M))*x*α + w*β @test mul!(copy(W), transform2(transform1(CS)), X, α, β) ≈ transform2(transform1(M))*X*α + W*β + @test mul!(copy(M), transform2(transform1(CS)), 2, α, β) ≈ transform2(transform1(M))*(2*α) + M*β end id = @inferred LinearMap(identity, identity, 10; issymmetric=true, ishermitian=true, isposdef=true) diff --git a/test/uniformscalingmap.jl b/test/uniformscalingmap.jl index 3a0f177a..840c04df 100644 --- a/test/uniformscalingmap.jl +++ b/test/uniformscalingmap.jl @@ -50,4 +50,6 @@ using Test, LinearMaps, LinearAlgebra @test X*Id isa LinearMap @test Matrix(Id*X) == X @test Matrix(X*Id) == X + @test Matrix(Id + X) == I + X + @test Matrix(X + Id) == X + I end