From 699c40789fd509a01cd86f634013c400f06a929e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Sun, 1 Oct 2023 10:18:49 +0200 Subject: [PATCH] Avoid allocation in InnerProduct checks --- src/spaces/vectorspaces.jl | 4 ++++ src/tensors/factorizations.jl | 41 ++++++++++++---------------------- src/tensors/linalg.jl | 25 ++++++++------------- src/tensors/vectorinterface.jl | 3 +-- 4 files changed, 28 insertions(+), 45 deletions(-) diff --git a/src/spaces/vectorspaces.jl b/src/spaces/vectorspaces.jl index f7a69973..bc57dabc 100644 --- a/src/spaces/vectorspaces.jl +++ b/src/spaces/vectorspaces.jl @@ -198,6 +198,10 @@ Return the type of inner product for vector spaces, which can be either InnerProductStyle(V::VectorSpace) = InnerProductStyle(typeof(V)) InnerProductStyle(::Type{<:VectorSpace}) = NoInnerProduct() +@noinline function throw_invalid_innerproduct(fname) + throw(ArgumentError("$fname requires Euclidean inner product")) +end + dual(V::VectorSpace) = dual(InnerProductStyle(V), V) dual(::EuclideanProduct, V::VectorSpace) = conj(V) diff --git a/src/tensors/factorizations.jl b/src/tensors/factorizations.jl index 15084cab..6e51a8df 100644 --- a/src/tensors/factorizations.jl +++ b/src/tensors/factorizations.jl @@ -262,28 +262,22 @@ LinearAlgebra.isposdef(t::AbstractTensorMap) = isposdef!(copy(t)) # only correct if Euclidean inner product #------------------------------------------------------------------------------------------ function leftorth!(t::AdjointTensorMap; alg::OFA=QRpos()) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("leftorth! only defined for Euclidean inner product spaces")) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:leftorth!) return map(adjoint, reverse(rightorth!(adjoint(t); alg=alg'))) end function rightorth!(t::AdjointTensorMap; alg::OFA=LQpos()) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("rightorth! only defined for Euclidean inner product spaces")) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:rightorth!) return map(adjoint, reverse(leftorth!(adjoint(t); alg=alg'))) end -function leftnull!(t::AdjointTensorMap; alg::OFA=QR(), - kwargs...) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("leftnull! only defined for Euclidean inner product spaces")) +function leftnull!(t::AdjointTensorMap; alg::OFA=QR(), kwargs...) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:leftnull!) return adjoint(rightnull!(adjoint(t); alg=alg', kwargs...)) end -function rightnull!(t::AdjointTensorMap; alg::OFA=LQ(), - kwargs...) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("rightnull! only defined for Euclidean inner product spaces")) +function rightnull!(t::AdjointTensorMap; alg::OFA=LQ(), kwargs...) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:rightnull!) return adjoint(leftnull!(adjoint(t); alg=alg', kwargs...)) end @@ -291,9 +285,7 @@ function tsvd!(t::AdjointTensorMap; trunc::TruncationScheme=NoTruncation(), p::Real=2, alg::Union{SVD,SDD}=SDD()) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("tsvd! only defined for Euclidean inner product spaces")) - + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:tsvd!) u, s, vt, err = tsvd!(adjoint(t); trunc=trunc, p=p, alg=alg) return adjoint(vt), adjoint(s), adjoint(u), err end @@ -303,8 +295,7 @@ function leftorth!(t::TensorMap; atol::Real=zero(float(real(scalartype(t)))), rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(scalartype(t)))) : eps(real(float(one(scalartype(t))))) * iszero(atol)) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("leftorth! only defined for Euclidean inner product spaces")) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:leftorth!) if !iszero(rtol) atol = max(atol, rtol * norm(t)) end @@ -340,8 +331,7 @@ function leftnull!(t::TensorMap; atol::Real=zero(float(real(scalartype(t)))), rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(scalartype(t)))) : eps(real(float(one(scalartype(t))))) * iszero(atol)) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("leftnull! only defined for Euclidean inner product spaces")) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:leftnull!) if !iszero(rtol) atol = max(atol, rtol * norm(t)) end @@ -365,8 +355,7 @@ function rightorth!(t::TensorMap; atol::Real=zero(float(real(scalartype(t)))), rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(scalartype(t)))) : eps(real(float(one(scalartype(t))))) * iszero(atol)) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("rightorth! only defined for Euclidean inner product spaces")) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:rightorth!) if !iszero(rtol) atol = max(atol, rtol * norm(t)) end @@ -402,8 +391,7 @@ function rightnull!(t::TensorMap; atol::Real=zero(float(real(scalartype(t)))), rtol::Real=(alg ∉ (SVD(), SDD())) ? zero(float(real(scalartype(t)))) : eps(real(float(one(scalartype(t))))) * iszero(atol)) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("rightnull! only defined for Euclidean inner product spaces")) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:rightnull!) if !iszero(rtol) atol = max(atol, rtol * norm(t)) end @@ -426,8 +414,7 @@ function tsvd!(t::TensorMap; trunc::TruncationScheme=NoTruncation(), p::Real=2, alg::Union{SVD,SDD}=SDD()) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("tsvd! only defined for Euclidean inner product spaces")) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:tsvd!) S = spacetype(t) I = sectortype(t) A = storagetype(t) @@ -500,8 +487,7 @@ end LinearAlgebra.eigen!(t::TensorMap) = ishermitian(t) ? eigh!(t) : eig!(t) function eigh!(t::TensorMap; kwargs...) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("eigh! only defined for Euclidean inner product spaces")) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:eigh!) domain(t) == codomain(t) || throw(SpaceMismatch("`eigh!` requires domain and codomain to be the same")) S = spacetype(t) @@ -527,6 +513,7 @@ function eigh!(t::TensorMap; kwargs...) end function eig!(t::TensorMap; kwargs...) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:eig!) domain(t) == codomain(t) || throw(SpaceMismatch("`eig!` requires domain and codomain to be the same")) S = spacetype(t) diff --git a/src/tensors/linalg.jl b/src/tensors/linalg.jl index c66be7e4..c421c98d 100644 --- a/src/tensors/linalg.jl +++ b/src/tensors/linalg.jl @@ -99,23 +99,19 @@ for a specific isomorphism, but the current choice is such that `unitary(cod, dom) == inv(unitary(dom, cod)) = adjoint(unitary(dom, cod))`. """ function unitary(cod::TensorSpace{S}, dom::TensorSpace{S}) where {S} - InnerProductStyle(S) === EuclideanProduct() || - throw(ArgumentError("unitary requires inner product spaces")) + InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:unitary) return isomorphism(cod, dom) end function unitary(P::TensorMapSpace{S}) where {S} - InnerProductStyle(S) === EuclideanProduct() || - throw(ArgumentError("unitary requires inner product spaces")) + InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:unitary) return isomorphism(P) end function unitary(A::Type{<:DenseMatrix}, P::TensorMapSpace{S}) where {S} - InnerProductStyle(S) === EuclideanProduct() || - throw(ArgumentError("unitary requires inner product spaces")) + InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:unitary) return isomorphism(A, P) end function unitary(A::Type{<:DenseMatrix}, cod::TensorSpace{S}, dom::TensorSpace{S}) where {S} - InnerProductStyle(S) === EuclideanProduct() || - throw(ArgumentError("unitary requires inner product spaces")) + InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:unitary) return isomorphism(A, cod, dom) end @@ -139,8 +135,7 @@ end function isometry(::Type{A}, cod::ProductSpace{S}, dom::ProductSpace{S}) where {A<:DenseMatrix,S<:ElementarySpace} - InnerProductStyle(S) === EuclideanProduct() || - throw(ArgumentError("isometries require Euclidean inner product")) + InnerProductStyle(S) === EuclideanProduct() || throw_invalid_innerproduct(:isometry) dom ≾ cod || throw(SpaceMismatch("codomain $cod and domain $dom do not allow for an isometric mapping")) t = TensorMap(s -> A(undef, s), cod, dom) @@ -167,10 +162,9 @@ function Base.fill!(t::AbstractTensorMap, value::Number) end return t end -function LinearAlgebra.adjoint!(tdst::AbstractTensorMap, - tsrc::AbstractTensorMap) - spacetype(tdst) === spacetype(tsrc) && InnerProductStyle(tdst) === EuclideanProduct() || - throw(ArgumentError("adjoint! requires Euclidean inner product spacetype")) +function LinearAlgebra.adjoint!(tdst::AbstractTensorMap{S}, + tsrc::AbstractTensorMap{S}) where {S} + InnerProductStyle(tdst) === EuclideanProduct() || throw_invalid_innerproduct(:adjoint!) space(tdst) == adjoint(space(tsrc)) || throw(SpaceMismatch("$(space(tdst)) ≠ adjoint($(space(tsrc)))")) for c in blocksectors(tdst) @@ -207,8 +201,7 @@ end LinearAlgebra.dot(t1::AbstractTensorMap, t2::AbstractTensorMap) = inner(t1, t2) function LinearAlgebra.norm(t::AbstractTensorMap, p::Real=2) - InnerProductStyle(t) === EuclideanProduct() || - throw(ArgumentError("norm requires Euclidean inner product")) + InnerProductStyle(t) === EuclideanProduct() || throw_invalid_innerproduct(:norm) return _norm(blocks(t), p, float(zero(real(scalartype(t))))) end function _norm(blockiter, p::Real, init::Real) diff --git a/src/tensors/vectorinterface.jl b/src/tensors/vectorinterface.jl index 9462fc90..b2563352 100644 --- a/src/tensors/vectorinterface.jl +++ b/src/tensors/vectorinterface.jl @@ -80,8 +80,7 @@ end #------- function VectorInterface.inner(tx::AbstractTensorMap, ty::AbstractTensorMap) space(tx) == space(ty) || throw(SpaceMismatch("$(space(tx)) ≠ $(space(ty))")) - InnerProductStyle(tx) === EuclideanProduct() || - throw(ArgumentError("dot requires Euclidean inner product")) + InnerProductStyle(tx) === EuclideanProduct() || throw_invalid_innerproduct(:inner) T = VectorInterface.promote_inner(tx, ty) s = zero(T) for c in blocksectors(tx)