diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 1d26db06..3f40d042 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -23,7 +23,8 @@ export ZNSpace, Z2Space, Z3Space, Z4Space, U1Space, CU1Space, SU2Space export Vect, Rep # space constructors export CompositeSpace, ProductSpace # composite spaces export FusionTree -export IndexSpace, TensorSpace, AbstractTensorMap, AbstractTensor, TensorMap, Tensor # tensors and tensor properties +export IndexSpace, TensorSpace, TensorMapSpace +export AbstractTensorMap, AbstractTensor, TensorMap, Tensor, TrivialTensorMap # tensors and tensor properties export TruncationScheme export SpaceMismatch, SectorMismatch, IndexError # error types diff --git a/src/auxiliary/deprecate.jl b/src/auxiliary/deprecate.jl index 7e282c64..0a669fa2 100644 --- a/src/auxiliary/deprecate.jl +++ b/src/auxiliary/deprecate.jl @@ -1,6 +1,4 @@ -import Base: eltype, transpose -@deprecate eltype(T::Type{<:AbstractTensorMap}) scalartype(T) -@deprecate eltype(t::AbstractTensorMap) scalartype(t) +import Base: transpose #! format: off @deprecate permute(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; copy::Bool=false) permute(t, (p1, p2); copy=copy) diff --git a/src/fusiontrees/manipulations.jl b/src/fusiontrees/manipulations.jl index fb10b8ac..cea269fa 100644 --- a/src/fusiontrees/manipulations.jl +++ b/src/fusiontrees/manipulations.jl @@ -256,7 +256,7 @@ function bendright(f₁::FusionTree{I,N₁}, f₂::FusionTree{I,N₂}) where {I< coeff₀ = sqrtdim(c) * isqrtdim(a) if f₁.isdual[N₁] coeff₀ *= conj(frobeniusschur(dual(b))) - end + end if FusionStyle(I) isa MultiplicityFreeFusion coeff = coeff₀ * Bsymbol(a, b, c) vertices2 = N₂ > 0 ? (f₂.vertices..., nothing) : () diff --git a/src/spaces/cartesianspace.jl b/src/spaces/cartesianspace.jl index 86e145a4..bca3f854 100644 --- a/src/spaces/cartesianspace.jl +++ b/src/spaces/cartesianspace.jl @@ -1,11 +1,11 @@ """ - struct CartesianSpace <: ElementarySpace{ℝ} + struct CartesianSpace <: ElementarySpace A real Euclidean space `ℝ^d`, which is therefore self-dual. `CartesianSpace` has no additonal structure and is completely characterised by its dimension `d`. This is the vector space that is implicitly assumed in most of matrix algebra. """ -struct CartesianSpace <: ElementarySpace{ℝ} +struct CartesianSpace <: ElementarySpace d::Int end CartesianSpace(d::Integer=0; dual=false) = CartesianSpace(Int(d)) @@ -28,8 +28,10 @@ function CartesianSpace(dims::AbstractDict; kwargs...) end end +field(::Type{CartesianSpace}) = ℝ InnerProductStyle(::Type{CartesianSpace}) = EuclideanProduct() +Base.conj(V::CartesianSpace) = V isdual(V::CartesianSpace) = false # convenience constructor diff --git a/src/spaces/complexspace.jl b/src/spaces/complexspace.jl index 01a77df7..44e3fd02 100644 --- a/src/spaces/complexspace.jl +++ b/src/spaces/complexspace.jl @@ -1,11 +1,11 @@ """ - struct ComplexSpace <: ElementarySpace{ℂ} + struct ComplexSpace <: ElementarySpace A standard complex vector space ℂ^d with Euclidean inner product and no additional structure. It is completely characterised by its dimension and whether its the normal space or its dual (which is canonically isomorphic to the conjugate space). """ -struct ComplexSpace <: ElementarySpace{ℂ} +struct ComplexSpace <: ElementarySpace d::Int dual::Bool end @@ -29,6 +29,7 @@ function ComplexSpace(dims::AbstractDict; kwargs...) end end +field(::Type{ComplexSpace}) = ℂ InnerProductStyle(::Type{ComplexSpace}) = EuclideanProduct() # convenience constructor diff --git a/src/spaces/deligne.jl b/src/spaces/deligne.jl index f846895a..f9cc51c8 100644 --- a/src/spaces/deligne.jl +++ b/src/spaces/deligne.jl @@ -12,41 +12,50 @@ The Deligne tensor product also works in the type domain and for sectors and ten group representations, we have `Rep[G₁] ⊠ Rep[G₂] == Rep[G₁ × G₂]`, i.e. these are the natural representation spaces of the direct product of two groups. """ -⊠(V₁::VectorSpace, V₂::VectorSpace) = (V₁ ⊠ one(V₂)) ⊗ (one(V₁) ⊠ V₂) +function ⊠(V₁::VectorSpace, V₂::VectorSpace) + field(V₁) == field(V₂) || throw_incompatible_fields(V₁, V₂) + return (V₁ ⊠ one(V₂)) ⊗ (one(V₁) ⊠ V₂) +end -# define deligne products with empty tensor product: just add a trivial sector of the type of the empty space to each of the sectors in the non-empty space -function ⊠(V::GradedSpace, P₀::ProductSpace{<:ElementarySpace{ℂ},0}) +# define deligne products with empty tensor product: just add a trivial sector of the type +# of the empty space to each of the sectors in the non-empty space +function ⊠(V::GradedSpace, P₀::ProductSpace{<:ElementarySpace,0}) + field(V) == field(P₀) || throw_incompatible_fields(V, P₀) I₁ = sectortype(V) I₂ = sectortype(P₀) return Vect[I₁ ⊠ I₂](ifelse(isdual(V), dual(c), c) ⊠ one(I₂) => dim(V, c) for c in sectors(V); dual=isdual(V)) end -function ⊠(P₀::ProductSpace{<:ElementarySpace{ℂ},0}, V::GradedSpace) +function ⊠(P₀::ProductSpace{<:ElementarySpace,0}, V::GradedSpace) + field(P₀) == field(V) || throw_incompatible_fields(P₀, V) I₁ = sectortype(P₀) I₂ = sectortype(V) return Vect[I₁ ⊠ I₂](one(I₁) ⊠ ifelse(isdual(V), dual(c), c) => dim(V, c) for c in sectors(V); dual=isdual(V)) end -function ⊠(V::ComplexSpace, P₀::ProductSpace{<:ElementarySpace{ℂ},0}) +function ⊠(V::ComplexSpace, P₀::ProductSpace{<:ElementarySpace,0}) + field(V) == field(P₀) || throw_incompatible_fields(V, P₀) I₂ = sectortype(P₀) return Vect[I₂](one(I₂) => dim(V); dual=isdual(V)) end -function ⊠(P₀::ProductSpace{<:ElementarySpace{ℂ},0}, V::ComplexSpace) +function ⊠(P₀::ProductSpace{<:ElementarySpace,0}, V::ComplexSpace) + field(P₀) == field(V) || throw_incompatible_fields(P₀, V) I₁ = sectortype(P₀) return Vect[I₁](one(I₁) => dim(V); dual=isdual(V)) end -function ⊠(P::ProductSpace{<:ElementarySpace{ℂ},0}, - P₀::ProductSpace{<:ElementarySpace{ℂ},0}) +function ⊠(P::ProductSpace{<:ElementarySpace,0}, P₀::ProductSpace{<:ElementarySpace,0}) + field(P) == field(P₀) || throw_incompatible_fields(P, P₀) I₁ = sectortype(P) I₂ = sectortype(P₀) return one(Vect[I₁ ⊠ I₂]) end -function ⊠(P::ProductSpace{<:ElementarySpace{ℂ}}, P₀::ProductSpace{<:ElementarySpace{ℂ},0}) +function ⊠(P::ProductSpace{<:ElementarySpace}, P₀::ProductSpace{<:ElementarySpace,0}) + field(P) == field(P₀) || throw_incompatible_fields(P, P₀) I₁ = sectortype(P) I₂ = sectortype(P₀) S = Vect[I₁ ⊠ I₂] @@ -54,10 +63,15 @@ function ⊠(P::ProductSpace{<:ElementarySpace{ℂ}}, P₀::ProductSpace{<:Eleme return ProductSpace{S,N}(map(V -> V ⊠ P₀, tuple(P...))) end -function ⊠(P₀::ProductSpace{<:ElementarySpace{ℂ},0}, P::ProductSpace{<:ElementarySpace{ℂ}}) +function ⊠(P₀::ProductSpace{<:ElementarySpace,0}, P::ProductSpace{<:ElementarySpace}) + field(P₀) == field(P) || throw_incompatible_fields(P₀, P) I₁ = sectortype(P₀) I₂ = sectortype(P) S = Vect[I₁ ⊠ I₂] N = length(P) return ProductSpace{S,N}(map(V -> P₀ ⊠ V, tuple(P...))) end + +@noinline function throw_incompatible_fields(P₁, P₂) + throw(ArgumentError("Deligne products require spaces over the same field: $(field(P₁)) ≠ $(field(P₂))")) +end diff --git a/src/spaces/generalspace.jl b/src/spaces/generalspace.jl index 3f725683..8f9e9ab3 100644 --- a/src/spaces/generalspace.jl +++ b/src/spaces/generalspace.jl @@ -1,11 +1,11 @@ """ - struct GeneralSpace{𝕜} <: ElementarySpace{𝕜} + struct GeneralSpace{𝕜} <: ElementarySpace A finite-dimensional space over an arbitrary field `𝕜` without additional structure. It is thus characterized by its dimension, and whether or not it is the dual and/or conjugate space. For a real field `𝕜`, the space and its conjugate are the same. """ -struct GeneralSpace{𝕜} <: ElementarySpace{𝕜} +struct GeneralSpace{𝕜} <: ElementarySpace d::Int dual::Bool conj::Bool @@ -29,6 +29,7 @@ isconj(V::GeneralSpace) = V.conj Base.axes(V::GeneralSpace) = Base.OneTo(dim(V)) +field(::Type{GeneralSpace{𝕜}}) where {𝕜} = 𝕜 InnerProductStyle(::Type{<:GeneralSpace}) = NoInnerProduct() dual(V::GeneralSpace{𝕜}) where {𝕜} = GeneralSpace{𝕜}(dim(V), !isdual(V), isconj(V)) diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl index 214ffbea..b57fcfab 100644 --- a/src/spaces/gradedspace.jl +++ b/src/spaces/gradedspace.jl @@ -1,5 +1,5 @@ """ - struct GradedSpace{I<:Sector, D} <: ElementarySpace{ℂ} + struct GradedSpace{I<:Sector, D} <: ElementarySpace dims::D dual::Bool end @@ -24,7 +24,7 @@ and should typically be of no concern. The concrete type `GradedSpace{I,D}` with correct `D` can be obtained as `Vect[I]`, or if `I == Irrep[G]` for some `G<:Group`, as `Rep[G]`. """ -struct GradedSpace{I<:Sector,D} <: ElementarySpace{ℂ} +struct GradedSpace{I<:Sector,D} <: ElementarySpace dims::D dual::Bool end @@ -84,6 +84,7 @@ Base.hash(V::GradedSpace, h::UInt) = hash(V.dual, hash(V.dims, h)) # Corresponding methods: # properties +field(::Type{<:GradedSpace}) = ℂ InnerProductStyle(::Type{<:GradedSpace}) = EuclideanProduct() function dim(V::GradedSpace) return reduce(+, dim(V, c) * dim(c) for c in sectors(V); diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index 3a036f59..abcfb5ef 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -43,13 +43,14 @@ function Base.getindex(W::TensorMapSpace{<:IndexSpace,N₁,N₂}, i) where {N₁ return i <= N₁ ? codomain(W)[i] : dual(domain(W)[i - N₁]) end -function →(dom::TensorSpace{S}, codom::TensorSpace{S}) where {S<:ElementarySpace} - return HomSpace(ProductSpace(codom), ProductSpace(dom)) +function ←(codom::ProductSpace{S}, dom::ProductSpace{S}) where {S<:ElementarySpace} + return HomSpace(codom, dom) end - -function ←(codom::TensorSpace{S}, dom::TensorSpace{S}) where {S<:ElementarySpace} +function ←(codom::S, dom::S) where {S<:ElementarySpace} return HomSpace(ProductSpace(codom), ProductSpace(dom)) end +←(codom::VectorSpace, dom::VectorSpace) = ←(promote(codom, dom)...) +→(dom::VectorSpace, codom::VectorSpace) = ←(codom, dom) function Base.show(io::IO, W::HomSpace) if length(W.codomain) == 1 diff --git a/src/spaces/productspace.jl b/src/spaces/productspace.jl index bc5094fd..bb3d7dc2 100644 --- a/src/spaces/productspace.jl +++ b/src/spaces/productspace.jl @@ -177,21 +177,11 @@ Base.hash(P::ProductSpace, h::UInt) = hash(P.spaces, h) # Default construction from product of spaces #--------------------------------------------- -⊗(V₁::S, V₂::S) where {S<:ElementarySpace} = ProductSpace((V₁, V₂)) -function ⊗(P1::ProductSpace{S}, V₂::S) where {S<:ElementarySpace} - return ProductSpace(tuple(P1.spaces..., V₂)) -end -function ⊗(V₁::S, P2::ProductSpace{S}) where {S<:ElementarySpace} - return ProductSpace(tuple(V₁, P2.spaces...)) -end +⊗(V::Vararg{S}) where {S<:ElementarySpace} = ProductSpace(V) +⊗(P::ProductSpace) = P function ⊗(P1::ProductSpace{S}, P2::ProductSpace{S}) where {S<:ElementarySpace} - return ProductSpace(tuple(P1.spaces..., P2.spaces...)) + return ProductSpace{S}(tuple(P1.spaces..., P2.spaces...)) end -⊗(P::ProductSpace{S,0}, ::ProductSpace{S,0}) where {S<:ElementarySpace} = P -⊗(P::ProductSpace{S}, ::ProductSpace{S,0}) where {S<:ElementarySpace} = P -⊗(::ProductSpace{S,0}, P::ProductSpace{S}) where {S<:ElementarySpace} = P -⊗(V::ElementarySpace) = ProductSpace((V,)) -⊗(P::ProductSpace) = P # unit element with respect to the monoidal structure of taking tensor products """ @@ -205,14 +195,12 @@ Base.one(V::VectorSpace) = one(typeof(V)) Base.one(::Type{<:ProductSpace{S}}) where {S<:ElementarySpace} = ProductSpace{S,0}(()) Base.one(::Type{S}) where {S<:ElementarySpace} = ProductSpace{S,0}(()) -Base.convert(::Type{<:ProductSpace}, V::ElementarySpace) = ProductSpace((V,)) Base.:^(V::ElementarySpace, N::Int) = ProductSpace{typeof(V),N}(ntuple(n -> V, N)) Base.:^(V::ProductSpace, N::Int) = ⊗(ntuple(n -> V, N)...) function Base.literal_pow(::typeof(^), V::ElementarySpace, p::Val{N}) where {N} return ProductSpace{typeof(V),N}(ntuple(n -> V, p)) end -Base.convert(::Type{S}, P::ProductSpace{S,0}) where {S<:ElementarySpace} = oneunit(S) -Base.convert(::Type{S}, P::ProductSpace{S}) where {S<:ElementarySpace} = fuse(P.spaces...) + fuse(P::ProductSpace{S,0}) where {S<:ElementarySpace} = oneunit(S) fuse(P::ProductSpace{S}) where {S<:ElementarySpace} = fuse(P.spaces...) @@ -255,3 +243,18 @@ Base.eltype(P::ProductSpace) = eltype(typeof(P)) Base.IteratorEltype(::Type{<:ProductSpace}) = Base.HasEltype() Base.IteratorSize(::Type{<:ProductSpace}) = Base.HasLength() + +Base.reverse(P::ProductSpace) = ProductSpace(reverse(P.spaces)) + +# Promotion and conversion +# ------------------------ +function Base.promote_rule(::Type{S}, ::Type{<:ProductSpace{S}}) where {S<:ElementarySpace} + return ProductSpace{S} +end + +# ProductSpace to ElementarySpace +Base.convert(::Type{S}, P::ProductSpace{S,0}) where {S<:ElementarySpace} = oneunit(S) +Base.convert(::Type{S}, P::ProductSpace{S}) where {S<:ElementarySpace} = fuse(P.spaces...) + +# ElementarySpace to ProductSpace +Base.convert(::Type{<:ProductSpace}, V::S) where {S<:ElementarySpace} = ⊗(V) diff --git a/src/spaces/vectorspaces.jl b/src/spaces/vectorspaces.jl index b213a6c1..bc57dabc 100644 --- a/src/spaces/vectorspaces.jl +++ b/src/spaces/vectorspaces.jl @@ -87,10 +87,10 @@ function isdual end # Hierarchy of elementary vector spaces #--------------------------------------- """ - abstract type ElementarySpace{𝕜} <: VectorSpace end + abstract type ElementarySpace <: VectorSpace end -Elementary finite-dimensional vector space over a field `𝕜` that can be used as the index -space corresponding to the indices of a tensor. ElementarySpace is a super type for all +Elementary finite-dimensional vector space over a field that can be used as the index +space corresponding to the indices of a tensor. ElementarySpace is a supertype for all vector spaces (objects) that can be associated with the individual indices of a tensor, as hinted to by its alias IndexSpace. @@ -100,10 +100,11 @@ complex conjugate of the dual space is obtained as `dual(conj(V)) === conj(dual( different spaces should be of the same type, so that a tensor can be defined as an element of a homogeneous tensor product of these spaces. """ -abstract type ElementarySpace{𝕜} <: VectorSpace end +abstract type ElementarySpace <: VectorSpace end const IndexSpace = ElementarySpace -field(::Type{<:ElementarySpace{𝕜}}) where {𝕜} = 𝕜 +field(V::ElementarySpace) = field(typeof(V)) +# field(::Type{<:ElementarySpace{𝕜}}) where {𝕜} = 𝕜 """ oneunit(V::S) where {S<:ElementarySpace} -> S @@ -123,7 +124,8 @@ spaces `V₁`, `V₂`, ... Note that all the individual spaces should have the s [`isdual`](@ref), as otherwise the direct sum is not defined. """ function ⊕ end -⊕(V₁, V₂, V₃, V₄...) = ⊕(⊕(V₁, V₂), V₃, V₄...) +⊕(V₁::VectorSpace, V₂::VectorSpace) = ⊕(promote(V₁, V₂)...) +⊕(V::Vararg{VectorSpace}) = foldl(⊕, V) """ ⊗(V₁::S, V₂::S, V₃::S...) where {S<:ElementarySpace} -> S @@ -136,7 +138,8 @@ The tensor product structure is preserved, see [`fuse`](@ref) for returning a si elementary space of type `S` that is isomorphic to this tensor product. """ function ⊗ end -⊗(V₁, V₂, V₃, V₄...) = ⊗(⊗(V₁, V₂), V₃, V₄...) +⊗(V₁::VectorSpace, V₂::VectorSpace) = ⊗(promote(V₁, V₂)...) +⊗(V::Vararg{VectorSpace}) = foldl(⊗, V) # convenience definitions: Base.:*(V₁::VectorSpace, V₂::VectorSpace) = ⊗(V₁, V₂) @@ -171,7 +174,10 @@ Return the conjugate space of `V`. This should satisfy `conj(conj(V)) == V`. For `field(V)==ℝ`, `conj(V) == V`. It is assumed that `typeof(V) == typeof(conj(V))`. """ -Base.conj(V::ElementarySpace{ℝ}) = V +function Base.conj(V::ElementarySpace) + @assert field(V) == ℝ "default conj only defined for Vector spaces over ℝ" + return V +end # trait to describe the inner product type of vector spaces abstract type InnerProductStyle end @@ -192,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) @@ -229,7 +239,7 @@ end abstract type CompositeSpace{S<:ElementarySpace} <: VectorSpace end Abstract type for composite spaces that are defined in terms of a number of elementary -vector spaces of a homogeneous type `S<:ElementarySpace{𝕜}`. +vector spaces of a homogeneous type `S<:ElementarySpace`. """ abstract type CompositeSpace{S<:ElementarySpace} <: VectorSpace end diff --git a/src/tensors/braidingtensor.jl b/src/tensors/braidingtensor.jl index fb21289c..594e4feb 100644 --- a/src/tensors/braidingtensor.jl +++ b/src/tensors/braidingtensor.jl @@ -32,7 +32,11 @@ function BraidingTensor(V1::S, V2::S, adjoint::Bool=false) where {S<:IndexSpace} return BraidingTensor{S,Matrix{ComplexF64}}(V1, V2, adjoint) end end - +function BraidingTensor(V::HomSpace{S}, adjoint::Bool=false) where {S<:IndexSpace} + domain(V) == reverse(codomain(V)) || + throw(SpaceMismatch("Cannot define a braiding on $V")) + return BraidingTensor(V[1], V[2], adjoint) +end function Base.adjoint(b::BraidingTensor{S,A}) where {S<:IndexSpace,A<:DenseMatrix} return BraidingTensor{S,A}(b.V1, b.V2, !b.adjoint) end 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 ac059ce6..9f1820ca 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) diff --git a/test/spaces.jl b/test/spaces.jl index 8d02c97f..3634ef6c 100644 --- a/test/spaces.jl +++ b/test/spaces.jl @@ -119,8 +119,10 @@ println("------------------------------------") @test @constinferred(oneunit(V)) == W == oneunit(typeof(V)) @test @constinferred(⊕(V, V)) == ℂ^(2d) @test_throws SpaceMismatch (⊕(V, V')) - @test_throws MethodError (⊕(ℝ^d, ℂ^d)) - @test_throws MethodError (⊗(ℝ^d, ℂ^d)) + promote_except = ErrorException("promotion of types $(typeof(ℝ^d)) and " * + "$(typeof(ℂ^d)) failed to change any arguments") + @test_throws promote_except (⊕(ℝ^d, ℂ^d)) + @test_throws promote_except (⊗(ℝ^d, ℂ^d)) @test @constinferred(⊕(V, V)) == ℂ^(2d) @test @constinferred(⊕(V, oneunit(V))) == ℂ^(d + 1) @test @constinferred(⊕(V, V, V, V)) == ℂ^(4d)