From 032832e7d887bf9ad94e19ade4aa6c2044a31039 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 24 Mar 2021 14:42:22 +0100 Subject: [PATCH] Allow wrapping of vectors, extend Kronecker product (#136) --- Project.toml | 2 +- docs/src/history.md | 8 ++++++++ src/LinearMaps.jl | 18 +++++++++--------- src/blockmap.jl | 10 ++++------ src/conversion.jl | 32 +++++++++++++++++--------------- src/kronecker.jl | 15 +++++++++------ src/show.jl | 2 +- src/wrappedmap.jl | 40 +++++++++++++++++++++++----------------- test/blockmap.jl | 4 +++- test/conversion.jl | 3 +++ test/kronecker.jl | 4 ++++ 11 files changed, 82 insertions(+), 56 deletions(-) diff --git a/Project.toml b/Project.toml index 3352797e..72528e48 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "LinearMaps" uuid = "7a12625a-238d-50fd-b39a-03d52299707e" -version = "3.2.4" +version = "3.3.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/docs/src/history.md b/docs/src/history.md index 52905427..94296a9f 100644 --- a/docs/src/history.md +++ b/docs/src/history.md @@ -1,5 +1,13 @@ # Version history +## What's new in v3.3 + +* `AbstractVector`s can now be wrapped by a `LinearMap` just like `AbstractMatrix`` + typed objects. Upon wrapping, there are not implicitly reshaped to matrices. This + feature might be helpful, for instance, in the lazy representation of rank-1 + operators `kron(LinearMap(u), v') == ⊗(u, v') == u ⊗ v'` for vectors `u` and `v`. + The action on vectors,`(u⊗v')*x`, is implemented optimally via `u*(v'x)`. + ## What's new in v3.2 * In-place left-multiplication `mul!(Y, X, A::LinearMap)` is now allowed for diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 6d43d10e..d3355e5d 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -34,9 +34,9 @@ MulStyle(::FiveArg, ::ThreeArg) = ThreeArg() MulStyle(::ThreeArg, ::ThreeArg) = ThreeArg() MulStyle(::LinearMap) = ThreeArg() # default @static if VERSION ≥ v"1.3.0-alpha.115" - MulStyle(::AbstractMatrix) = FiveArg() + MulStyle(::AbstractVecOrMat) = FiveArg() else - MulStyle(::AbstractMatrix) = ThreeArg() + MulStyle(::AbstractVecOrMat) = ThreeArg() end MulStyle(A::LinearMap, As::LinearMap...) = MulStyle(MulStyle(A), MulStyle(As...)) @@ -68,8 +68,8 @@ function check_dim_mul(C, A, B) return nothing end -# conversion of AbstractMatrix to LinearMap -convert_to_lmaps_(A::AbstractMatrix) = LinearMap(A) +# conversion of AbstractVecOrMat to LinearMap +convert_to_lmaps_(A::AbstractVecOrMat) = LinearMap(A) convert_to_lmaps_(A::LinearMap) = A convert_to_lmaps() = () convert_to_lmaps(A) = (convert_to_lmaps_(A),) @@ -231,8 +231,8 @@ function _generic_mapmat_mul!(Y, A, X, α=true, β=false) return Y end -_unsafe_mul!(y, A::MapOrMatrix, x) = mul!(y, A, x) -_unsafe_mul!(y, A::AbstractMatrix, x, α, β) = mul!(y, A, x, α, β) +_unsafe_mul!(y, A::MapOrVecOrMat, x) = mul!(y, A, x) +_unsafe_mul!(y, A::AbstractVecOrMat, x, α, β) = mul!(y, A, x, α, β) function _unsafe_mul!(y::AbstractVecOrMat, A::LinearMap, x::AbstractVector, α, β) return _generic_mapvec_mul!(y, A, x, α, β) end @@ -261,11 +261,11 @@ include("show.jl") # show methods for LinearMap objects """ LinearMap(A::LinearMap; kwargs...)::WrappedMap - LinearMap(A::AbstractMatrix; kwargs...)::WrappedMap + LinearMap(A::AbstractVecOrMat; kwargs...)::WrappedMap LinearMap(J::UniformScaling, M::Int)::UniformScalingMap LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)::FunctionMap -Construct a linear map object, either from an existing `LinearMap` or `AbstractMatrix` `A`, +Construct a linear map object, either from an existing `LinearMap` or `AbstractVecOrMat` `A`, with the purpose of redefining its properties via the keyword arguments `kwargs`; a `UniformScaling` object `J` with specified (square) dimension `M`; from a `Number` object to lazily represent filled matrices; or @@ -293,7 +293,7 @@ For the function-based constructor, there is one more keyword argument: 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(A::MapOrVecOrMat; kwargs...) = WrappedMap(A; kwargs...) LinearMap(J::UniformScaling, M::Int) = UniformScalingMap(J.λ, M) LinearMap(f, M::Int; kwargs...) = LinearMap{Float64}(f, M; kwargs...) LinearMap(f, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, M, N; kwargs...) diff --git a/src/blockmap.jl b/src/blockmap.jl index 1e5b5a88..cc67069b 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -38,15 +38,15 @@ function rowcolranges(maps, rows) rowstart = 1 for (i, row) in enumerate(rows) mapind += 1 - rowend = rowstart + size(maps[mapind], 1) - 1 + rowend = rowstart + Int(size(maps[mapind], 1))::Int - 1 rowranges = Base.setindex(rowranges, rowstart:rowend, i) colstart = 1 - colend = size(maps[mapind], 2) + colend = Int(size(maps[mapind], 2))::Int colranges = Base.setindex(colranges, colstart:colend, mapind) for colind in 2:row mapind += 1 colstart = colend + 1 - colend += size(maps[mapind], 2) + colend += Int(size(maps[mapind], 2))::Int colranges = Base.setindex(colranges, colstart:colend, mapind) end rowstart = rowend + 1 @@ -226,9 +226,7 @@ function check_dim(A, dim, n) return nothing 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::AbstractVecOrMat) = (check_dim(A, dim, n); LinearMap(A)) 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) = () diff --git a/src/conversion.jl b/src/conversion.jl index 8588a5c3..309150f4 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -47,10 +47,10 @@ SparseArrays.SparseMatrixCSC(A::LinearMap) = sparse(A) # special cases # ScaledMap -Base.Matrix{T}(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) where {T} = +Base.Matrix{T}(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) where {T} = convert(Matrix{T}, A.λ * A.lmap.lmap) -SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) = - convert(SparseMatrixCSC, A.λ * A.lmap.lmap) +SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) = + A.λ * sparse(A.lmap.lmap) # UniformScalingMap Base.Matrix{T}(J::UniformScalingMap) where {T} = Matrix{T}(J.λ*I, size(J)) @@ -59,6 +59,8 @@ Base.convert(::Type{AbstractMatrix}, J::UniformScalingMap) = Diagonal(fill(J.λ, # WrappedMap Base.Matrix{T}(A::WrappedMap) where {T} = Matrix{T}(A.lmap) Base.convert(::Type{T}, A::WrappedMap) where {T<:Matrix} = convert(T, A.lmap) +Base.Matrix{T}(A::VectorMap) where {T} = copyto!(Matrix{eltype(T)}(undef, size(A)), A.lmap) +Base.convert(::Type{T}, A::VectorMap) where {T<:Matrix} = T(A) Base.convert(::Type{AbstractMatrix}, A::WrappedMap) = convert(AbstractMatrix, A.lmap) SparseArrays.sparse(A::WrappedMap) = sparse(A.lmap) Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap) @@ -70,19 +72,19 @@ for (T, t) in ((AdjointMap, adjoint), (TransposeMap, transpose)) end # LinearCombination -function Base.Matrix{T}(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}}) where {T} +function Base.Matrix{T}(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}}) where {T} maps = ΣA.maps mats = map(A->getfield(A, :lmap), maps) return Matrix{T}(sum(mats)) end -function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}}) +function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}}) maps = ΣA.maps mats = map(A->getfield(A, :lmap), maps) return convert(SparseMatrixCSC, sum(mats)) end # CompositeMap -function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, LinearMap}}) where {T} +function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, LinearMap}}) where {T} B, A = AB.maps require_one_based_indexing(B) Y = Matrix{eltype(AB)}(undef, size(AB)) @@ -91,36 +93,36 @@ function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, LinearMap}}) end return Y end -for ((TA, fieldA), (TB, fieldB)) in (((MatrixMap, :lmap), (MatrixMap, :lmap)), - ((MatrixMap, :lmap), (UniformScalingMap, :λ)), - ((UniformScalingMap, :λ), (MatrixMap, :lmap))) +for ((TA, fieldA), (TB, fieldB)) in (((VecOrMatMap, :lmap), (VecOrMatMap, :lmap)), + ((VecOrMatMap, :lmap), (UniformScalingMap, :λ)), + ((UniformScalingMap, :λ), (VecOrMatMap, :lmap))) @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{VecOrMatMap, VecOrMatMap}}) 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{VecOrMatMap, VecOrMatMap}}) 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{VecOrMatMap, UniformScalingMap}}) where {T} A, J = λA.maps return convert(Matrix{T}, J.λ*A.lmap) end -function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}}) +function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}}) A, J = λA.maps return convert(SparseMatrixCSC, J.λ*A.lmap) end -function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}}) where {T} +function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}}) where {T} J, A = Aλ.maps return convert(Matrix{T}, A.lmap*J.λ) end -function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}}) +function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}}) J, A = Aλ.maps return convert(SparseMatrixCSC, A.lmap*J.λ) end diff --git a/src/kronecker.jl b/src/kronecker.jl index a627b9e4..4443b90f 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -64,18 +64,18 @@ Base.kron(A::KroneckerMap, B::ScaledMap) = kron(A, B.lmap) * B.λ # generic definitions 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)) +Base.kron(A::AbstractVecOrMat, B::LinearMap) = kron(LinearMap(A), B) +Base.kron(A::LinearMap, B::AbstractVecOrMat) = 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))::AbstractVecOrMat), Val(k-1)) # yields (:A1, :A2, :A3, ..., :A(k-1)) L = :($(Symbol(:A, k))::LinearMap) # yields :Ak::LinearMap 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...) = + @eval Base.kron($(Is...), $L, As::MapOrVecOrMat...) = kron($(mapargs...), $(Symbol(:A, k)), convert_to_lmaps(As...)...) end @@ -139,7 +139,7 @@ end !isone(A.λ) && rmul!(y, A.λ) return y end -@inline function _kronmul!(y, B, x, A::MatrixMap, _) +@inline function _kronmul!(y, B, x, A::VecOrMatMap, _) ma, na = size(A) mb, nb = size(B) X = reshape(x, (nb, na)) @@ -153,10 +153,13 @@ end elseif nb*ma <= mb*na _unsafe_mul!(Y, B, X * At) else - _unsafe_mul!(Y, Matrix(B*X), At) + _unsafe_mul!(Y, Matrix(B * X), At) end return y end +const VectorMap{T} = WrappedMap{T,<:AbstractVector} +const AdjOrTransVectorMap{T} = WrappedMap{T,<:LinearAlgebra.AdjOrTransAbsVec} +@inline _kronmul!(y, B::AdjOrTransVectorMap, x, a::VectorMap, _) = mul!(y, a.lmap, B.lmap * x) ################# # multiplication with vectors diff --git a/src/show.jl b/src/show.jl index 514ac39d..e0c25856 100644 --- a/src/show.jl +++ b/src/show.jl @@ -8,7 +8,7 @@ end Base.show(io::IO, A::LinearMap) = print(io, map_show(io, A, 0)) map_show(io::IO, A::LinearMap, i) = ' '^i * map_summary(A) * _show(io, A, i) -map_show(io::IO, A::AbstractMatrix, i) = ' '^i * summary(A) +map_show(io::IO, A::AbstractVecOrMat, i) = ' '^i * summary(A) _show(io::IO, ::LinearMap, _) = "" function _show(io::IO, A::FunctionMap{T,F,Nothing}, _) where {T,F} "($(A.f); ismutating=$(A._ismutating), issymmetric=$(A._issymmetric), ishermitian=$(A._ishermitian), isposdef=$(A._isposdef))" diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 1c8fd1ec..3c6ee41c 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -1,43 +1,49 @@ -struct WrappedMap{T, A<:MapOrMatrix} <: LinearMap{T} +struct WrappedMap{T, A<:MapOrVecOrMat} <: LinearMap{T} lmap::A _issymmetric::Bool _ishermitian::Bool _isposdef::Bool end -function WrappedMap(lmap::MapOrMatrix{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} WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef) end +WrappedMap(lmap::MapOrMatrix{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kwargs...) +function WrappedMap{T}(lmap::AbstractVector; + issymmetric::Bool = false, + ishermitian::Bool = false, + isposdef::Bool = false) where {T} + WrappedMap{T, typeof(lmap)}(lmap, + length(lmap) == 1 && isreal(T), + length(lmap) == 1 && isreal(T), + length(lmap) == 1 && isposdef(first(lmap))) +end +WrappedMap(lmap::AbstractVector{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kwargs...) -const MatrixMap{T} = WrappedMap{T,<:AbstractMatrix} +const VecOrMatMap{T} = WrappedMap{T,<:AbstractVecOrMat} -MulStyle(A::WrappedMap) = MulStyle(A.lmap) +MulStyle(A::VecOrMatMap) = MulStyle(A.lmap) -LinearAlgebra.transpose(A::MatrixMap{T}) where {T} = +LinearAlgebra.transpose(A::VecOrMatMap{T}) where {T} = WrappedMap{T}(transpose(A.lmap); issymmetric = A._issymmetric, ishermitian = A._ishermitian, isposdef = A._isposdef) -LinearAlgebra.adjoint(A::MatrixMap{T}) where {T} = +LinearAlgebra.adjoint(A::VecOrMatMap{T}) where {T} = WrappedMap{T}(adjoint(A.lmap); issymmetric = A._issymmetric, ishermitian = A._ishermitian, isposdef = A._isposdef) -Base.:(==)(A::MatrixMap, B::MatrixMap) = +Base.:(==)(A::VecOrMatMap, B::VecOrMatMap) = (eltype(A) == eltype(B) && A.lmap == B.lmap && A._issymmetric == B._issymmetric && A._ishermitian == B._ishermitian && A._isposdef == B._isposdef) # properties Base.size(A::WrappedMap) = size(A.lmap) +Base.size(A::WrappedMap{<:Any,<:AbstractVector}) = (Int(length(A.lmap))::Int, 1) LinearAlgebra.issymmetric(A::WrappedMap) = A._issymmetric LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef @@ -61,9 +67,9 @@ for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractM end end -mul!(Y::AbstractMatrix, X::AbstractMatrix, A::MatrixMap) = mul!(Y, X, A.lmap) +mul!(Y::AbstractMatrix, X::AbstractMatrix, A::VecOrMatMap) = mul!(Y, X, A.lmap) # the following method is needed for disambiguation with left-multiplication -mul!(Y::AbstractMatrix, X::TransposeAbsVecOrMat, A::MatrixMap) = mul!(Y, X, A.lmap) +mul!(Y::AbstractMatrix, X::TransposeAbsVecOrMat, A::VecOrMatMap) = mul!(Y, X, A.lmap) if VERSION ≥ v"1.3.0-alpha.115" for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix)) @@ -87,14 +93,14 @@ if VERSION ≥ v"1.3.0-alpha.115" end end - mul!(X::AbstractMatrix, Y::AbstractMatrix, A::MatrixMap, α::Number, β::Number) = + mul!(X::AbstractMatrix, Y::AbstractMatrix, A::VecOrMatMap, α::Number, β::Number) = mul!(X, Y, A.lmap, α, β) # the following method is needed for disambiguation with left-multiplication - function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::AbstractMatrix{<:RealOrComplex}, A::MatrixMap{<:RealOrComplex}, + function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::AbstractMatrix{<:RealOrComplex}, A::VecOrMatMap{<:RealOrComplex}, α::RealOrComplex, β::RealOrComplex) return mul!(Y, X, A.lmap, α, β) end - function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::TransposeAbsVecOrMat{<:RealOrComplex}, A::MatrixMap{<:RealOrComplex}, + function mul!(Y::AbstractMatrix{<:RealOrComplex}, X::TransposeAbsVecOrMat{<:RealOrComplex}, A::VecOrMatMap{<:RealOrComplex}, α::RealOrComplex, β::RealOrComplex) return mul!(Y, X, A.lmap, α, β) end diff --git a/test/blockmap.jl b/test/blockmap.jl index 10fc1cd7..fa647c44 100644 --- a/test/blockmap.jl +++ b/test/blockmap.jl @@ -26,8 +26,10 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools, Interactive @test (@which [I I I A11 A11 A11]).module != LinearMaps @test (@which hcat(I, I, I)).module != LinearMaps @test (@which hcat(I, I, I, LinearMap(A11), A11, A11)).module == LinearMaps + maps = @inferred LinearMaps.promote_to_lmaps(ntuple(i->10, 7), 1, 1, I, I, I, LinearMap(A11), A11, A11, v) + @inferred LinearMaps.rowcolranges(maps, (7,)) L = @inferred hcat(I, I, I, LinearMap(A11), A11, A11, v) - @test L == [I I I LinearMap(A11) LinearMap(A11) LinearMap(A11) LinearMap(reshape(v, :, 1))] + @test L == [I I I LinearMap(A11) LinearMap(A11) LinearMap(A11) LinearMap(v)] x = rand(elty, 61) @test L isa LinearMaps.BlockMap{elty} @test L * x ≈ A * x diff --git a/test/conversion.jl b/test/conversion.jl index 038cc378..d7a8c8a0 100644 --- a/test/conversion.jl +++ b/test/conversion.jl @@ -10,6 +10,7 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, Quaternions β = rand() M = @inferred LinearMap(A) N = @inferred LinearMap(M) + U = @inferred LinearMap(v) @test Matrix(M) == A @test Array(M) == A @@ -25,6 +26,8 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, Quaternions @test convert(AbstractMatrix, M') isa Adjoint @test convert(Matrix, M*3I) == A*3 @test convert(Matrix, M+M) == A + A + @test all(Matrix(U) .== v) + @test convert(Matrix{ComplexF32}, U) isa Matrix{ComplexF32} # UniformScalingMap J = LinearMap(α*I, 10) diff --git a/test/kronecker.jl b/test/kronecker.jl index 143fac4d..29eaa0b9 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -22,6 +22,10 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays @test @inferred size(LK, i) == size(K, i) end @test LK isa LinearMaps.KroneckerMap{ComplexF64} + L = ones(3) ⊗ ones(ComplexF64, 4)' + v = rand(4) + @test Matrix(L) == ones(3,4) + @test L*v ≈ fill(sum(v), 3) for transform in (identity, transpose, adjoint) @test Matrix(transform(LK)) ≈ transform(Matrix(LK)) ≈ transform(K)