From 2c58a409d325e43a683b8a9eea80c807cca3b96f Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Sat, 5 Dec 2020 01:03:43 +0100 Subject: [PATCH 01/15] review part 1: some spaces and linebreaks, in/outtype -> In/Out for brevity --- src/LinearMaps.jl | 42 +++++++++++++------------ src/left.jl | 2 ++ src/scaledmap.jl | 8 ++--- src/transpose.jl | 66 +++++++++++++++++++--------------------- src/uniformscalingmap.jl | 64 ++++++++++++++++---------------------- src/wrappedmap.jl | 63 ++++++++++++++++++++++---------------- 6 files changed, 124 insertions(+), 121 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index f8cc42e6..f0206e66 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -16,8 +16,8 @@ end abstract type LinearMap{T} end -const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}} -const RealOrComplex = Union{Real,Complex} +const MapOrMatrix{T} = Union{LinearMap{T}, AbstractMatrix{T}} +const RealOrComplex = Union{Real, Complex} Base.eltype(::LinearMap{T}) where {T} = T @@ -45,7 +45,8 @@ LinearAlgebra.ishermitian(::LinearMap) = false # default assumptions LinearAlgebra.isposdef(::LinearMap) = false # default assumptions Base.ndims(::LinearMap) = 2 -Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions")) +Base.size(A::LinearMap, n) = + (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions")) Base.length(A::LinearMap) = size(A)[1] * size(A)[2] """ @@ -115,8 +116,9 @@ end mul!(Y::AbstractVecOrMat, A::LinearMap, B::AbstractVector) -> Y mul!(Y::AbstractMatrix, A::LinearMap, B::AbstractMatrix) -> Y -Calculates the action of the linear map `A` on the vector or matrix `B` and stores the result in `Y`, -overwriting the existing value of `Y`. Note that `Y` must not be aliased with either `A` or `B`. +Calculates the action of the linear map `A` on the vector or matrix `B` and stores the +result in `Y`, overwriting the existing value of `Y`. Note that `Y` must not be aliased +with either `A` or `B`. ## Examples ```jldoctest; setup=(using LinearAlgebra, LinearMaps) @@ -144,26 +146,26 @@ end mul!(C::AbstractVecOrMat, A::LinearMap, B::AbstractVector, α, β) -> C mul!(C::AbstractMatrix, A::LinearMap, B::AbstractMatrix, α, β) -> C -Combined inplace multiply-add ``A B α + C β``. The result is stored in `C` by overwriting it. -Note that `C` must not be aliased with either `A` or `B`. +Combined inplace multiply-add ``A B α + C β``. The result is stored in `C` by overwriting +it. Note that `C` must not be aliased with either `A` or `B`. ## Examples ```jldoctest; setup=(using LinearAlgebra, LinearMaps) julia> A=LinearMap([1.0 2.0; 3.0 4.0]); B=[1.0, 1.0]; C=[1.0, 3.0]; - + julia> mul!(C, A, B, 100.0, 10.0) === C true - + julia> C 2-element Array{Float64,1}: 310.0 730.0 julia> A=LinearMap([1.0 2.0; 3.0 4.0]); B=[1.0 1.0; 1.0 1.0]; C=[1.0 2.0; 3.0 4.0]; - + julia> mul!(C, A, B, 100.0, 10.0) === C true - + julia> C 2×2 Array{Float64,2}: 310.0 320.0 @@ -257,15 +259,17 @@ with the purpose of redefining its properties via the keyword arguments `kwargs` a `UniformScaling` object `J` with specified (square) dimension `M`; or from a function or callable object `f`. In the latter case, one also needs to specify the size of the equivalent matrix representation `(M, N)`, i.e., for functions `f` acting -on length `N` vectors and producing length `M` vectors (with default value `N=M`). Preferably, -also the `eltype` `T` of the corresponding matrix representation needs to be specified, i.e. -whether the action of `f` on a vector will be similar to, e.g., multiplying by numbers of type `T`. -If not specified, the devault value `T=Float64` will be assumed. Optionally, a corresponding -function `fc` can be specified that implements the adjoint (=transpose in the real case) of `f`. +on length `N` vectors and producing length `M` vectors (with default value `N=M`). +Preferably, also the `eltype` `T` of the corresponding matrix representation needs to be +specified, i.e. whether the action of `f` on a vector will be similar to, e.g., multiplying +by numbers of type `T`. If not specified, the devault value `T=Float64` will be assumed. +Optionally, a corresponding function `fc` can be specified that implements the adjoint +(=transpose in the real case) of `f`. The keyword arguments and their default values for the function-based constructor are: * `issymmetric::Bool = false` : whether `A` or `f` acts as a symmetric matrix -* `ishermitian::Bool = issymmetric & T<:Real` : whether `A` or `f` acts as a Hermitian matrix +* `ishermitian::Bool = issymmetric & T<:Real` : whether `A` or `f` acts as a Hermitian + matrix * `isposdef::Bool = false` : whether `A` or `f` acts as a positive definite matrix. For existing linear maps or matrices `A`, the default values will be taken by calling `issymmetric`, `ishermitian` and `isposdef` on the existing object `A`. @@ -274,8 +278,8 @@ For the function-based constructor, there is one more keyword argument: * `ismutating::Bool` : flags whether the function acts as a mutating matrix multiplication `f(y,x)` where the result vector `y` is the first argument (in case of `true`), or as a normal matrix multiplication that is called as `y=f(x)` (in case of `false`). - The default value is guessed by looking at the number of arguments of the first occurrence - of `f` in the method table. + The default value is guessed by looking at the number of arguments of the first + occurrence of `f` in the method table. """ LinearMap(A::MapOrMatrix; kwargs...) = WrappedMap(A; kwargs...) LinearMap(J::UniformScaling, M::Int) = UniformScalingMap(J.λ, M) diff --git a/src/left.jl b/src/left.jl index 4f6bddff..9ad42d57 100644 --- a/src/left.jl +++ b/src/left.jl @@ -50,6 +50,7 @@ end function mul!(x::AbstractMatrix, y::AdjointAbsVecOrMat, A::LinearMap, α::Number, β::Number) check_dim_mul(x, y, A) _unsafe_mul!(x', conj(α)*A', y', true, conj(β)) + # why not?: _unsafe_mul!(x', A', y', conj(α), conj(β)) return x end @@ -62,5 +63,6 @@ end function mul!(x::AbstractMatrix, y::TransposeAbsVecOrMat, A::LinearMap, α::Number, β::Number) check_dim_mul(x, y, A) _unsafe_mul!(transpose(x), α*transpose(A), transpose(y), true, β) + # why not?: _unsafe_mul!(transpose(x), transpose(A), transpose(y), α, β) return x end diff --git a/src/scaledmap.jl b/src/scaledmap.jl index ef80f59b..167c6e7e 100644 --- a/src/scaledmap.jl +++ b/src/scaledmap.jl @@ -24,7 +24,7 @@ Base.adjoint(A::ScaledMap) = conj(A.λ) * adjoint(A.lmap) # comparison (sufficient, not necessary) Base.:(==)(A::ScaledMap, B::ScaledMap) = - (eltype(A) == eltype(B) && A.lmap == B.lmap) && A.λ == B.λ + eltype(A) == eltype(B) && A.lmap == B.lmap && A.λ == B.λ # scalar multiplication and division Base.:(*)(α::RealOrComplex, A::LinearMap) = ScaledMap(α, A) @@ -43,12 +43,12 @@ Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B) Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ # multiplication with vectors/matrices -for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) +for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) @eval begin - function _unsafe_mul!(y::$outtype, A::ScaledMap, x::$intype) + function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In) return _unsafe_mul!(y, A.lmap, x, A.λ, false) end - function _unsafe_mul!(y::$outtype, A::ScaledMap, x::$intype, α::Number, β::Number) + function _unsafe_mul!(y::$Out, A::ScaledMap, x::$In, α::Number, β::Number) return _unsafe_mul!(y, A.lmap, x, A.λ * α, β) end end diff --git a/src/transpose.jl b/src/transpose.jl index e4e30361..3e840f73 100644 --- a/src/transpose.jl +++ b/src/transpose.jl @@ -52,46 +52,44 @@ Base.:(==)(A::LinearMap, B::AdjointMap) = ishermitian(A) && B.lmap == A # multiplication with vector/matrices # # TransposeMap -function _unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector) - issymmetric(A.lmap) ? _unsafe_mul!(y, A.lmap, x) : error("transpose not implemented for $(A.lmap)") -end -function _unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix) - issymmetric(A.lmap) ? _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x) -end -function _unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector, α::Number, β::Number) - issymmetric(A.lmap) ? _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β) -end -function _unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix, α::Number, β::Number) - issymmetric(A.lmap) ? _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β) -end +_unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector) = + issymmetric(A.lmap) ? + _unsafe_mul!(y, A.lmap, x) : error("transpose not implemented for $(A.lmap)") +_unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix) = + issymmetric(A.lmap) ? + _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x) +_unsafe_mul!(y::AbstractVecOrMat, A::TransposeMap, x::AbstractVector, α::Number, β::Number)= + issymmetric(A.lmap) ? + _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β) +_unsafe_mul!(y::AbstractMatrix, A::TransposeMap, x::AbstractMatrix, α::Number, β::Number) = + issymmetric(A.lmap) ? + _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β) # # AdjointMap -function _unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector) - ishermitian(A.lmap) ? _unsafe_mul!(y, A.lmap, x) : error("adjoint not implemented for $(A.lmap)") -end -function _unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix) - ishermitian(A.lmap) ? _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x) -end -function _unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector, α::Number, β::Number) - ishermitian(A.lmap) ? _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β) -end -function _unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix, α::Number, β::Number) - ishermitian(A.lmap) ? _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β) -end +_unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector) = + ishermitian(A.lmap) ? + _unsafe_mul!(y, A.lmap, x) : error("adjoint not implemented for $(A.lmap)") +_unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix) = + ishermitian(A.lmap) ? + _unsafe_mul!(y, A.lmap, x) : _generic_mapmat_mul!(y, A, x) +_unsafe_mul!(y::AbstractVecOrMat, A::AdjointMap, x::AbstractVector, α::Number, β::Number) = + ishermitian(A.lmap) ? + _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapvec_mul!(y, A, x, α, β) +_unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix, α::Number, β::Number) = + ishermitian(A.lmap) ? + _unsafe_mul!(y, A.lmap, x, α, β) : _generic_mapmat_mul!(y, A, x, α, β) + # # ConjugateMap -for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) +const ConjugateMap = AdjointMap{<:Any, <:TransposeMap} +# canonical order of adjoint followed by transpose +LinearAlgebra.transpose(A::AdjointMap) = adjoint(transpose(A.lmap)) +for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) @eval begin - function _unsafe_mul!(y::$outtype, Ac::AdjointMap{<:Any,<:TransposeMap}, x::$intype) + function _unsafe_mul!(y::$Out, Ac::ConjugateMap, x::$In) return _conjmul!(y, Ac.lmap.lmap, x) end - function _unsafe_mul!(y::$outtype, Ac::AdjointMap{<:Any,<:TransposeMap}, x::$intype, α::Number, β::Number) + function _unsafe_mul!(y::$Out, Ac::ConjugateMap, x::$In, α::Number, β::Number) return _conjmul!(y, Ac.lmap.lmap, x, α, β) end - function _unsafe_mul!(y::$outtype, At::TransposeMap{<:Any,<:AdjointMap}, x::$intype) - return _conjmul!(y, At.lmap.lmap, x) - end - function _unsafe_mul!(y::$outtype, At::TransposeMap{<:Any,<:AdjointMap}, x::$intype, α::Number, β::Number) - return _conjmul!(y, At.lmap.lmap, x, α, β) - end end end @@ -99,7 +97,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..462920ab 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -7,9 +7,11 @@ struct UniformScalingMap{T} <: LinearMap{T} end end UniformScalingMap(λ::Number, M::Int, N::Int) = - (M == N ? UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square")) + (M == N ? + UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square")) UniformScalingMap(λ::T, sz::Dims{2}) where {T} = - (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square")) + (sz[1] == sz[2] ? + UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square")) MulStyle(::UniformScalingMap) = FiveArg() @@ -42,13 +44,13 @@ Base.:(*)(A::UniformScalingMap, x::AbstractVector) = length(x) == A.M ? A.λ * x : throw(DimensionMismatch("*")) # multiplication with vector/matrix -for (intype, outtype) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) +for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) @eval begin - function _unsafe_mul!(y::$outtype, J::UniformScalingMap, x::$intype) + function _unsafe_mul!(y::$Out, J::UniformScalingMap, x::$In) _scaling!(y, J.λ, x, true, false) return y end - function _unsafe_mul!(y::$outtype, J::UniformScalingMap, x::$intype, + function _unsafe_mul!(y::$Out, J::UniformScalingMap, x::$In, α::Number, β::Number) _scaling!(y, J.λ, x, α, β) return y @@ -58,41 +60,27 @@ end function _scaling!(y, λ, x, α, β) if (iszero(α) || iszero(λ)) - iszero(β) && (fill!(y, zero(eltype(y))); return y) + iszero(β) && return fill!(y, zero(eltype(y))) isone(β) && return y - # β != 0, 1 - rmul!(y, β) - return y + return rmul!(y, β) + elseif isone(α) && isone(λ) + iszero(β) && return copyto!(y, x) + isone(β) && return y .+= x + return y .= y .* β .+ x elseif isone(α) - if iszero(β) - isone(λ) && return copyto!(y, x) - y .= λ .* x - return y - elseif isone(β) - isone(λ) && return y .+= x - y .+= λ .* x - return y - else # β != 0, 1 - isone(λ) && (axpby!(one(eltype(x)), x, β, y); return y) - y .= y .* β .+ λ .* x - return y - end - else # α != 0, 1 - if iszero(β) - isone(λ) && (y .= x .* α; return y) - y .= λ .* x .* α - return y - elseif isone(β) - isone(λ) && (axpby!(α, x, β, y); return y) - y .+= λ .* x .* α - return y - else # β != 0, 1 - isone(λ) && (y .= y .* β .+ x .* α; return y) - y .= y .* β .+ λ .* x .* α - return y - end # β-cases - end # α-cases -end # function _scaling! + iszero(β) && return y .= λ .* x + isone(β) && return y .+= λ .* x + return y .= y .* β .+ λ .* x + elseif isone(λ) + iszero(β) && return y .= x .* α + isone(β) && return y .+= x .* α + return y .= y .* β .+ x .* α + else + iszero(β) && return y .= λ .* x .* α + isone(β) && return y .+= λ .* x .* α + return y .= y .* β .+ λ .* x .* α + end +end # combine LinearMap and UniformScaling objects in linear combinations Base.:(+)(A₁::LinearMap, A₂::UniformScaling) = A₁ + UniformScalingMap(A₂.λ, size(A₁, 1)) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 9a5f77ef..d149af46 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,13 +22,19 @@ 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) @@ -38,37 +44,42 @@ 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) From b9036aa7ceee4ceffbb8fa4e1de7b2a825a5e7da Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Sun, 6 Dec 2020 00:58:52 +0100 Subject: [PATCH 02/15] review part 2 --- src/LinearMaps.jl | 2 ++ src/composition.jl | 49 ++++++++++++++++++++++++---------------- src/functionmap.jl | 22 +++++++++++------- src/linearcombination.jl | 29 ++++++++++++------------ src/scaledmap.jl | 4 ++-- src/show.jl | 2 +- src/transpose.jl | 1 + test/composition.jl | 2 +- 8 files changed, 66 insertions(+), 45 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index f0206e66..25b6929b 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -235,6 +235,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 diff --git a/src/composition.jl b/src/composition.jl index 844382b8..564818f8 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -1,18 +1,18 @@ -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) + 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)) @@ -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,30 +118,40 @@ 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) - 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 +const CompositeMap1 = CompositeMap{<:Any,<:Tuple{LinearMap}} +const CompositeMap2 = CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}} +_unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap1, x::AbstractVector) = + _unsafe_mul!(y, A.maps[1], x) +# function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap2, x::AbstractVector) +# _compositemul!(y, A, x, similar(y, size(A.maps[1], 1))) +# end +# function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) +# _compositemul!(y, A, x, similar(y, size(A.maps[1], 1)), similar(y, size(A.maps[2], 1))) +# end +_unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) = + _compositemul!(y, A, x) -function _compositemul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, z::AbstractVector) +function _compositemul!(y::AbstractVecOrMat, A::CompositeMap2, x::AbstractVector) + # was there any advantage in having this z an argument, this is not a public method? + z = similar(y, size(A.maps[1], 1)) # no size checking, will be done by individual maps _unsafe_mul!(z, A.maps[1], x) _unsafe_mul!(y, A.maps[2], z) return y end -function _compositemul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{Vararg{LinearMap}}}, x::AbstractVector, source::AbstractVector, dest::AbstractVector) +function _compositemul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) + # was there any advantage in having these arguments, this is not a public method? + source = similar(y, size(A.maps[1], 1)) + dest = similar(y, size(A.maps[2], 1)) # no size checking, will be done by individual maps N = length(A.maps) _unsafe_mul!(source, A.maps[1], x) diff --git a/src/functionmap.jl b/src/functionmap.jl index 8a7f97aa..585302d8 100644 --- a/src/functionmap.jl +++ b/src/functionmap.jl @@ -17,13 +17,16 @@ 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) +Base.parent(A::FunctionMap) = (A.f, A.fc) # is this meanigful/useful? LinearAlgebra.issymmetric(A::FunctionMap) = A._issymmetric LinearAlgebra.ishermitian(A::FunctionMap) = A._ishermitian LinearAlgebra.isposdef(A::FunctionMap) = A._isposdef @@ -31,6 +34,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 +47,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 +73,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 +111,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 +132,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/linearcombination.jl b/src/linearcombination.jl index 13fa9abc..1906f70b 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) @@ -19,9 +19,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 +60,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 +110,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 167c6e7e..d39f388e 100644 --- a/src/scaledmap.jl +++ b/src/scaledmap.jl @@ -15,8 +15,8 @@ ScaledMap(λ::S, lmap::A) where {S<:RealOrComplex,A<:LinearMap} = 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.issymmetric(A::ScaledMap) = isreal(A.λ) && issymmetric(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) diff --git a/src/show.jl b/src/show.jl index 5b3ec08e..298361f1 100644 --- a/src/show.jl +++ b/src/show.jl @@ -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 3e840f73..70c8b6d7 100644 --- a/src/transpose.jl +++ b/src/transpose.jl @@ -82,6 +82,7 @@ _unsafe_mul!(y::AbstractMatrix, A::AdjointMap, x::AbstractMatrix, α::Number, β 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::$Out, Ac::ConjugateMap, x::$In) diff --git a/test/composition.jl b/test/composition.jl index 3a538254..d1780275 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)) From 04b923d228050b7a95172fd51c70b5f8aa5b5e05 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Sun, 6 Dec 2020 22:13:33 +0100 Subject: [PATCH 03/15] review part 2b: new compositemul attempt --- src/composition.jl | 37 ++++++++++++++++++------------------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/src/composition.jl b/src/composition.jl index 564818f8..a19b3e6c 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -127,31 +127,30 @@ LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} = Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps) # multiplication with vectors -const CompositeMap1 = CompositeMap{<:Any,<:Tuple{LinearMap}} -const CompositeMap2 = CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}} -_unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap1, x::AbstractVector) = - _unsafe_mul!(y, A.maps[1], x) -# function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap2, x::AbstractVector) -# _compositemul!(y, A, x, similar(y, size(A.maps[1], 1))) -# end -# function _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) -# _compositemul!(y, A, x, similar(y, size(A.maps[1], 1)), similar(y, size(A.maps[2], 1))) -# end _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) = _compositemul!(y, A, x) - -function _compositemul!(y::AbstractVecOrMat, A::CompositeMap2, x::AbstractVector) +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 _compositemul!(y::AbstractVecOrMat, + A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, + source = similar(y, size(A.maps[1], 1), + dest = nothing) # was there any advantage in having this z an argument, this is not a public method? - z = similar(y, size(A.maps[1], 1)) # no size checking, will be done by individual maps - _unsafe_mul!(z, A.maps[1], x) - _unsafe_mul!(y, A.maps[2], z) + _unsafe_mul!(source, A.maps[1], x) + _unsafe_mul!(y, A.maps[2], source) return y end -function _compositemul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) - # was there any advantage in having these arguments, this is not a public method? - source = similar(y, size(A.maps[1], 1)) - dest = similar(y, size(A.maps[2], 1)) +function _compositemul!(y::AbstractVecOrMat, + A::CompositeMap, + x::AbstractVector, + source = similar(y, size(A.maps[1], 1), + dest = similar(y, size(A.maps[2], 1)) # no size checking, will be done by individual maps N = length(A.maps) _unsafe_mul!(source, A.maps[1], x) From 59312cd814e5eac3839f4974bcec3c851624b705 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Sun, 6 Dec 2020 22:16:53 +0100 Subject: [PATCH 04/15] typo (don't push without testing locally :D ) --- src/composition.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/composition.jl b/src/composition.jl index a19b3e6c..bea2bcca 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -138,7 +138,7 @@ function _compositemul!(y::AbstractVecOrMat, end function _compositemul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, - source = similar(y, size(A.maps[1], 1), + source = similar(y, size(A.maps[1], 1)), dest = nothing) # was there any advantage in having this z an argument, this is not a public method? # no size checking, will be done by individual maps @@ -149,8 +149,8 @@ end function _compositemul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector, - source = similar(y, size(A.maps[1], 1), - dest = similar(y, size(A.maps[2], 1)) + source = similar(y, size(A.maps[1], 1)), + dest = similar(y, size(A.maps[2], 1))) # no size checking, will be done by individual maps N = length(A.maps) _unsafe_mul!(source, A.maps[1], x) From e7bb9f5ef47ace5717016490edbf733f2c90cea9 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Sun, 6 Dec 2020 22:20:47 +0100 Subject: [PATCH 05/15] review part 2c: add comment --- src/composition.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/composition.jl b/src/composition.jl index bea2bcca..dc1f8d5b 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -6,6 +6,7 @@ struct CompositeMap{T, As<:LinearMapTuple} <: LinearMap{T} check_dim_mul(maps[n], maps[n-1]) end 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 From b1208fb5d6164d865e9ec8eea7ceff12d5ad2b6e Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Mon, 7 Dec 2020 00:39:18 +0100 Subject: [PATCH 06/15] review part 3: blockmap --- src/LinearMaps.jl | 1 + src/blockmap.jl | 175 +++++++++++++++++++++++++++------------------- test/blockmap.jl | 4 +- 3 files changed, 105 insertions(+), 75 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 25b6929b..c18be429 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -16,6 +16,7 @@ end abstract type LinearMap{T} end +const MapOrVecOrMat{T} = Union{LinearMap{T}, AbstractVecOrMat{T}} const MapOrMatrix{T} = Union{LinearMap{T}, AbstractMatrix{T}} const RealOrComplex = Union{Real, Complex} diff --git a/src/blockmap.jl b/src/blockmap.jl index d31281b9..e877562a 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -1,19 +1,26 @@ -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...) @@ -27,24 +34,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 +80,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,7 +92,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")) + errmsg = "hcat of only UniformScaling objects cannot determine the linear map size" + nrows == -1 && throw(ArgumentError(errmsg)) + # does this mean we still have type piracy? I vagely remember this being discussed return BlockMap{T}(promote_to_lmaps(fill(nrows, nbc), 1, 1, As...), (nbc,)) end @@ -124,8 +135,9 @@ 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")) - + errmsg = "vcat of only UniformScaling objects cannot determine the linear map size" + ncols == -1 && throw(ArgumentError(errmsg)) + # same question about type piracy return BlockMap{T}(promote_to_lmaps(fill(ncols, nbr), 1, 2, As...), ntuple(i->1, nbr)) end @@ -160,10 +172,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: @@ -173,13 +187,14 @@ function Base.hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling na = size(As[j+k], 1) ni >= 0 && ni != na && throw(DimensionMismatch("mismatch in number of rows")) + # here we check all maps, in hcat and vcat we only find the first map + # which is not UniformScaling to set the size. + # Should we also include this check in hcat and vcat? ni = na 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 +233,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 +252,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 +275,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 +290,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 +318,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 +331,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 +355,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 +380,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 +410,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 +429,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 +441,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 +472,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 @@ -468,19 +495,21 @@ 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/test/blockmap.jl b/test/blockmap.jl index 79b37aca..7e330971 100644 --- a/test/blockmap.jl +++ b/test/blockmap.jl @@ -12,7 +12,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive @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 +193,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 From f000f83579d47533d7a92ee3b2ccd21da91fcdb2 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Tue, 8 Dec 2020 00:15:42 +0100 Subject: [PATCH 07/15] review part 3b: change comments --- src/blockmap.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/blockmap.jl b/src/blockmap.jl index e877562a..88e281b5 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -92,10 +92,9 @@ function Base.hcat(As::Union{LinearMap, UniformScaling, AbstractVecOrMat}...) break end end - errmsg = "hcat of only UniformScaling objects cannot determine the linear map size" - nrows == -1 && throw(ArgumentError(errmsg)) - # does this mean we still have type piracy? I vagely remember this being discussed - 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 ############ @@ -135,10 +134,10 @@ function Base.vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...) break end end - errmsg = "vcat of only UniformScaling objects cannot determine the linear map size" - ncols == -1 && throw(ArgumentError(errmsg)) - # same question about type piracy - 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 ############ From ba51f26432d49262929b4d590bf0c66399e11bf3 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Tue, 8 Dec 2020 00:15:58 +0100 Subject: [PATCH 08/15] review part 4:kronecker --- src/kronecker.jl | 106 ++++++++++++++++++++++++++-------------------- test/kronecker.jl | 2 +- 2 files changed, 62 insertions(+), 46 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 5837bdc0..48514f0f 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,30 @@ 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, Cs::LinearMap...) = + KroneckerMap{promote_type(eltype(A), eltype(B), map(eltype, Cs)...)}(tuple(A, B, Cs...)) + # could this just be a recursive definition: kron(A, B, C) = kron(kron(A, B), C) ? 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,10 +85,10 @@ 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)...) @@ -94,8 +98,8 @@ 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,11 +122,11 @@ 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)) @@ -134,24 +138,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 +172,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 +193,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,7 +258,8 @@ 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)) @@ -256,10 +269,13 @@ 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/test/kronecker.jl b/test/kronecker.jl index f9031a62..75a56d6f 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -9,7 +9,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays 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() From 7d3300b830abd5cd124bdd71d1429d1c877b2b70 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Tue, 8 Dec 2020 23:30:39 +0100 Subject: [PATCH 09/15] review part 4b: recursive kron --- src/kronecker.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 48514f0f..6145729a 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -52,9 +52,8 @@ 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...)) - # could this just be a recursive definition: kron(A, B, C) = kron(kron(A, B), C) ? +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 From 122baed0aef9fab38ffd8b40f631ddc43c240e10 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Tue, 8 Dec 2020 23:31:07 +0100 Subject: [PATCH 10/15] review part 5: conversion spaces and questions --- src/conversion.jl | 54 +++++++++++++++++++++++++++++++++-------------- 1 file changed, 38 insertions(+), 16 deletions(-) diff --git a/src/conversion.jl b/src/conversion.jl index 21c0e932..e99a929e 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,29 @@ 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} +# Are all these methods below really useful? +# The conversion to matrices is not something that belongs to ordinary use of LinearMaps. +# It's rather for debugging purposes, where efficiency is not the primary concern. +# Also, I don't think they are really that much more efficient +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 +98,48 @@ 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) + # if you want less allocations: lmul!(J.λ, convert(Matrix{T}, 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) + # probably also works with sparse: lmul!(J.λ, convert(SparseMatrixCSC, 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.λ) + # if you want less allocations: rmul!(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.λ) + # probably also works with sparse: rmul!(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)...) +# seems like code duplication: +# convert(AbstractMatrix, A) will anyway fall back to Matrix{eltype(A)}(A) + function SparseArrays.sparse(A::BlockMap) return hvcat( A.rows, @@ -140,13 +155,17 @@ 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)...) +# seems like code duplication: +# convert(AbstractMatrix, A) will anyway fall back to Matrix{eltype(A)}(A) 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))) @@ -159,6 +178,9 @@ function Base.convert(::Type{AbstractMatrix}, L::KroneckerSumMap) IB = Diagonal(ones(Bool, size(B, 1))) return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B)) end +# seems like code duplication: +# convert(AbstractMatrix, L) will anyway fall back to Matrix{eltype(L)}(L) + function SparseArrays.sparse(L::KroneckerSumMap) A, B = L.maps IA = sparse(Diagonal(ones(Bool, size(A, 1)))) From cb63bbcbaa8d5f85e2a6e84d091b507194e1dffb Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Tue, 8 Dec 2020 23:33:31 +0100 Subject: [PATCH 11/15] review part 6: show, linebreak for 92-character rule --- src/show.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/show.jl b/src/show.jl index 298361f1..936fa958 100644 --- a/src/show.jl +++ b/src/show.jl @@ -26,8 +26,8 @@ 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) From 025cd53f723906a79855ecdd58850a7b83971c28 Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Wed, 9 Dec 2020 14:00:11 +0100 Subject: [PATCH 12/15] review part 6b: remove comments --- src/conversion.jl | 20 ++------------------ 1 file changed, 2 insertions(+), 18 deletions(-) diff --git a/src/conversion.jl b/src/conversion.jl index e99a929e..ba80130f 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -82,10 +82,6 @@ function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{Matrix end # CompositeMap -# Are all these methods below really useful? -# The conversion to matrices is not something that belongs to ordinary use of LinearMaps. -# It's rather for debugging purposes, where efficiency is not the primary concern. -# Also, I don't think they are really that much more efficient function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, LinearMap}}) where {T} B, A = AB.maps require_one_based_indexing(B) @@ -115,31 +111,24 @@ end function Base.Matrix{T}(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}}) where {T} A, J = λA.maps return convert(Matrix{T}, J.λ*A.lmap) - # if you want less allocations: lmul!(J.λ, convert(Matrix{T}, A.lmap)) end function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}}) A, J = λA.maps return convert(SparseMatrixCSC, J.λ*A.lmap) - # probably also works with sparse: lmul!(J.λ, convert(SparseMatrixCSC, A.lmap)) - end function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}}) where {T} J, A = Aλ.maps return convert(Matrix{T}, A.lmap*J.λ) - # if you want less allocations: rmul!(convert(Matrix{T}, A.lmap), J.λ) end function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}}) J, A = Aλ.maps return convert(SparseMatrixCSC, A.lmap*J.λ) - # probably also works with sparse: rmul!(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)...) -# seems like code duplication: -# convert(AbstractMatrix, A) will anyway fall back to Matrix{eltype(A)}(A) - +Base.convert(::Type{AbstractMatrix}, A::BlockMap) = + hvcat(A.rows, convert.(AbstractMatrix, A.maps)...) function SparseArrays.sparse(A::BlockMap) return hvcat( A.rows, @@ -157,8 +146,6 @@ end Base.Matrix{T}(A::KroneckerMap) where {T} = kron(convert.(Matrix{T}, A.maps)...) Base.convert(::Type{AbstractMatrix}, A::KroneckerMap) = kron(convert.(AbstractMatrix, A.maps)...) -# seems like code duplication: -# convert(AbstractMatrix, A) will anyway fall back to Matrix{eltype(A)}(A) function SparseArrays.sparse(A::KroneckerMap) return kron( convert(SparseMatrixCSC, first(A.maps)), @@ -178,9 +165,6 @@ function Base.convert(::Type{AbstractMatrix}, L::KroneckerSumMap) IB = Diagonal(ones(Bool, size(B, 1))) return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B)) end -# seems like code duplication: -# convert(AbstractMatrix, L) will anyway fall back to Matrix{eltype(L)}(L) - function SparseArrays.sparse(L::KroneckerSumMap) A, B = L.maps IA = sparse(Diagonal(ones(Bool, size(A, 1)))) From 94bab043d5d69ac5e59077448552a0a57ee1cbeb Mon Sep 17 00:00:00 2001 From: Jutho Haegeman Date: Wed, 9 Dec 2020 15:42:47 +0100 Subject: [PATCH 13/15] clean up final comments and suggestions --- src/blockmap.jl | 3 --- src/composition.jl | 1 - src/functionmap.jl | 2 +- src/scaledmap.jl | 2 +- 4 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/blockmap.jl b/src/blockmap.jl index 88e281b5..1f88c50b 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -186,9 +186,6 @@ function Base.hvcat(rows::Tuple{Vararg{Int}}, na = size(As[j+k], 1) ni >= 0 && ni != na && throw(DimensionMismatch("mismatch in number of rows")) - # here we check all maps, in hcat and vcat we only find the first map - # which is not UniformScaling to set the size. - # Should we also include this check in hcat and vcat? ni = na end end diff --git a/src/composition.jl b/src/composition.jl index dc1f8d5b..c359aa5b 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -141,7 +141,6 @@ function _compositemul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, source = similar(y, size(A.maps[1], 1)), dest = nothing) - # was there any advantage in having this z an argument, this is not a public method? # no size checking, will be done by individual maps _unsafe_mul!(source, A.maps[1], x) _unsafe_mul!(y, A.maps[2], source) diff --git a/src/functionmap.jl b/src/functionmap.jl index 585302d8..b665e8b7 100644 --- a/src/functionmap.jl +++ b/src/functionmap.jl @@ -26,7 +26,7 @@ FunctionMap{T}(f, fc, M::Int; kwargs...) where {T} = # properties Base.size(A::FunctionMap) = (A.M, A.N) -Base.parent(A::FunctionMap) = (A.f, A.fc) # is this meanigful/useful? +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 diff --git a/src/scaledmap.jl b/src/scaledmap.jl index d39f388e..298cbbbe 100644 --- a/src/scaledmap.jl +++ b/src/scaledmap.jl @@ -15,7 +15,7 @@ ScaledMap(λ::S, lmap::A) where {S<:RealOrComplex,A<:LinearMap} = 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) = isreal(A.λ) && issymmetric(A.lmap) +LinearAlgebra.issymmetric(A::ScaledMap) = issymmetric(A.lmap) LinearAlgebra.ishermitian(A::ScaledMap) = isreal(A.λ) && ishermitian(A.lmap) LinearAlgebra.isposdef(A::ScaledMap) = isposdef(A.λ) && isposdef(A.lmap) From 10b2531d0a76098c5a440995433a2d070e2d835d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 9 Dec 2020 21:08:02 +0100 Subject: [PATCH 14/15] remove parent, remove a few comments --- docs/src/types.md | 6 ------ src/blockmap.jl | 4 ---- src/composition.jl | 3 --- src/functionmap.jl | 1 - src/kronecker.jl | 5 +---- src/left.jl | 2 -- src/linearcombination.jl | 1 - src/scaledmap.jl | 1 - src/show.jl | 2 +- src/transpose.jl | 2 -- src/uniformscalingmap.jl | 2 +- src/wrappedmap.jl | 2 +- test/blockmap.jl | 2 -- test/composition.jl | 1 - test/functionmap.jl | 1 - test/kronecker.jl | 2 -- test/linearcombination.jl | 1 - test/linearmaps.jl | 1 - test/scaledmap.jl | 1 - test/transpose.jl | 2 -- test/uniformscalingmap.jl | 2 +- test/wrappedmap.jl | 2 +- 22 files changed, 6 insertions(+), 40 deletions(-) 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/blockmap.jl b/src/blockmap.jl index 1f88c50b..1e5b5a88 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -24,8 +24,6 @@ BlockMap{T}(maps::As, rows::Rs) where {T, As<:LinearMapTuple, Rs} = MulStyle(A::BlockMap) = MulStyle(A.maps...) -Base.parent(A::BlockMap) = A.maps - """ rowcolranges(maps, rows) @@ -485,8 +483,6 @@ 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) diff --git a/src/composition.jl b/src/composition.jl index c359aa5b..9f27418c 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -18,7 +18,6 @@ CompositeMap{T}(maps::As) where {T, As<:LinearMapTuple} = CompositeMap{T, As}(ma # 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), @@ -141,7 +140,6 @@ function _compositemul!(y::AbstractVecOrMat, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector, source = similar(y, size(A.maps[1], 1)), dest = nothing) - # no size checking, will be done by individual maps _unsafe_mul!(source, A.maps[1], x) _unsafe_mul!(y, A.maps[2], source) return y @@ -151,7 +149,6 @@ function _compositemul!(y::AbstractVecOrMat, x::AbstractVector, source = similar(y, size(A.maps[1], 1)), dest = similar(y, size(A.maps[2], 1))) - # no size checking, will be done by individual maps N = length(A.maps) _unsafe_mul!(source, A.maps[1], x) for n in 2:N-1 diff --git a/src/functionmap.jl b/src/functionmap.jl index b665e8b7..6c80b3a9 100644 --- a/src/functionmap.jl +++ b/src/functionmap.jl @@ -26,7 +26,6 @@ FunctionMap{T}(f, fc, M::Int; kwargs...) where {T} = # 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 diff --git a/src/kronecker.jl b/src/kronecker.jl index 6145729a..289a689b 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -91,8 +91,6 @@ Base.:(^)(A::MapOrMatrix, ::KronPower{p}) where {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) @@ -128,7 +126,7 @@ end 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 @@ -262,7 +260,6 @@ Base.:(^)(A::MapOrMatrix, ::KronSumPower{p}) where {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) diff --git a/src/left.jl b/src/left.jl index 9ad42d57..4f6bddff 100644 --- a/src/left.jl +++ b/src/left.jl @@ -50,7 +50,6 @@ end function mul!(x::AbstractMatrix, y::AdjointAbsVecOrMat, A::LinearMap, α::Number, β::Number) check_dim_mul(x, y, A) _unsafe_mul!(x', conj(α)*A', y', true, conj(β)) - # why not?: _unsafe_mul!(x', A', y', conj(α), conj(β)) return x end @@ -63,6 +62,5 @@ end function mul!(x::AbstractMatrix, y::TransposeAbsVecOrMat, A::LinearMap, α::Number, β::Number) check_dim_mul(x, y, A) _unsafe_mul!(transpose(x), α*transpose(A), transpose(y), true, β) - # why not?: _unsafe_mul!(transpose(x), transpose(A), transpose(y), α, β) return x end diff --git a/src/linearcombination.jl b/src/linearcombination.jl index 1906f70b..ffaccd11 100644 --- a/src/linearcombination.jl +++ b/src/linearcombination.jl @@ -18,7 +18,6 @@ MulStyle(A::LinearCombination) = MulStyle(A.maps...) # basic methods Base.size(A::LinearCombination) = size(A.maps[1]) -Base.parent(A::LinearCombination) = A.maps # following conditions are sufficient but not necessary LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps) LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps) diff --git a/src/scaledmap.jl b/src/scaledmap.jl index 298cbbbe..23cb7534 100644 --- a/src/scaledmap.jl +++ b/src/scaledmap.jl @@ -13,7 +13,6 @@ 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) = isreal(A.λ) && ishermitian(A.lmap) diff --git a/src/show.jl b/src/show.jl index 936fa958..3b5ca5c2 100644 --- a/src/show.jl +++ b/src/show.jl @@ -21,7 +21,7 @@ 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) diff --git a/src/transpose.jl b/src/transpose.jl index 70c8b6d7..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) diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 462920ab..350812bf 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -17,7 +17,7 @@ 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) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index d149af46..4631c246 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -38,7 +38,7 @@ Base.:(==)(A::MatrixMap, B::MatrixMap) = # 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 diff --git a/test/blockmap.jl b/test/blockmap.jl index 7e330971..10fc1cd7 100644 --- a/test/blockmap.jl +++ b/test/blockmap.jl @@ -7,7 +7,6 @@ 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} @@ -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 d1780275..cb2d42a1 100644 --- a/test/composition.jl +++ b/test/composition.jl @@ -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 75a56d6f..6713610b 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -8,7 +8,6 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays LA = LinearMap(A) LB = LinearMap(B) LK = @inferred kron(LA, LB) - @test parent(LK) == (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) @@ -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..b86c9873 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) 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 From 85e350e0ff98f34f7837d41d8da3f2d34a0833b8 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 9 Dec 2020 21:31:31 +0100 Subject: [PATCH 15/15] improve coverage, minor kronecker fix --- src/LinearMaps.jl | 9 --------- src/kronecker.jl | 2 +- test/kronecker.jl | 2 +- test/transpose.jl | 2 ++ 4 files changed, 4 insertions(+), 11 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index c18be429..cc9b146d 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -50,15 +50,6 @@ 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) diff --git a/src/kronecker.jl b/src/kronecker.jl index 289a689b..f6beefe5 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -84,7 +84,7 @@ Construct a lazy representation of the `k`-th Kronecker power """ ⊗(k::Integer) = KronPower(k) -⊗(A, B, Cs...) = kron(convert_to_lmaps(A, B, Cs...)) +⊗(A, B, Cs...) = kron(convert_to_lmaps(A, B, Cs...)...) Base.:(^)(A::MapOrMatrix, ::KronPower{p}) where {p} = kron(ntuple(n -> convert_to_lmaps_(A), Val(p))...) diff --git a/test/kronecker.jl b/test/kronecker.jl index 6713610b..7acd6e78 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -23,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) diff --git a/test/transpose.jl b/test/transpose.jl index b86c9873..ef240505 100644 --- a/test/transpose.jl +++ b/test/transpose.jl @@ -27,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