diff --git a/docs/src/types.md b/docs/src/types.md index 4afb2088..7dafc147 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -14,12 +14,6 @@ Abstract supertype LinearMaps.LinearMap ``` -Unwrapping function - -```@docs -Base.parent -``` - ### `FunctionMap` Type for wrapping an arbitrary function that is supposed to implement the diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index f8cc42e6..cc9b146d 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -16,8 +16,9 @@ end abstract type LinearMap{T} end -const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}} -const RealOrComplex = Union{Real,Complex} +const MapOrVecOrMat{T} = Union{LinearMap{T}, AbstractVecOrMat{T}} +const MapOrMatrix{T} = Union{LinearMap{T}, AbstractMatrix{T}} +const RealOrComplex = Union{Real, Complex} Base.eltype(::LinearMap{T}) where {T} = T @@ -45,18 +46,10 @@ 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] -""" - parent(A::LinearMap) - -Return the underlying "parent map". This parent map is what was passed as an argument to -the specific `LinearMap` constructor, including implicit constructors and up to implicit -promotion to a `LinearMap` subtype. The fallback is to return the input itself. -""" -Base.parent(A::LinearMap) = A - # check dimension consistency for multiplication A*B _iscompatible((A, B)) = size(A, 2) == size(B, 1) function check_dim_mul(A, B) @@ -115,8 +108,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 +138,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 @@ -233,6 +227,8 @@ function _unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix, α, β return _generic_mapmat_mul!(y, A, x, α, β) end +const LinearMapTuple = Tuple{Vararg{LinearMap}} + include("left.jl") # left multiplication by a transpose or adjoint vector 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 @@ -257,15 +253,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 +272,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) diff --git a/src/blockmap.jl b/src/blockmap.jl index d31281b9..1e5b5a88 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -1,24 +1,29 @@ -struct BlockMap{T,As<:Tuple{Vararg{LinearMap}},Rs<:Tuple{Vararg{Int}},Rranges<:Tuple{Vararg{UnitRange{Int}}},Cranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T} +struct BlockMap{T, + As<:LinearMapTuple, + Rs<:Tuple{Vararg{Int}}, + Rranges<:Tuple{Vararg{UnitRange{Int}}}, + Cranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T} maps::As rows::Rs rowranges::Rranges colranges::Cranges - function BlockMap{T,R,S}(maps::R, rows::S) where {T, R<:Tuple{Vararg{LinearMap}}, S<:Tuple{Vararg{Int}}} - for n in eachindex(maps) - A = maps[n] - @assert promote_type(T, eltype(A)) == T "eltype $(eltype(A)) cannot be promoted to $T in BlockMap constructor" + function BlockMap{T,As,Rs}(maps::As, rows::Rs) where + {T, As<:LinearMapTuple, Rs<:Tuple{Vararg{Int}}} + for TA in Base.Generator(eltype, maps) + promote_type(T, TA) == T || + error("eltype $TA cannot be promoted to $T in BlockMap constructor") end rowranges, colranges = rowcolranges(maps, rows) - return new{T,R,S,typeof(rowranges),typeof(colranges)}(maps, rows, rowranges, colranges) + Rranges, Cranges = typeof(rowranges), typeof(colranges) + return new{T, As, Rs, Rranges, Cranges}(maps, rows, rowranges, colranges) end end -BlockMap{T}(maps::As, rows::S) where {T,As<:Tuple{Vararg{LinearMap}},S} = BlockMap{T,As,S}(maps, rows) +BlockMap{T}(maps::As, rows::Rs) where {T, As<:LinearMapTuple, Rs} = + BlockMap{T, As, Rs}(maps, rows) MulStyle(A::BlockMap) = MulStyle(A.maps...) -Base.parent(A::BlockMap) = A.maps - """ rowcolranges(maps, rows) @@ -27,24 +32,26 @@ map in `maps`, according to its position in a virtual matrix representation of t block linear map obtained from `hvcat(rows, maps...)`. """ function rowcolranges(maps, rows) - rowranges = () - colranges = () + rowranges = ntuple(n->1:0, Val(length(rows))) + colranges = ntuple(n->1:0, Val(length(maps))) mapind = 0 rowstart = 1 - for row in rows - xinds = vcat(1, map(a -> size(a, 2), maps[mapind+1:mapind+row])...) - cumsum!(xinds, xinds) + for (i, row) in enumerate(rows) mapind += 1 rowend = rowstart + size(maps[mapind], 1) - 1 - rowranges = (rowranges..., rowstart:rowend) - colranges = (colranges..., xinds[1]:xinds[2]-1) + rowranges = Base.setindex(rowranges, rowstart:rowend, i) + colstart = 1 + colend = size(maps[mapind], 2) + colranges = Base.setindex(colranges, colstart:colend, mapind) for colind in 2:row - mapind +=1 - colranges = (colranges..., xinds[colind]:xinds[colind+1]-1) + mapind += 1 + colstart = colend + 1 + colend += size(maps[mapind], 2) + colranges = Base.setindex(colranges, colstart:colend, mapind) end rowstart = rowend + 1 end - return rowranges::NTuple{length(rows), UnitRange{Int}}, colranges::NTuple{length(maps), UnitRange{Int}} + return rowranges, colranges end Base.size(A::BlockMap) = (last(last(A.rowranges)), last(last(A.colranges))) @@ -71,7 +78,7 @@ julia> L * ones(Int, 6) 6 ``` """ -function Base.hcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...) +function Base.hcat(As::Union{LinearMap, UniformScaling, AbstractVecOrMat}...) T = promote_type(map(eltype, As)...) nbc = length(As) @@ -83,8 +90,9 @@ function Base.hcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...) break end end - nrows == -1 && throw(ArgumentError("hcat of only UniformScaling objects cannot determine the linear map size")) - return BlockMap{T}(promote_to_lmaps(fill(nrows, nbc), 1, 1, As...), (nbc,)) + @assert nrows != -1 + # this should not happen, function should only be called with at least one LinearMap + return BlockMap{T}(promote_to_lmaps(ntuple(i->nrows, nbc), 1, 1, As...), (nbc,)) end ############ @@ -124,9 +132,10 @@ function Base.vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...) break end end - ncols == -1 && throw(ArgumentError("vcat of only UniformScaling objects cannot determine the linear map size")) - - return BlockMap{T}(promote_to_lmaps(fill(ncols, nbr), 1, 2, As...), ntuple(i->1, nbr)) + @assert ncols != -1 + # this should not happen, function should only be called with at least one LinearMap + rows = ntuple(i->1, nbr) + return BlockMap{T}(promote_to_lmaps(ntuple(i->ncols, nbr), 1, 2, As...), rows) end ############ @@ -160,10 +169,12 @@ julia> L * ones(Int, 6) """ Base.hvcat -function Base.hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...) +function Base.hvcat(rows::Tuple{Vararg{Int}}, + As::Union{LinearMap, UniformScaling, AbstractVecOrMat}...) nr = length(rows) T = promote_type(map(eltype, As)...) - sum(rows) == length(As) || throw(ArgumentError("mismatch between row sizes and number of arguments")) + sum(rows) == length(As) || + throw(ArgumentError("mismatch between row sizes and number of arguments")) n = fill(-1, length(As)) j = 0 for i in 1:nr # infer UniformScaling sizes from row counts, if possible: @@ -177,9 +188,7 @@ function Base.hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling end end if ni >= 0 - for k = 1:rows[i] - n[j+k] = ni - end + n[j .+ (1:rows[i])] .= ni end j += rows[i] end @@ -218,7 +227,8 @@ function check_dim(A, dim, n) end promote_to_lmaps_(n::Int, dim, A::AbstractMatrix) = (check_dim(A, dim, n); LinearMap(A)) -promote_to_lmaps_(n::Int, dim, A::AbstractVector) = (check_dim(A, dim, n); LinearMap(reshape(A, length(A), 1))) +promote_to_lmaps_(n::Int, dim, A::AbstractVector) = + (check_dim(A, dim, n); LinearMap(reshape(A, length(A), 1))) promote_to_lmaps_(n::Int, dim, J::UniformScaling) = UniformScalingMap(J.λ, n) promote_to_lmaps_(n::Int, dim, A::LinearMap) = (check_dim(A, dim, n); A) promote_to_lmaps(n, k, dim) = () @@ -236,16 +246,18 @@ function isblocksquare(A::BlockMap) return all(==(N), rows) end +symindex(i, N) = ((k, l) = divrem(i-1, N); return k + l * N + 1) + # the following rules are sufficient but not necessary function LinearAlgebra.issymmetric(A::BlockMap) isblocksquare(A) || return false N = length(A.rows) maps = A.maps - symindex = vec(permutedims(reshape(collect(1:N*N), N, N))) for i in 1:N*N - if (i == symindex[i] && !issymmetric(maps[i])) + isym = symindex(i, N) + if (i == isym && !issymmetric(maps[i])) return false - elseif (maps[i] != transpose(maps[symindex[i]])) + elseif (maps[i] != transpose(maps[isym])) return false end end @@ -257,11 +269,11 @@ function LinearAlgebra.ishermitian(A::BlockMap) isblocksquare(A) || return false N = length(A.rows) maps = A.maps - symindex = vec(permutedims(reshape(collect(1:N*N), N, N))) for i in 1:N*N - if (i == symindex[i] && !ishermitian(maps[i])) + isym = symindex(i, N) + if (i == isym && !ishermitian(maps[i])) return false - elseif (maps[i] != adjoint(maps[symindex[i]])) + elseif (maps[i] != adjoint(maps[isym])) return false end end @@ -272,7 +284,8 @@ end # comparison of BlockMap objects, sufficient but not necessary ############ -Base.:(==)(A::BlockMap, B::BlockMap) = (eltype(A) == eltype(B) && A.maps == B.maps && A.rows == B.rows) +Base.:(==)(A::BlockMap, B::BlockMap) = + (eltype(A) == eltype(B) && A.maps == B.maps && A.rows == B.rows) ############ # multiplication helper functions @@ -299,7 +312,7 @@ function ___blockmul!(y, A, x, α, β, ::Nothing) mapind += 1 _unsafe_mul!(yrow, maps[mapind], selectdim(x, 1, xinds[mapind]), α, β) for _ in 2:row - mapind +=1 + mapind += 1 _unsafe_mul!(yrow, maps[mapind], selectdim(x, 1, xinds[mapind]), α, true) end end @@ -312,15 +325,17 @@ function ___blockmul!(y, A, x, α, β, z) yrow = selectdim(y, 1, yi) zrow = selectdim(z, 1, yi) mapind += 1 + xrow = selectdim(x, 1, xinds[mapind]) if MulStyle(maps[mapind]) === ThreeArg() && !iszero(β) !isone(β) && rmul!(yrow, β) - muladd!(ThreeArg(), yrow, maps[mapind], selectdim(x, 1, xinds[mapind]), α, zrow) + muladd!(ThreeArg(), yrow, maps[mapind], xrow, α, zrow) else - _unsafe_mul!(yrow, maps[mapind], selectdim(x, 1, xinds[mapind]), α, β) + _unsafe_mul!(yrow, maps[mapind], xrow, α, β) end for _ in 2:row mapind +=1 - muladd!(MulStyle(maps[mapind]), yrow, maps[mapind], selectdim(x, 1, xinds[mapind]), α, zrow) + xrow = selectdim(x, 1, xinds[mapind]) + muladd!(MulStyle(maps[mapind]), yrow, maps[mapind], xrow, α, zrow) end end return y @@ -334,19 +349,21 @@ function _transblockmul!(y, A::BlockMap, x, α, β, transform) return rmul!(y, β) else # first block row (rowind = 1) of A, meaning first block column of A', fill all of y - xcol = selectdim(x, 1, first(xinds)) + xrow = selectdim(x, 1, first(xinds)) for rowind in 1:first(rows) - _unsafe_mul!(selectdim(y, 1, yinds[rowind]), transform(maps[rowind]), xcol, α, β) + 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 for (row, xi) in zip(Base.tail(rows), Base.tail(xinds)) - xcol = selectdim(x, 1, xi) + xrow = selectdim(x, 1, xi) for _ in 1:row mapind +=1 - _unsafe_mul!(selectdim(y, 1, yinds[mapind]), transform(maps[mapind]), xcol, α, true) + yrow = selectdim(y, 1, yinds[mapind]) + _unsafe_mul!(yrow, transform(maps[mapind]), xrow, α, true) end end end @@ -357,27 +374,26 @@ end # 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::BlockMap, x::$intype) + function _unsafe_mul!(y::$Out, A::BlockMap, x::$In) require_one_based_indexing(y, x) return _blockmul!(y, A, x, true, false) end - function _unsafe_mul!(y::$outtype, A::BlockMap, x::$intype, - α::Number, β::Number) + function _unsafe_mul!(y::$Out, A::BlockMap, x::$In, α::Number, β::Number) require_one_based_indexing(y, x) return _blockmul!(y, A, x, α, β) end end - for (maptype, transform) in ((:(TransposeMap{<:Any,<:BlockMap}), :transpose), (:(AdjointMap{<:Any,<:BlockMap}), :adjoint)) + for (MT, transform) in ((:TransposeMap, :transpose), (:AdjointMap, :adjoint)) @eval begin - function _unsafe_mul!(y::$outtype, wrapA::$maptype, x::$intype) + MapType = $MT{<:Any, <:BlockMap} + function _unsafe_mul!(y::$Out, wrapA::MapType, x::$In) require_one_based_indexing(y, x) return _transblockmul!(y, wrapA.lmap, x, true, false, $transform) end - function _unsafe_mul!(y::$outtype, wrapA::$maptype, x::$intype, - α::Number, β::Number) + function _unsafe_mul!(y::$Out, wrapA::MapType, x::$In, α::Number, β::Number) require_one_based_indexing(y, x) return _transblockmul!(y, wrapA.lmap, x, α, β, $transform) end @@ -388,15 +404,16 @@ end ############ # BlockDiagonalMap ############ - -struct BlockDiagonalMap{T,As<:Tuple{Vararg{LinearMap}},Ranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T} +struct BlockDiagonalMap{T, + As<:LinearMapTuple, + Ranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T} maps::As rowranges::Ranges colranges::Ranges - function BlockDiagonalMap{T,As}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} - for n in eachindex(maps) - A = maps[n] - @assert promote_type(T, eltype(A)) == T "eltype $(eltype(A)) cannot be promoted to $T in BlockDiagonalMap constructor" + function BlockDiagonalMap{T, As}(maps::As) where {T, As<:LinearMapTuple} + for TA in Base.Generator(eltype, maps) + promote_type(T, TA) == T || + error("eltype $TA cannot be promoted to $T in BlockDiagonalMap constructor") end # row ranges inds = vcat(1, size.(maps, 1)...) @@ -406,11 +423,11 @@ struct BlockDiagonalMap{T,As<:Tuple{Vararg{LinearMap}},Ranges<:Tuple{Vararg{Unit inds[2:end] .= size.(maps, 2) cumsum!(inds, inds) colranges = ntuple(i -> inds[i]:inds[i+1]-1, Val(length(maps))) - return new{T,As,typeof(rowranges)}(maps, rowranges, colranges) + return new{T, As, typeof(rowranges)}(maps, rowranges, colranges) end end -BlockDiagonalMap{T}(maps::As) where {T,As<:Tuple{Vararg{LinearMap}}} = +BlockDiagonalMap{T}(maps::As) where {T, As<:LinearMapTuple} = BlockDiagonalMap{T,As}(maps) BlockDiagonalMap(maps::LinearMap...) = BlockDiagonalMap{promote_type(map(eltype, maps)...)}(maps) @@ -418,21 +435,25 @@ BlockDiagonalMap(maps::LinearMap...) = # since the below methods are more specific than the Base method, # they would redefine Base/SparseArrays behavior for k in 1:8 # is 8 sufficient? - Is = ntuple(n->:($(Symbol(:A,n))::AbstractVecOrMat), Val(k-1)) + Is = ntuple(n->:($(Symbol(:A, n))::AbstractVecOrMat), Val(k-1)) # yields (:A1, :A2, :A3, ..., :A(k-1)) - L = :($(Symbol(:A,k))::LinearMap) + L = :($(Symbol(:A, k))::LinearMap) # yields :Ak - mapargs = ntuple(n -> :($(Symbol(:A,n))), Val(k-1)) + mapargs = ntuple(n ->:($(Symbol(:A, n))), Val(k-1)) # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1))) @eval begin - function SparseArrays.blockdiag($(Is...), $L, As::Union{LinearMap,AbstractVecOrMat}...) - return BlockDiagonalMap(convert_to_lmaps($(mapargs...))..., $(Symbol(:A,k)), convert_to_lmaps(As...)...) + function SparseArrays.blockdiag($(Is...), $L, As::MapOrVecOrMat...) + return BlockDiagonalMap(convert_to_lmaps($(mapargs...))..., + $(Symbol(:A, k)), + convert_to_lmaps(As...)...) end - function Base.cat($(Is...), $L, As::Union{LinearMap,AbstractVecOrMat}...; dims::Dims{2}) + function Base.cat($(Is...), $L, As::MapOrVecOrMat...; dims::Dims{2}) if dims == (1,2) - return BlockDiagonalMap(convert_to_lmaps($(mapargs...))..., $(Symbol(:A,k)), convert_to_lmaps(As...)...) + return BlockDiagonalMap(convert_to_lmaps($(mapargs...))..., + $(Symbol(:A, k)), + convert_to_lmaps(As...)...) else throw(ArgumentError("dims keyword in cat of LinearMaps must be (1,2)")) end @@ -445,7 +466,7 @@ end Construct a (lazy) representation of the diagonal concatenation of the arguments. To avoid fallback to the generic `SparseArrays.blockdiag`, there must be a `LinearMap` -object among the first 8 arguments. +object among the first 8 arguments. """ SparseArrays.blockdiag @@ -462,25 +483,25 @@ Base.size(A::BlockDiagonalMap) = (last(A.rowranges[end]), last(A.colranges[end]) MulStyle(A::BlockDiagonalMap) = MulStyle(A.maps...) -Base.parent(A::BlockDiagonalMap) = A.maps - LinearAlgebra.issymmetric(A::BlockDiagonalMap) = all(issymmetric, A.maps) LinearAlgebra.ishermitian(A::BlockDiagonalMap{<:Real}) = all(issymmetric, A.maps) LinearAlgebra.ishermitian(A::BlockDiagonalMap) = all(ishermitian, A.maps) -LinearAlgebra.adjoint(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}(map(adjoint, A.maps)) -LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}(map(transpose, A.maps)) +LinearAlgebra.adjoint(A::BlockDiagonalMap{T}) where {T} = + BlockDiagonalMap{T}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} = + BlockDiagonalMap{T}(map(transpose, A.maps)) -Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) = (eltype(A) == eltype(B) && A.maps == B.maps) +Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) = + (eltype(A) == eltype(B) && A.maps == B.maps) -for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) +for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) @eval begin - function _unsafe_mul!(y::$outtype, A::BlockDiagonalMap, x::$intype) + function _unsafe_mul!(y::$Out, A::BlockDiagonalMap, x::$In) require_one_based_indexing(y, x) return _blockscaling!(y, A, x, true, false) end - function _unsafe_mul!(y::$outtype, A::BlockDiagonalMap, x::$intype, - α::Number, β::Number) + function _unsafe_mul!(y::$Out, A::BlockDiagonalMap, x::$In, α::Number, β::Number) require_one_based_indexing(y, x) return _blockscaling!(y, A, x, α, β) end diff --git a/src/composition.jl b/src/composition.jl index 844382b8..9f27418c 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -1,23 +1,23 @@ -struct CompositeMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} +struct CompositeMap{T, As<:LinearMapTuple} <: LinearMap{T} maps::As # stored in order of application to vector function CompositeMap{T, As}(maps::As) where {T, As} N = length(maps) for n in 2:N check_dim_mul(maps[n], maps[n-1]) end - for n in eachindex(maps) - A = maps[n] - @assert promote_type(T, eltype(A)) == T "eltype $(eltype(A)) cannot be promoted to $T in CompositeMap constructor" + for TA in Base.Generator(eltype, maps) + # like lazy map; could use Base.Iterators.map in Julia >= 1.6 + promote_type(T, TA) == T || + error("eltype $TA cannot be promoted to $T in CompositeMap constructor") end new{T, As}(maps) end end -CompositeMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = CompositeMap{T, As}(maps) +CompositeMap{T}(maps::As) where {T, As<:LinearMapTuple} = CompositeMap{T, As}(maps) # 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 -Base.parent(A::CompositeMap) = A.maps # the following rules are sufficient but not necessary for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose), @@ -27,7 +27,8 @@ for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose), LinearAlgebra.$f(A::CompositeMap) = $_f(A.maps) $_f(maps::Tuple{}) = true $_f(maps::Tuple{<:LinearMap}) = $f(maps[1]) - $_f(maps::Tuple{Vararg{<:LinearMap}}) = maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps))) + $_f(maps::Tuple{Vararg{<:LinearMap}}) = + maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps))) # since the introduction of ScaledMap, the following cases cannot occur # function $_f(maps::Tuple{Vararg{<:LinearMap}}) # length(maps) >= 2 # if maps[1] isa UniformScalingMap{<:RealOrComplex} @@ -117,31 +118,37 @@ function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap) end # special transposition behavior -LinearAlgebra.transpose(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(transpose, reverse(A.maps))) -LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(adjoint, reverse(A.maps))) +LinearAlgebra.transpose(A::CompositeMap{T}) where {T} = + CompositeMap{T}(map(transpose, reverse(A.maps))) +LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} = + CompositeMap{T}(map(adjoint, reverse(A.maps))) # comparison of CompositeMap objects Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps) # multiplication with vectors -function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap}}, x::AbstractVector) +_unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) = + _compositemul!(y, A, x) +function _compositemul!(y::AbstractVecOrMat, + A::CompositeMap{<:Any,<:Tuple{LinearMap}}, + x::AbstractVector, + source = nothing, + dest = nothing) return _unsafe_mul!(y, A.maps[1], x) end -function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector) - _compositemul!(y, A, x, similar(y, size(A.maps[1], 1))) -end -function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector) - _compositemul!(y, A, x, similar(y, size(A.maps[1], 1)), similar(y, size(A.maps[2], 1))) -end - -function _compositemul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) - # no size checking, will be done by individual maps - _unsafe_mul!(z, A.maps[1], x) - _unsafe_mul!(y, A.maps[2], z) +function _compositemul!(y::AbstractVecOrMat, + A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, + source = similar(y, size(A.maps[1], 1)), + dest = nothing) + _unsafe_mul!(source, A.maps[1], x) + _unsafe_mul!(y, A.maps[2], source) return y end -function _compositemul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) - # no size checking, will be done by individual maps +function _compositemul!(y::AbstractVecOrMat, + A::CompositeMap, + x::AbstractVector, + source = similar(y, size(A.maps[1], 1)), + dest = similar(y, size(A.maps[2], 1))) N = length(A.maps) _unsafe_mul!(source, A.maps[1], x) for n in 2:N-1 diff --git a/src/conversion.jl b/src/conversion.jl index 21c0e932..ba80130f 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -47,8 +47,10 @@ SparseArrays.SparseMatrixCSC(A::LinearMap) = sparse(A) # special cases # ScaledMap -Base.Matrix{T}(A::ScaledMap{<:Any,<:Any,<:MatrixMap}) where {T} = convert(Matrix{T}, A.λ*A.lmap.lmap) -SparseArrays.sparse(A::ScaledMap{<:Any,<:Any,<:MatrixMap}) = convert(SparseMatrixCSC, A.λ*A.lmap.lmap) +Base.Matrix{T}(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) where {T} = + convert(Matrix{T}, A.λ * A.lmap.lmap) +SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) = + convert(SparseMatrixCSC, A.λ * A.lmap.lmap) # UniformScalingMap Base.Matrix{T}(J::UniformScalingMap) where {T} = Matrix{T}(J.λ*I, size(J)) @@ -62,25 +64,25 @@ SparseArrays.sparse(A::WrappedMap) = sparse(A.lmap) Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap) # TransposeMap & AdjointMap -for (TT, T) in ((AdjointMap, adjoint), (TransposeMap, transpose)) - @eval Base.convert(::Type{AbstractMatrix}, A::$TT) = $T(convert(AbstractMatrix, A.lmap)) - @eval SparseArrays.sparse(A::$TT) = $T(convert(SparseMatrixCSC, A.lmap)) +for (T, t) in ((AdjointMap, adjoint), (TransposeMap, transpose)) + @eval Base.convert(::Type{AbstractMatrix}, A::$T) = $t(convert(AbstractMatrix, A.lmap)) + @eval SparseArrays.sparse(A::$T) = $t(convert(SparseMatrixCSC, A.lmap)) end # LinearCombination -function Base.Matrix{T}(ΣA::LinearCombination{<:Any,<:Tuple{Vararg{MatrixMap}}}) where {T} +function Base.Matrix{T}(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}}) 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{MatrixMap}}}) +function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}}) maps = ΣA.maps mats = map(A->getfield(A, :lmap), maps) return convert(SparseMatrixCSC, sum(mats)) end # CompositeMap -function Base.Matrix{T}(AB::CompositeMap{<:Any,<:Tuple{MatrixMap,LinearMap}}) where {T} +function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, LinearMap}}) where {T} B, A = AB.maps require_one_based_indexing(B) Y = Matrix{eltype(AB)}(undef, size(AB)) @@ -92,39 +94,41 @@ end for ((TA, fieldA), (TB, fieldB)) in (((MatrixMap, :lmap), (MatrixMap, :lmap)), ((MatrixMap, :lmap), (UniformScalingMap, :λ)), ((UniformScalingMap, :λ), (MatrixMap, :lmap))) - @eval function Base.convert(::Type{AbstractMatrix}, AB::CompositeMap{<:Any,<:Tuple{$TB,$TA}}) + @eval function Base.convert(::Type{AbstractMatrix}, + AB::CompositeMap{<:Any,<:Tuple{$TB,$TA}}) B, A = AB.maps return A.$fieldA*B.$fieldB end end -function Base.Matrix{T}(AB::CompositeMap{<:Any,<:Tuple{MatrixMap,MatrixMap}}) where {T} +function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, MatrixMap}}) where {T} B, A = AB.maps return convert(Matrix{T}, A.lmap*B.lmap) end -function SparseArrays.sparse(AB::CompositeMap{<:Any,<:Tuple{MatrixMap,MatrixMap}}) +function SparseArrays.sparse(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, MatrixMap}}) B, A = AB.maps return convert(SparseMatrixCSC, A.lmap*B.lmap) end -function Base.Matrix{T}(λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}}) where {T} +function Base.Matrix{T}(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}}) where {T} A, J = λA.maps return convert(Matrix{T}, J.λ*A.lmap) end -function SparseArrays.sparse(λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}}) +function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}}) A, J = λA.maps return convert(SparseMatrixCSC, J.λ*A.lmap) end -function Base.Matrix{T}(Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}}) where {T} +function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}}) where {T} J, A = Aλ.maps return convert(Matrix{T}, A.lmap*J.λ) end -function SparseArrays.sparse(Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}}) +function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}}) J, A = Aλ.maps return convert(SparseMatrixCSC, 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)...) +Base.convert(::Type{AbstractMatrix}, A::BlockMap) = + hvcat(A.rows, convert.(AbstractMatrix, A.maps)...) function SparseArrays.sparse(A::BlockMap) return hvcat( A.rows, @@ -140,13 +144,15 @@ end # KroneckerMap & KroneckerSumMap Base.Matrix{T}(A::KroneckerMap) where {T} = kron(convert.(Matrix{T}, A.maps)...) -Base.convert(::Type{AbstractMatrix}, A::KroneckerMap) = kron(convert.(AbstractMatrix, A.maps)...) +Base.convert(::Type{AbstractMatrix}, A::KroneckerMap) = + kron(convert.(AbstractMatrix, A.maps)...) function SparseArrays.sparse(A::KroneckerMap) return kron( convert(SparseMatrixCSC, first(A.maps)), convert.(AbstractMatrix, Base.tail(A.maps))... ) end + function Base.Matrix{T}(L::KroneckerSumMap) where {T} A, B = L.maps IA = Diagonal(ones(Bool, size(A, 1))) diff --git a/src/functionmap.jl b/src/functionmap.jl index 8a7f97aa..6c80b3a9 100644 --- a/src/functionmap.jl +++ b/src/functionmap.jl @@ -17,13 +17,15 @@ function FunctionMap{T}(f::F1, fc::F2, M::Int, N::Int; end # additional constructors -FunctionMap{T}(f, M::Int; kwargs...) where {T} = FunctionMap{T}(f, nothing, M, M; kwargs...) -FunctionMap{T}(f, M::Int, N::Int; kwargs...) where {T} = FunctionMap{T}(f, nothing, M, N; kwargs...) -FunctionMap{T}(f, fc, M::Int; kwargs...) where {T} = FunctionMap{T}(f, fc, M, M; kwargs...) +FunctionMap{T}(f, M::Int; kwargs...) where {T} = + FunctionMap{T}(f, nothing, M, M; kwargs...) +FunctionMap{T}(f, M::Int, N::Int; kwargs...) where {T} = + FunctionMap{T}(f, nothing, M, N; kwargs...) +FunctionMap{T}(f, fc, M::Int; kwargs...) where {T} = + FunctionMap{T}(f, fc, M, M; kwargs...) # properties Base.size(A::FunctionMap) = (A.M, A.N) -Base.parent(A::FunctionMap) = (A.f, A.fc) LinearAlgebra.issymmetric(A::FunctionMap) = A._issymmetric LinearAlgebra.ishermitian(A::FunctionMap) = A._ishermitian LinearAlgebra.isposdef(A::FunctionMap) = A._isposdef @@ -31,6 +33,9 @@ ismutating(A::FunctionMap) = A._ismutating _ismutating(f) = first(methods(f)).nargs == 3 # multiplication with vector +const TransposeFunctionMap = TransposeMap{<:Any, <:FunctionMap} +const AdjointFunctionMap = AdjointMap{<:Any, <:FunctionMap} + function Base.:(*)(A::FunctionMap, x::AbstractVector) length(x) == A.N || throw(DimensionMismatch()) if ismutating(A) @@ -41,7 +46,7 @@ function Base.:(*)(A::FunctionMap, x::AbstractVector) end return y end -function Base.:(*)(A::AdjointMap{<:Any,<:FunctionMap}, x::AbstractVector) +function Base.:(*)(A::AdjointFunctionMap, x::AbstractVector) Afun = A.lmap ishermitian(Afun) && return Afun*x length(x) == size(A, 2) || throw(DimensionMismatch()) @@ -67,7 +72,7 @@ function Base.:(*)(A::AdjointMap{<:Any,<:FunctionMap}, x::AbstractVector) error("adjoint not implemented for $(A.lmap)") end end -function Base.:(*)(A::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector) +function Base.:(*)(A::TransposeFunctionMap, x::AbstractVector) Afun = A.lmap (issymmetric(Afun) || (isreal(A) && ishermitian(Afun))) && return Afun*x length(x) == size(A, 2) || throw(DimensionMismatch()) @@ -105,7 +110,7 @@ function _unsafe_mul!(y::AbstractVecOrMat, A::FunctionMap, x::AbstractVector) return y end -function _unsafe_mul!(y::AbstractVecOrMat, At::TransposeMap{<:Any,<:FunctionMap}, x::AbstractVector) +function _unsafe_mul!(y::AbstractVecOrMat, At::TransposeFunctionMap, x::AbstractVector) A = At.lmap (issymmetric(A) || (isreal(A) && ishermitian(A))) && return _unsafe_mul!(y, A, x) if A.fc !== nothing @@ -126,7 +131,7 @@ function _unsafe_mul!(y::AbstractVecOrMat, At::TransposeMap{<:Any,<:FunctionMap} end end -function _unsafe_mul!(y::AbstractVecOrMat, Ac::AdjointMap{<:Any,<:FunctionMap}, x::AbstractVector) +function _unsafe_mul!(y::AbstractVecOrMat, Ac::AdjointFunctionMap, x::AbstractVector) A = Ac.lmap ishermitian(A) && return _unsafe_mul!(y, A, x) if A.fc !== nothing diff --git a/src/kronecker.jl b/src/kronecker.jl index 5837bdc0..f6beefe5 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -1,16 +1,14 @@ -struct KroneckerMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} +struct KroneckerMap{T, As<:LinearMapTuple} <: LinearMap{T} maps::As - function KroneckerMap{T, As}(maps::As) where {T, As} - for n in eachindex(maps) - A = maps[n] - @assert promote_type(T, eltype(A)) == T "eltype $(eltype(A)) cannot be promoted to $T in KroneckerMap constructor" + function KroneckerMap{T}(maps::LinearMapTuple) where {T} + for TA in Base.Generator(eltype, maps) + promote_type(T, TA) == T || + error("eltype $TA cannot be promoted to $T in KroneckerMap constructor") end - return new{T,As}(maps) + return new{T,typeof(maps)}(maps) end end -KroneckerMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = KroneckerMap{T, As}(maps) - """ kron(A::LinearMap, B::LinearMap)::KroneckerMap kron(A, B, Cs...)::KroneckerMap @@ -46,24 +44,29 @@ julia> Matrix(Δ) 0 1 1 -4 ``` """ -Base.kron(A::LinearMap, B::LinearMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}((A, B)) -Base.kron(A::LinearMap, B::KroneckerMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A, B.maps...)) -Base.kron(A::KroneckerMap, B::LinearMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B)) -Base.kron(A::KroneckerMap, B::KroneckerMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B.maps...)) -Base.kron(A::LinearMap, B::LinearMap, Cs::LinearMap...) = KroneckerMap{promote_type(eltype(A), eltype(B), map(eltype, Cs)...)}(tuple(A, B, Cs...)) +Base.kron(A::LinearMap, B::LinearMap) = + KroneckerMap{promote_type(eltype(A), eltype(B))}((A, B)) +Base.kron(A::LinearMap, B::KroneckerMap) = + KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A, B.maps...)) +Base.kron(A::KroneckerMap, B::LinearMap) = + KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B)) +Base.kron(A::KroneckerMap, B::KroneckerMap) = + KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B.maps...)) +Base.kron(A::LinearMap, B::LinearMap, C::LinearMap, Ds::LinearMap...) = + kron(kron(A, B), C, Ds...) Base.kron(A::AbstractMatrix, B::LinearMap) = kron(LinearMap(A), B) Base.kron(A::LinearMap, B::AbstractMatrix) = kron(A, LinearMap(B)) # promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker product for k in 3:8 # is 8 sufficient? - Is = ntuple(n->:($(Symbol(:A,n))::AbstractMatrix), Val(k-1)) + Is = ntuple(n->:($(Symbol(:A, n))::AbstractMatrix), Val(k-1)) # yields (:A1, :A2, :A3, ..., :A(k-1)) - L = :($(Symbol(:A,k))::LinearMap) + L = :($(Symbol(:A, k))::LinearMap) # yields :Ak::LinearMap - mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1)) + mapargs = ntuple(n -> :(LinearMap($(Symbol(:A, n)))), Val(k-1)) # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1))) @eval Base.kron($(Is...), $L, As::MapOrMatrix...) = - kron($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...) + kron($(mapargs...), $(Symbol(:A, k)), convert_to_lmaps(As...)...) end struct KronPower{p} @@ -81,21 +84,19 @@ Construct a lazy representation of the `k`-th Kronecker power """ ⊗(k::Integer) = KronPower(k) -⊗(A, B, Cs...) = KroneckerMap{promote_type(eltype(A), eltype(B), map(eltype, Cs)...)}(convert_to_lmaps(A, B, Cs...)) +⊗(A, B, Cs...) = kron(convert_to_lmaps(A, B, Cs...)...) Base.:(^)(A::MapOrMatrix, ::KronPower{p}) where {p} = - ⊗(ntuple(n -> convert_to_lmaps_(A), Val(p))...) + kron(ntuple(n -> convert_to_lmaps_(A), Val(p))...) Base.size(A::KroneckerMap) = map(*, size.(A.maps)...) -Base.parent(A::KroneckerMap) = A.maps - LinearAlgebra.issymmetric(A::KroneckerMap) = all(issymmetric, A.maps) LinearAlgebra.ishermitian(A::KroneckerMap{<:Real}) = issymmetric(A) LinearAlgebra.ishermitian(A::KroneckerMap) = all(ishermitian, A.maps) -LinearAlgebra.adjoint(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(adjoint, A.maps)) -LinearAlgebra.transpose(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(transpose, A.maps)) +LinearAlgebra.adjoint(A::KroneckerMap) = KroneckerMap{eltype(A)}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::KroneckerMap) = KroneckerMap{eltype(A)}(map(transpose, A.maps)) Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) @@ -118,14 +119,14 @@ Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps end return y end -@inline function _kronmul!(y, B, X, At::Union{MatrixMap,UniformScalingMap}, T) +@inline function _kronmul!(y, B, X, At::Union{MatrixMap, UniformScalingMap}, T) na, ma = size(At) mb, nb = size(B) Y = reshape(y, (mb, ma)) - if nb*ma < mb*na + if nb*ma < mb*na _unsafe_mul!(Y, B, Matrix(X*At)) else - _unsafe_mul!(Y, Matrix(B*X), parent(At)) + _unsafe_mul!(Y, Matrix(B*X), _parent(At)) end return y end @@ -134,24 +135,30 @@ end # multiplication with vectors ################# -function _unsafe_mul!(y::AbstractVecOrMat, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} +const KroneckerMap2{T} = KroneckerMap{T, <:Tuple{LinearMap, LinearMap}} + +function _unsafe_mul!(y::AbstractVecOrMat, L::KroneckerMap2, x::AbstractVector) require_one_based_indexing(y) A, B = L.maps - X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); issymmetric=false, ishermitian=false, isposdef=false) - _kronmul!(y, B, X, transpose(A), T) + X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); + issymmetric = false, ishermitian = false, isposdef = false) + _kronmul!(y, B, X, transpose(A), eltype(L)) return y end -function _unsafe_mul!(y::AbstractVecOrMat, L::KroneckerMap{T}, x::AbstractVector) where {T} +function _unsafe_mul!(y::AbstractVecOrMat, L::KroneckerMap, x::AbstractVector) require_one_based_indexing(y) A = first(L.maps) B = kron(Base.tail(L.maps)...) - X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); issymmetric=false, ishermitian=false, isposdef=false) - _kronmul!(y, B, X, transpose(A), T) + X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); + issymmetric = false, ishermitian = false, isposdef = false) + _kronmul!(y, B, X, transpose(A), eltype(L)) return y end # mixed-product rule, prefer the right if possible # (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ) -function _unsafe_mul!(y::AbstractVecOrMat, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector) +function _unsafe_mul!(y::AbstractVecOrMat, + L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, + x::AbstractVector) require_one_based_indexing(y) B, A = L.maps if length(A.maps) == length(B.maps) && all(_iscompatible, zip(A.maps, B.maps)) @@ -162,8 +169,10 @@ function _unsafe_mul!(y::AbstractVecOrMat, L::CompositeMap{<:Any,<:Tuple{Kroneck return y end # mixed-product rule, prefer the right if possible -# (A₁ ⊗ B₁)*(A₂⊗B₂)*...*(Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ) -function _unsafe_mul!(y::AbstractVecOrMat, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T} +# (A₁⊗B₁) * (A₂⊗B₂) * ... * (Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ) +function _unsafe_mul!(y::AbstractVecOrMat, + L::CompositeMap{T, <:Tuple{Vararg{KroneckerMap2}}}, + x::AbstractVector) where {T} require_one_based_indexing(y) As = map(AB -> AB.maps[1], L.maps) Bs = map(AB -> AB.maps[2], L.maps) @@ -181,20 +190,20 @@ end ############### # KroneckerSumMap ############### -struct KroneckerSumMap{T, As<:Tuple{LinearMap,LinearMap}} <: LinearMap{T} +struct KroneckerSumMap{T, As<:Tuple{LinearMap, LinearMap}} <: LinearMap{T} maps::As - function KroneckerSumMap{T, As}(maps::As) where {T, As} - for n in eachindex(maps) - A = maps[n] - size(A, 1) == size(A, 2) || throw(ArgumentError("operators need to be square in Kronecker sums")) - @assert promote_type(T, eltype(A)) == T "eltype $(eltype(A)) cannot be promoted to $T in KroneckerSumMap constructor" + function KroneckerSumMap{T}(maps::Tuple{LinearMap,LinearMap}) where {T} + A1, A2 = maps + (size(A1, 1) == size(A1, 2) && size(A2, 1) == size(A2, 2)) || + throw(ArgumentError("operators need to be square in Kronecker sums")) + for TA in Base.Generator(eltype, maps) + promote_type(T, TA) == T || + error("eltype $TA cannot be promoted to $T in KroneckerSumMap constructor") end - return new{T,As}(maps) + return new{T, typeof(maps)}(maps) end end -KroneckerSumMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = KroneckerSumMap{T, As}(maps) - """ kronsum(A, B)::KroneckerSumMap kronsum(A, B, Cs...)::KroneckerSumMap @@ -246,20 +255,23 @@ where `A` can be a square `AbstractMatrix` or a `LinearMap`. ⊕(a, b, c...) = kronsum(a, b, c...) -Base.:(^)(A::MapOrMatrix, ::KronSumPower{p}) where {p} = kronsum(ntuple(n -> convert_to_lmaps_(A), Val(p))...) +Base.:(^)(A::MapOrMatrix, ::KronSumPower{p}) where {p} = + kronsum(ntuple(n->convert_to_lmaps_(A), Val(p))...) Base.size(A::KroneckerSumMap, i) = prod(size.(A.maps, i)) Base.size(A::KroneckerSumMap) = (size(A, 1), size(A, 2)) -Base.parent(A::KroneckerSumMap) = A.maps LinearAlgebra.issymmetric(A::KroneckerSumMap) = all(issymmetric, A.maps) LinearAlgebra.ishermitian(A::KroneckerSumMap{<:Real}) = all(issymmetric, A.maps) LinearAlgebra.ishermitian(A::KroneckerSumMap) = all(ishermitian, A.maps) -LinearAlgebra.adjoint(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(map(adjoint, A.maps)) -LinearAlgebra.transpose(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(map(transpose, A.maps)) +LinearAlgebra.adjoint(A::KroneckerSumMap) = + KroneckerSumMap{eltype(A)}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::KroneckerSumMap) = + KroneckerSumMap{eltype(A)}(map(transpose, A.maps)) -Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = (eltype(A) == eltype(B) && A.maps == B.maps) +Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = + (eltype(A) == eltype(B) && A.maps == B.maps) function _unsafe_mul!(y::AbstractVecOrMat, L::KroneckerSumMap, x::AbstractVector) A, B = L.maps diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 13fa9abc..ffaccd11 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -1,4 +1,4 @@ -struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} +struct LinearCombination{T, As<:LinearMapTuple} <: LinearMap{T} maps::As function LinearCombination{T, As}(maps::As) where {T, As} N = length(maps) @@ -18,10 +18,10 @@ MulStyle(A::LinearCombination) = MulStyle(A.maps...) # basic methods Base.size(A::LinearCombination) = size(A.maps[1]) -Base.parent(A::LinearCombination) = A.maps -LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps) # sufficient but not necessary -LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps) # sufficient but not necessary -LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps) # sufficient but not necessary +# following conditions are sufficient but not necessary +LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps) +LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps) +LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps) # adding linear maps """ @@ -59,24 +59,24 @@ end Base.:(-)(A₁::LinearMap, A₂::LinearMap) = +(A₁, -A₂) # comparison of LinearCombination objects, sufficient but not necessary -Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A) == eltype(B) && A.maps == B.maps) +Base.:(==)(A::LinearCombination, B::LinearCombination) = + (eltype(A) == eltype(B) && A.maps == B.maps) # special transposition behavior -LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A)}(map(transpose, A.maps)) -LinearAlgebra.adjoint(A::LinearCombination) = LinearCombination{eltype(A)}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::LinearCombination) = + LinearCombination{eltype(A)}(map(transpose, A.maps)) +LinearAlgebra.adjoint(A::LinearCombination) = + LinearCombination{eltype(A)}(map(adjoint, A.maps)) # 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::LinearCombination, x::$intype) - + function _unsafe_mul!(y::$Out, A::LinearCombination, x::$In) _unsafe_mul!(y, first(A.maps), x) _mul!(MulStyle(A), y, A, x, true) return y end - - function _unsafe_mul!(y::$outtype, A::LinearCombination, x::$intype, - α::Number, β::Number) + function _unsafe_mul!(y::$Out, A::LinearCombination, x::$In, α::Number, β::Number) if iszero(α) # trivial cases iszero(β) && return fill!(y, zero(eltype(y))) isone(β) && return y @@ -109,7 +109,7 @@ function _mul!(::ThreeArg, y, A::LinearCombination, x, α) return __mul!(y, Base.tail(A.maps), x, α, similar(y)) end -__mul!(y, As::Tuple{Vararg{LinearMap}}, x, α, z) = +__mul!(y, As::LinearMapTuple, x, α, z) = __mul!(__mul!(y, first(As), x, α, z), Base.tail(As), x, α, z) __mul!(y, A::Tuple{LinearMap}, x, α, z) = __mul!(y, first(A), x, α, z) __mul!(y, A::LinearMap, x, α, z) = muladd!(MulStyle(A), y, A, x, α, z) diff --git a/src/scaledmap.jl b/src/scaledmap.jl index ef80f59b..23cb7534 100644 --- a/src/scaledmap.jl +++ b/src/scaledmap.jl @@ -13,10 +13,9 @@ ScaledMap(λ::S, lmap::A) where {S<:RealOrComplex,A<:LinearMap} = # basic methods Base.size(A::ScaledMap) = size(A.lmap) -Base.parent(A::ScaledMap) = (A.λ, A.lmap) Base.isreal(A::ScaledMap) = isreal(A.λ) && isreal(A.lmap) LinearAlgebra.issymmetric(A::ScaledMap) = issymmetric(A.lmap) -LinearAlgebra.ishermitian(A::ScaledMap) = ishermitian(A.lmap) +LinearAlgebra.ishermitian(A::ScaledMap) = isreal(A.λ) && ishermitian(A.lmap) LinearAlgebra.isposdef(A::ScaledMap) = isposdef(A.λ) && isposdef(A.lmap) Base.transpose(A::ScaledMap) = A.λ * transpose(A.lmap) @@ -24,7 +23,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 +42,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 diff --git a/src/show.jl b/src/show.jl index 5b3ec08e..3b5ca5c2 100644 --- a/src/show.jl +++ b/src/show.jl @@ -21,13 +21,13 @@ function _show(io::IO, A::Union{CompositeMap,LinearCombination,KroneckerMap,Kron " with $n map" * (n>1 ? "s" : "") * ":\n" * print_maps(io, A.maps, i+2) end function _show(io::IO, A::Union{AdjointMap,TransposeMap,WrappedMap}, i) - " of\n" * map_show(io, parent(A), i+2) + " of\n" * map_show(io, A.lmap, i+2) end function _show(io::IO, A::BlockMap, i) nrows = length(A.rows) n = length(A.maps) - " with $n block map" * (n>1 ? "s" : "") * " in $nrows block row" * (nrows>1 ? "s" : "") * - '\n' * print_maps(io, A.maps, i+2) + " with $n block map" * (n>1 ? "s" : "") * + " in $nrows block row" * (nrows>1 ? "s" : "") * '\n' * print_maps(io, A.maps, i+2) end function _show(io::IO, A::BlockDiagonalMap, i) n = length(A.maps) @@ -45,7 +45,7 @@ function _show_typeof(A::LinearMap{T}) where {T} split(string(typeof(A)), '{')[1] * '{' * string(T) * '}' end -function print_maps(io::IO, maps::Tuple{Vararg{LinearMap}}, k) +function print_maps(io::IO, maps::LinearMapTuple, k) n = length(maps) str = "" if get(io, :limit, true) && n > 10 diff --git a/src/transpose.jl b/src/transpose.jl index e4e30361..001cfd09 100644 --- a/src/transpose.jl +++ b/src/transpose.jl @@ -11,8 +11,6 @@ MulStyle(A::Union{TransposeMap,AdjointMap}) = MulStyle(A.lmap) LinearAlgebra.transpose(A::TransposeMap) = A.lmap LinearAlgebra.adjoint(A::AdjointMap) = A.lmap -Base.parent(A::Union{AdjointMap,TransposeMap}) = A.lmap - """ transpose(A::LinearMap) @@ -52,46 +50,45 @@ 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)) +LinearAlgebra.transpose(A::ConjugateMap) = adjoint(A.lmap.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 @@ -99,7 +96,7 @@ end _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 diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 32d83081..350812bf 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -7,15 +7,17 @@ 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() # properties Base.size(A::UniformScalingMap) = (A.M, A.M) -Base.parent(A::UniformScalingMap) = A.λ +_parent(A::UniformScalingMap) = A.λ Base.isreal(A::UniformScalingMap) = isreal(A.λ) LinearAlgebra.issymmetric(::UniformScalingMap) = true LinearAlgebra.ishermitian(A::UniformScalingMap) = isreal(A) @@ -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)) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 9a5f77ef..4631c246 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -5,15 +5,15 @@ struct WrappedMap{T, A<:MapOrMatrix} <: LinearMap{T} _isposdef::Bool end function WrappedMap(lmap::MapOrMatrix{T}; - issymmetric::Bool = issymmetric(lmap), - ishermitian::Bool = ishermitian(lmap), - isposdef::Bool = isposdef(lmap)) where {T} + issymmetric::Bool = issymmetric(lmap), + ishermitian::Bool = ishermitian(lmap), + isposdef::Bool = isposdef(lmap)) where {T} WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef) end function WrappedMap{T}(lmap::MapOrMatrix; - issymmetric::Bool = issymmetric(lmap), - ishermitian::Bool = ishermitian(lmap), - isposdef::Bool = isposdef(lmap)) where {T} + issymmetric::Bool = issymmetric(lmap), + ishermitian::Bool = ishermitian(lmap), + isposdef::Bool = isposdef(lmap)) where {T} WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef) end @@ -22,53 +22,64 @@ const MatrixMap{T} = WrappedMap{T,<:AbstractMatrix} MulStyle(A::WrappedMap) = MulStyle(A.lmap) LinearAlgebra.transpose(A::MatrixMap{T}) where {T} = - WrappedMap{T}(transpose(A.lmap); issymmetric=A._issymmetric, ishermitian=A._ishermitian, isposdef=A._isposdef) + WrappedMap{T}(transpose(A.lmap); + issymmetric = A._issymmetric, + ishermitian = A._ishermitian, + isposdef = A._isposdef) LinearAlgebra.adjoint(A::MatrixMap{T}) where {T} = - WrappedMap{T}(adjoint(A.lmap); issymmetric=A._issymmetric, ishermitian=A._ishermitian, isposdef=A._isposdef) + WrappedMap{T}(adjoint(A.lmap); + issymmetric = A._issymmetric, + ishermitian = A._ishermitian, + isposdef = A._isposdef) Base.:(==)(A::MatrixMap, B::MatrixMap) = - (eltype(A)==eltype(B) && A.lmap==B.lmap && A._issymmetric==B._issymmetric && - A._ishermitian==B._ishermitian && A._isposdef==B._isposdef) + (eltype(A) == eltype(B) && A.lmap == B.lmap && A._issymmetric == B._issymmetric && + A._ishermitian == B._ishermitian && A._isposdef == B._isposdef) # properties Base.size(A::WrappedMap) = size(A.lmap) -Base.parent(A::WrappedMap) = A.lmap +_parent(A::WrappedMap) = A.lmap LinearAlgebra.issymmetric(A::WrappedMap) = A._issymmetric LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef # multiplication with vectors & matrices -for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) +for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) @eval begin - _unsafe_mul!(y::$outtype, A::WrappedMap, x::$intype) = _unsafe_mul!(y, A.lmap, x) - function _unsafe_mul!(y::$outtype, At::TransposeMap{<:Any,<:WrappedMap}, x::$intype) + _unsafe_mul!(y::$Out, A::WrappedMap, x::$In) = _unsafe_mul!(y, A.lmap, x) + function _unsafe_mul!(y::$Out, At::TransposeMap{<:Any,<:WrappedMap}, x::$In) A = At.lmap return (issymmetric(A) || (isreal(A) && ishermitian(A))) ? _unsafe_mul!(y, A.lmap, x) : _unsafe_mul!(y, transpose(A.lmap), x) end - function _unsafe_mul!(y::$outtype, Ac::AdjointMap{<:Any,<:WrappedMap}, x::$intype) + function _unsafe_mul!(y::$Out, Ac::AdjointMap{<:Any,<:WrappedMap}, x::$In) A = Ac.lmap - return ishermitian(A) ? _unsafe_mul!(y, A.lmap, x) : _unsafe_mul!(y, adjoint(A.lmap), x) + return ishermitian(A) ? + _unsafe_mul!(y, A.lmap, x) : + _unsafe_mul!(y, adjoint(A.lmap), x) end end end if VERSION ≥ v"1.3.0-alpha.115" - for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) + for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) @eval begin - function _unsafe_mul!(y::$outtype, A::WrappedMap, x::$intype, α::Number, β::Number) + function _unsafe_mul!(y::$Out, A::WrappedMap, x::$In, α::Number, β::Number) return _unsafe_mul!(y, A.lmap, x, α, β) end - function _unsafe_mul!(y::$outtype, At::TransposeMap{<:Any,<:WrappedMap}, x::$intype, α::Number, β::Number) + function _unsafe_mul!(y::$Out, At::TransposeMap{<:Any,<:WrappedMap}, x::$In, + α::Number, β::Number) A = At.lmap return (issymmetric(A) || (isreal(A) && ishermitian(A))) ? _unsafe_mul!(y, A.lmap, x, α, β) : _unsafe_mul!(y, transpose(A.lmap), x, α, β) end - function _unsafe_mul!(y::$outtype, Ac::AdjointMap{<:Any,<:WrappedMap}, x::$intype, α::Number, β::Number) + function _unsafe_mul!(y::$Out, Ac::AdjointMap{<:Any,<:WrappedMap}, x::$In, α::Number, β::Number) A = Ac.lmap - return ishermitian(A) ? _unsafe_mul!(y, A.lmap, x, α, β) : _unsafe_mul!(y, adjoint(A.lmap), x, α, β) + return ishermitian(A) ? + _unsafe_mul!(y, A.lmap, x, α, β) : + _unsafe_mul!(y, adjoint(A.lmap), x, α, β) end end end @@ -83,10 +94,10 @@ Base.:(-)(A₁::AbstractMatrix, A₂::LinearMap) = -(WrappedMap(A₁), A₂) """ *(A::LinearMap, X::AbstractMatrix)::CompositeMap -Return the `CompositeMap` `A*LinearMap(X)`, interpreting the matrix `X` as a linear operator, -rather than a collection of column vectors. To compute the action of `A` on each column of `X`, -call `Matrix(A*X)` or use the in-place multiplication `mul!(Y, A, X[, α, β])` with an appropriately sized, -preallocated matrix `Y`. +Return the `CompositeMap` `A*LinearMap(X)`, interpreting the matrix `X` as a linear +operator, rather than a collection of column vectors. To compute the action of `A` on each +column of `X`, call `Matrix(A*X)` or use the in-place multiplication `mul!(Y, A, X[, α, β])` +with an appropriately sized, preallocated matrix `Y`. ## Examples ```jldoctest; setup=(using LinearMaps) @@ -101,8 +112,8 @@ Base.:(*)(A₁::LinearMap, A₂::AbstractMatrix) = *(A₁, WrappedMap(A₂)) """ *(X::AbstractMatrix, A::LinearMap)::CompositeMap -Return the `CompositeMap` `LinearMap(X)*A`, interpreting the matrix `X` as a linear operator. -To compute the right-action of `A` on each row of `X`, call `Matrix(X*A)`. +Return the `CompositeMap` `LinearMap(X)*A`, interpreting the matrix `X` as a linear +operator. To compute the right-action of `A` on each row of `X`, call `Matrix(X*A)`. ## Examples ```jldoctest; setup=(using LinearMaps) diff --git a/test/blockmap.jl b/test/blockmap.jl index 79b37aca..10fc1cd7 100644 --- a/test/blockmap.jl +++ b/test/blockmap.jl @@ -7,12 +7,11 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive A12 = rand(elty, 10, n2) v = rand(elty, 10) L = @inferred hcat(LinearMap(A11), LinearMap(A12)) - @test parent(L) == (LinearMap(A11), LinearMap(A12)) @test occursin("10×$(10+n2) LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L)) @test @inferred(LinearMaps.MulStyle(L)) === matrixstyle @test L isa LinearMaps.BlockMap{elty} if elty <: Complex - @test_throws AssertionError LinearMaps.BlockMap{Float64}((LinearMap(A11), LinearMap(A12)), (2,)) + @test_throws ErrorException LinearMaps.BlockMap{Float64}((LinearMap(A11), LinearMap(A12)), (2,)) end A = [A11 A12] x = rand(10+n2) @@ -193,7 +192,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse: if elty <: Complex - @test_throws AssertionError LinearMaps.BlockDiagonalMap{Float64}((L1, L2, L3, L2, L1)) + @test_throws ErrorException LinearMaps.BlockDiagonalMap{Float64}((L1, L2, L3, L2, L1)) end Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...)) @test (@which blockdiag(sparse.((M1, M2, M3, M2, M1))...)).module != LinearMaps @@ -201,7 +200,6 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive x = randn(elty, size(Md, 2)) Bd = @inferred blockdiag(L1, L2, L3, L2, L1) @test @inferred(LinearMaps.MulStyle(Bd)) === matrixstyle - @test parent(Bd) == (L1, L2, L3, L2, L1) @test occursin("25×39 LinearMaps.BlockDiagonalMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), Bd)) @test Matrix(Bd) == Md @test convert(AbstractMatrix, Bd) isa SparseMatrixCSC diff --git a/test/composition.jl b/test/composition.jl index 3a538254..cb2d42a1 100644 --- a/test/composition.jl +++ b/test/composition.jl @@ -6,7 +6,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays FCM = LinearMaps.CompositeMap{ComplexF64}((FC,)) L = LowerTriangular(ones(10,10)) @test_throws DimensionMismatch F * LinearMap(rand(2,2)) - @test_throws AssertionError LinearMaps.CompositeMap{Float64}((FC, LinearMap(rand(10,10)))) + @test_throws ErrorException LinearMaps.CompositeMap{Float64}((FC, LinearMap(rand(10,10)))) A = 2 * rand(ComplexF64, (10, 10)) .- 1 B = rand(size(A)...) H = LinearMap(Hermitian(A'A)) @@ -20,7 +20,6 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays @test @inferred (A * F) * v == @inferred A * (F * v) @test @inferred A * (F * F) * v == @inferred A * (F * (F * v)) F2 = F*F - @test parent(F2) == (F, F) FC2 = FC*FC F4 = FC2 * F2 @test occursin("10×10 LinearMaps.CompositeMap{$(eltype(F4))}", sprint((t, s) -> show(t, "text/plain", s), F4)) diff --git a/test/functionmap.jl b/test/functionmap.jl index 05ffd9fc..e73a31f8 100644 --- a/test/functionmap.jl +++ b/test/functionmap.jl @@ -18,7 +18,6 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools U = Matrix(MyFT) # will be a unitary matrix @test @inferred U'U ≈ Matrix{eltype(U)}(I, N, N) @test occursin("$N×$N LinearMaps.FunctionMap{$(eltype(MyFT))}", sprint((t, s) -> show(t, "text/plain", s), MyFT)) - @test parent(LinearMap{ComplexF64}(myft, N)) === (myft, nothing) CS = @inferred LinearMap(cumsum, 2) @test size(CS) == (2, 2) diff --git a/test/kronecker.jl b/test/kronecker.jl index f9031a62..7acd6e78 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -8,8 +8,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays LA = LinearMap(A) LB = LinearMap(B) LK = @inferred kron(LA, LB) - @test parent(LK) == (LA, LB) - @test_throws AssertionError LinearMaps.KroneckerMap{Float64}((LA, LB)) + @test_throws ErrorException LinearMaps.KroneckerMap{Float64}((LA, LB)) @test occursin("6×6 LinearMaps.KroneckerMap{$(eltype(LK))}", sprint((t, s) -> show(t, "text/plain", s), LK)) @test @inferred size(LK) == size(K) @test LinearMaps.MulStyle(LK) === LinearMaps.ThreeArg() @@ -24,7 +23,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays @test Matrix(transform(LinearMap(LK))) ≈ transform(Matrix(LK)) ≈ transform(K) end @test kron(LA, kron(LA, B)) == kron(LA, LA, LB) - @test kron(kron(LA, LB), kron(LA, LB)) == kron(LA, LB, LA, LB) + @test kron(kron(LA, LB), kron(LA, LB)) == kron(LA, LB, LA, LB) == ⊗(LA, LB, LA, LB) @test kron(A, A, A) ≈ Matrix(@inferred kron(LA, LA, LA)) ≈ Matrix(@inferred LA^⊗(3)) ≈ Matrix(@inferred A^⊗(3)) LAs = LinearMap(sparse(A)) K = @inferred kron(A, A, A, LAs) @@ -70,7 +69,6 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays LA = LinearMap(A) LB = LinearMap(B) KS = @inferred kronsum(LA, B) - @test parent(KS) == (LA, LB) @test occursin("6×6 LinearMaps.KroneckerSumMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), KS)) @test_throws ArgumentError kronsum(LA, [B B]) # non-square map KSmat = kron(A, Matrix(I, 2, 2)) + kron(Matrix(I, 3, 3), B) diff --git a/test/linearcombination.jl b/test/linearcombination.jl index 9d170a87..b1c0b2f2 100644 --- a/test/linearcombination.jl +++ b/test/linearcombination.jl @@ -15,7 +15,6 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools @test_throws AssertionError LinearMaps.LinearCombination{Float64}((CS!, CS!)) @test occursin("10×10 LinearMaps.LinearCombination{$(eltype(L))}", sprint((t, s) -> show(t, "text/plain", s), L)) @test occursin("10×10 LinearMaps.LinearCombination{$(eltype(L))}", sprint((t, s) -> show(t, "text/plain", s), L+CS!)) - @test parent(L) == ntuple(_ -> CS!, 10) @test mul!(u, L, v) ≈ n * cumsum(v) b = @benchmarkable mul!($u, $L, $v, 2, 2) @test run(b, samples=5).allocs <= 1 diff --git a/test/linearmaps.jl b/test/linearmaps.jl index 9801d61c..7164dbef 100644 --- a/test/linearmaps.jl +++ b/test/linearmaps.jl @@ -48,7 +48,6 @@ LinearAlgebra.mul!(y::AbstractVector, A::Union{SimpleFunctionMap,SimpleComplexFu @testset "new LinearMap type" begin F = SimpleFunctionMap(cumsum, 10) - @test parent(F) === F FC = SimpleComplexFunctionMap(cumsum, 10) @test @inferred ndims(F) == 2 @test @inferred size(F, 1) == 10 diff --git a/test/scaledmap.jl b/test/scaledmap.jl index 2948b908..1105093c 100644 --- a/test/scaledmap.jl +++ b/test/scaledmap.jl @@ -9,7 +9,6 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools α = float(π) B = @inferred α * A @test occursin("7×7 LinearMaps.ScaledMap{Float64} with scale: $α", sprint((t, s) -> show(t, "text/plain", s), B)) - @test parent(B) == (α, A) x = rand(N) @test @inferred size(B) == size(A) diff --git a/test/transpose.jl b/test/transpose.jl index 08e25f87..ef240505 100644 --- a/test/transpose.jl +++ b/test/transpose.jl @@ -10,8 +10,6 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays end @test occursin("10×10 LinearMaps.TransposeMap{$(eltype(CS))}", sprint((t, s) -> show(t, "text/plain", s), transpose(CS))) @test occursin("10×10 LinearMaps.AdjointMap{$(eltype(CS))}", sprint((t, s) -> show(t, "text/plain", s), adjoint(CS))) - @test parent(transpose(CS)) === CS - @test parent(adjoint(CS)) === CS @test !(transpose(CS) == adjoint(CS)) @test !(adjoint(CS) == transpose(CS)) M = Matrix(CS) @@ -29,6 +27,8 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays @test mul!(copy(W), transform(CS), X, α, β) ≈ transform(M)*X*α + W*β end for transform1 in (adjoint, transpose), transform2 in (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 diff --git a/test/uniformscalingmap.jl b/test/uniformscalingmap.jl index 068b2024..8b158701 100644 --- a/test/uniformscalingmap.jl +++ b/test/uniformscalingmap.jl @@ -11,7 +11,7 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools w = similar(v) Id = @inferred LinearMap(I, 10) @test occursin("10×10 LinearMaps.UniformScalingMap{Bool}", sprint((t, s) -> show(t, "text/plain", s), Id)) - @test parent(Id) == true + @test LinearMaps._parent(Id) == true @test_throws ErrorException LinearMaps.UniformScalingMap(1, 10, 20) @test_throws ErrorException LinearMaps.UniformScalingMap(1, (10, 20)) @test size(Id) == (10, 10) diff --git a/test/wrappedmap.jl b/test/wrappedmap.jl index 5cb26d5c..485ea99b 100644 --- a/test/wrappedmap.jl +++ b/test/wrappedmap.jl @@ -7,7 +7,7 @@ using Test, LinearMaps, LinearAlgebra SB = B'B + I L = @inferred LinearMap{Float64}(A) @test occursin("10×20 LinearMaps.WrappedMap{Float64}", sprint((t, s) -> show(t, "text/plain", s), L)) - @test parent(L) === A + @test LinearMaps._parent(L) === A MA = @inferred LinearMap(SA) MB = @inferred LinearMap(SB) @test eltype(Matrix{Complex{Float32}}(LinearMap(A))) <: Complex