From ce5aae5035b0286c80b41bf52d5dc31d6fc0dc17 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 15 Dec 2025 13:14:56 -0500 Subject: [PATCH 1/9] add `similar_diagonal` --- src/tensors/abstracttensor.jl | 66 +++++++++++++++++++++++++++++++++++ src/tensors/diagonal.jl | 10 +++--- 2 files changed, 72 insertions(+), 4 deletions(-) diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl index 7eeb3af9d..5e5b2b83f 100644 --- a/src/tensors/abstracttensor.jl +++ b/src/tensors/abstracttensor.jl @@ -491,6 +491,8 @@ The structure may be specified either as a single `HomSpace` argument or as `cod By default, this will result in `TensorMap{T}(undef, V)` when custom objects do not specialize this method. + +See also [`similar_diagonal`](@ref). """ Base.similar(::AbstractTensorMap, args...) function Base.similar( @@ -543,6 +545,70 @@ function Base.similar( return TensorMap{scalartype(TT)}(undef, cod, dom) end +# similar diagonal +# ---------------- +# The implementation is again written for similar_diagonal(t, TorA, V::TensorSpace) -> DiagonalTensorMap +# and all other methods are just filling in default arguments +@doc """ + similar_diagonal(t::AbstractTensorMap, [AorT=storagetype(t)], [V=space(t)]) + similar_diagonal(t::AbstractTensorMap, [AorT=storagetype(t)], codomain, domain) + +Creates an uninitialized mutable diagonal tensor with the given scalar or storagetype `AorT` and +structure `V` or `codomain ← domain`, based on the source tensormap. The second and third +arguments are both optional, defaulting to the given tensor's `storagetype` and `space`. +The structure may be specified either as a single `HomSpace` argument or as `codomain` and +`domain`. + +By default, this will result in `DiagonalTensorMap{T}(undef, V)` when custom objects do not +specialize this method. Furthermore, the method will throw if the provided space is not compatible +with a diagonal structure. + +See also [`Base.similar`](@ref). +""" similar_diagonal(::AbstractTensorMap, args...) + +# 4 arguments +function similar_diagonal(t::AbstractTensorMap, ::Type{T}, cod::TensorSpace, dom::TensorSpace) where {T} + length(cod) == length(dom) == 1 && cod == dom || + throw(ArgumentError("requested space is not square")) + return similar_diagonal(t, T, cod) +end + +# 3 arguments +function similar_diagonal(t::AbstractTensorMap, ::Type{T}, V::TensorMapSpace) where {T} + numout(V) == numin(V) == 1 && domain(V) == codomain(V) || + throw(ArgumentError("requested space is not square")) + return similar_diagonal(t, T, codomain(V)) +end +function similar_diagonal(t::AbstractTensorMap, ::Type{T}, V::ProductSpace) where {T} + length(V) == 1 || throw(ArgumentError()) + return similar_diagonal(t, T, only(V.spaces)) +end +function similar_diagonal(t::AbstractTensorMap, ::Type{TorA}, V::ElementarySpace) where {TorA} + if TorA <: Number + T = TorA + A = similarstoragetype(t, T) + elseif TorA <: DenseVector + A = TorA + T = scalartype(A) + else + throw(ArgumentError("Type $TorA not supported for similar")) + end + + return DiagonalTensorMap{T, spacetype(V), A}(undef, V) +end + +# 2 arguments +similar_diagonal(t::AbstractTensorMap, ::Type{T}) where {T} = + similar_diagonal(t, T, space(t)) +similar_diagonal(t::AbstractTensorMap, P::TensorMapSpace) = + similar_diagonal(t, similarstoragetype(t), P) +similar_diagonal(t::AbstractTensorMap, P::TensorSpace) = + similar_diagonal(t, similarstoragetype(t), P) + +# 1 argument +similar_diagonal(t::AbstractTensorMap) = + similar_diagonal(t, similarstoragetype(t), space(t)) + # Equality and approximality #---------------------------- function Base.:(==)(t1::AbstractTensorMap, t2::AbstractTensorMap) diff --git a/src/tensors/diagonal.jl b/src/tensors/diagonal.jl index 356f5dbf4..e35c2ffdd 100644 --- a/src/tensors/diagonal.jl +++ b/src/tensors/diagonal.jl @@ -78,10 +78,12 @@ function DiagonalTensorMap(t::AbstractTensorMap{T, S, 1, 1}) where {T, S} return d end -Base.similar(d::DiagonalTensorMap) = DiagonalTensorMap(similar(d.data), d.domain) -function Base.similar(d::DiagonalTensorMap, ::Type{T}) where {T <: Number} - return DiagonalTensorMap(similar(d.data, T), d.domain) -end +Base.similar(d::DiagonalTensorMap) = similar_diagonal(d) +Base.similar(d::DiagonalTensorMap, ::Type{T}) where {T} = similar_diagonal(d, T) + +similar_diagonal(d::DiagonalTensorMap) = DiagonalTensorMap(similar(d.data), d.domain) +similar_diagonal(d::DiagonalTensorMap, ::Type{T}) where {T <: Number} = + DiagonalTensorMap(similar(d.data, T), d.domain) # TODO: more constructors needed? From 69e6055958ff4ac00c0a979c574a992b9c79cff5 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 15 Dec 2025 13:24:56 -0500 Subject: [PATCH 2/9] update factorization --- src/factorizations/diagonal.jl | 2 +- src/factorizations/factorizations.jl | 2 +- src/factorizations/matrixalgebrakit.jl | 7 +++---- src/factorizations/truncation.jl | 4 ++-- 4 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/factorizations/diagonal.jl b/src/factorizations/diagonal.jl index c54bef00e..5784395cb 100644 --- a/src/factorizations/diagonal.jl +++ b/src/factorizations/diagonal.jl @@ -72,7 +72,7 @@ function MAK.initialize_output( V_cod = fuse(codomain(t)) V_dom = fuse(domain(t)) U = similar(t, codomain(t) ← V_cod) - S = DiagonalTensorMap{real(scalartype(t))}(undef, V_cod ← V_dom) + S = similar_diagonal(t, real(scalartype(t)), V_cod ← V_dom) Vᴴ = similar(t, V_dom ← domain(t)) return U, S, Vᴴ end diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 73fbc45eb..78dedd327 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -6,7 +6,7 @@ module Factorizations export copy_oftype, factorisation_scalartype, one!, truncspace using ..TensorKit -using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, blocktype, foreachblock, one! +using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, blocktype, foreachblock, one!, similar_diagonal using LinearAlgebra: LinearAlgebra, BlasFloat, Diagonal, svdvals, svdvals!, eigen, eigen!, isposdef, isposdef!, ishermitian diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index 4564a8137..d8ad2c0c2 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -74,7 +74,7 @@ end function MAK.initialize_output(::typeof(svd_compact!), t::AbstractTensorMap, ::AbstractAlgorithm) V_cod = V_dom = infimum(fuse(codomain(t)), fuse(domain(t))) U = similar(t, codomain(t) ← V_cod) - S = DiagonalTensorMap{real(scalartype(t))}(undef, V_cod) + S = similar_diagonal(t, real(scalartype(t)), V_cod) Vᴴ = similar(t, V_dom ← domain(t)) return U, S, Vᴴ end @@ -89,8 +89,7 @@ end # ------------------------ function MAK.initialize_output(::typeof(eigh_full!), t::AbstractTensorMap, ::AbstractAlgorithm) V_D = fuse(domain(t)) - T = real(scalartype(t)) - D = DiagonalTensorMap{T}(undef, V_D) + D = similar_diagonal(t, real(scalartype(t)), V_D) V = similar(t, codomain(t) ← V_D) return D, V end @@ -98,7 +97,7 @@ end function MAK.initialize_output(::typeof(eig_full!), t::AbstractTensorMap, ::AbstractAlgorithm) V_D = fuse(domain(t)) Tc = complex(scalartype(t)) - D = DiagonalTensorMap{Tc}(undef, V_D) + D = similar_diagonal(t, Tc, V_D) V = similar(t, Tc, codomain(t) ← V_D) return D, V end diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index dd2e663c6..8e1899616 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -67,7 +67,7 @@ function MAK.truncate( Ũ = similar(U, codomain(U) ← V_truncated) truncate_domain!(Ũ, U, ind) - S̃ = DiagonalTensorMap{scalartype(S)}(undef, V_truncated) + S̃ = similar_diagonal(S, V_truncated) truncate_diagonal!(S̃, S, ind) Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ)) truncate_codomain!(Ṽᴴ, Vᴴ, ind) @@ -132,7 +132,7 @@ for f! in (:eig_trunc!, :eigh_trunc!) ind = MAK.findtruncated(diagview(D), strategy) V_truncated = truncate_space(space(D, 1), ind) - D̃ = DiagonalTensorMap{scalartype(D)}(undef, V_truncated) + D = similar_diagonal(D, V_truncated) truncate_diagonal!(D̃, D, ind) Ṽ = similar(V, codomain(V) ← V_truncated) From 32571c8a60c4e1a151066968d983d38daf394367 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 15 Dec 2025 13:21:35 -0500 Subject: [PATCH 3/9] also update sectorvector --- src/factorizations/matrixalgebrakit.jl | 9 ++++++--- src/tensors/sectorvector.jl | 5 +++++ 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/factorizations/matrixalgebrakit.jl b/src/factorizations/matrixalgebrakit.jl index d8ad2c0c2..90776892c 100644 --- a/src/factorizations/matrixalgebrakit.jl +++ b/src/factorizations/matrixalgebrakit.jl @@ -82,7 +82,8 @@ end function MAK.initialize_output(::typeof(svd_vals!), t::AbstractTensorMap, alg::AbstractAlgorithm) V_cod = infimum(fuse(codomain(t)), fuse(domain(t))) T = real(scalartype(t)) - return SectorVector{T}(undef, V_cod) + A = similarstoragetype(t, T) + return SectorVector{T, sectortype(t), A}(undef, V_cod) end # Eigenvalue decomposition @@ -105,13 +106,15 @@ end function MAK.initialize_output(::typeof(eigh_vals!), t::AbstractTensorMap, alg::AbstractAlgorithm) V_D = fuse(domain(t)) T = real(scalartype(t)) - return SectorVector{T}(undef, V_D) + A = similarstoragetype(t, T) + return SectorVector{T, sectortype(t), A}(undef, V_D) end function MAK.initialize_output(::typeof(eig_vals!), t::AbstractTensorMap, alg::AbstractAlgorithm) V_D = fuse(domain(t)) Tc = complex(scalartype(t)) - return SectorVector{Tc}(undef, V_D) + A = similarstoragetype(t, Tc) + return SectorVector{Tc, sectortype(t), A}(undef, V_D) end # QR decomposition diff --git a/src/tensors/sectorvector.jl b/src/tensors/sectorvector.jl index a2a081058..f28642705 100644 --- a/src/tensors/sectorvector.jl +++ b/src/tensors/sectorvector.jl @@ -16,6 +16,11 @@ function SectorVector{T}(::UndefInitializer, V::ElementarySpace) where {T} structure = diagonalblockstructure(V ← V) return SectorVector(data, structure) end +function SectorVector{T, I, A}(::UndefInitializer, V::ElementarySpace) where {T, I, A} + data = A(undef, reduceddim(V)) + structure = diagonalblockstructure(V ← V) + return SectorVector{T, I, A}(data, structure) +end Base.parent(v::SectorVector) = v.data From c6b1eb92d525480bc87bdcbb3c5c0790097e74e1 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 16 Dec 2025 07:14:17 -0500 Subject: [PATCH 4/9] remove ambiguity --- src/tensors/sectorvector.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tensors/sectorvector.jl b/src/tensors/sectorvector.jl index f28642705..4c914d6de 100644 --- a/src/tensors/sectorvector.jl +++ b/src/tensors/sectorvector.jl @@ -16,7 +16,7 @@ function SectorVector{T}(::UndefInitializer, V::ElementarySpace) where {T} structure = diagonalblockstructure(V ← V) return SectorVector(data, structure) end -function SectorVector{T, I, A}(::UndefInitializer, V::ElementarySpace) where {T, I, A} +function SectorVector{T, I, A}(::UndefInitializer, V::ElementarySpace) where {T, I, A <: AbstractVector{T}} data = A(undef, reduceddim(V)) structure = diagonalblockstructure(V ← V) return SectorVector{T, I, A}(data, structure) From 85f5758752916d9cb0275d133afe4fa7a0cef6e5 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 16 Dec 2025 17:35:44 -0500 Subject: [PATCH 5/9] Update src/tensors/abstracttensor.jl Co-authored-by: Jutho --- src/tensors/abstracttensor.jl | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl index 5e5b2b83f..454dfc560 100644 --- a/src/tensors/abstracttensor.jl +++ b/src/tensors/abstracttensor.jl @@ -597,17 +597,16 @@ function similar_diagonal(t::AbstractTensorMap, ::Type{TorA}, V::ElementarySpace return DiagonalTensorMap{T, spacetype(V), A}(undef, V) end -# 2 arguments -similar_diagonal(t::AbstractTensorMap, ::Type{T}) where {T} = - similar_diagonal(t, T, space(t)) -similar_diagonal(t::AbstractTensorMap, P::TensorMapSpace) = - similar_diagonal(t, similarstoragetype(t), P) -similar_diagonal(t::AbstractTensorMap, P::TensorSpace) = - similar_diagonal(t, similarstoragetype(t), P) +similar_diagonal(t::AbstractTensorMap) = similar_diagonal(t, similarstoragetype(t), _diagspace(t)) +similar_diagonal(t::AbstractTensorMap, V::ElementarySpace) = similar_diagonal(t, similarstoragetype(t), V) +similar_diagonal(t::AbstractTensorMap, T::Type) = similar_diagonal(t, T, _diagspace(t)) -# 1 argument -similar_diagonal(t::AbstractTensorMap) = - similar_diagonal(t, similarstoragetype(t), space(t)) +function _diagspace(t) + cod, dom = codomain(t), domain(t) + length(cod) == 1 && cod == dom || + throw(ArgumentError("space does not support a diagonal TensorMap")) + return only(cod) +end # Equality and approximality #---------------------------- From 637e890dbdc8c159956119bff3ecaca4298eea56 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 16 Dec 2025 17:34:44 -0500 Subject: [PATCH 6/9] fix import --- src/factorizations/factorizations.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 78dedd327..02f1712eb 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -6,9 +6,12 @@ module Factorizations export copy_oftype, factorisation_scalartype, one!, truncspace using ..TensorKit -using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, blocktype, foreachblock, one!, similar_diagonal +using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, + blocktype, foreachblock, one!, + similar_diagonal, similarstoragetype -using LinearAlgebra: LinearAlgebra, BlasFloat, Diagonal, svdvals, svdvals!, eigen, eigen!, +using LinearAlgebra: LinearAlgebra, BlasFloat, Diagonal, + svdvals, svdvals!, eigen, eigen!, isposdef, isposdef!, ishermitian using TensorOperations: Index2Tuple From a0f61532c405d33eafb7e2308dc81d8404b02294 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 16 Dec 2025 17:41:56 -0500 Subject: [PATCH 7/9] purge some more constructors --- src/factorizations/diagonal.jl | 2 +- src/tensors/abstracttensor.jl | 28 +++++----------------------- 2 files changed, 6 insertions(+), 24 deletions(-) diff --git a/src/factorizations/diagonal.jl b/src/factorizations/diagonal.jl index 5784395cb..a101b1f21 100644 --- a/src/factorizations/diagonal.jl +++ b/src/factorizations/diagonal.jl @@ -72,7 +72,7 @@ function MAK.initialize_output( V_cod = fuse(codomain(t)) V_dom = fuse(domain(t)) U = similar(t, codomain(t) ← V_cod) - S = similar_diagonal(t, real(scalartype(t)), V_cod ← V_dom) + S = similar_diagonal(t, real(scalartype(t)), V_cod) Vᴴ = similar(t, V_dom ← domain(t)) return U, S, Vᴴ end diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl index 454dfc560..737a875d5 100644 --- a/src/tensors/abstracttensor.jl +++ b/src/tensors/abstracttensor.jl @@ -547,17 +547,15 @@ end # similar diagonal # ---------------- -# The implementation is again written for similar_diagonal(t, TorA, V::TensorSpace) -> DiagonalTensorMap +# The implementation is again written for similar_diagonal(t, TorA, V::ElementarySpace) -> DiagonalTensorMap # and all other methods are just filling in default arguments @doc """ - similar_diagonal(t::AbstractTensorMap, [AorT=storagetype(t)], [V=space(t)]) - similar_diagonal(t::AbstractTensorMap, [AorT=storagetype(t)], codomain, domain) + similar_diagonal(t::AbstractTensorMap, [AorT=storagetype(t)], [V::ElementarySpace]) Creates an uninitialized mutable diagonal tensor with the given scalar or storagetype `AorT` and -structure `V` or `codomain ← domain`, based on the source tensormap. The second and third -arguments are both optional, defaulting to the given tensor's `storagetype` and `space`. -The structure may be specified either as a single `HomSpace` argument or as `codomain` and -`domain`. +structure `V ← V`, based on the source tensormap. The second argument is optional and defaults +to the given tensor's `storagetype`, while the third argument can only be omitted for square +input tensors of space `V ← V`, to conform with the diagonal structure. By default, this will result in `DiagonalTensorMap{T}(undef, V)` when custom objects do not specialize this method. Furthermore, the method will throw if the provided space is not compatible @@ -566,23 +564,7 @@ with a diagonal structure. See also [`Base.similar`](@ref). """ similar_diagonal(::AbstractTensorMap, args...) -# 4 arguments -function similar_diagonal(t::AbstractTensorMap, ::Type{T}, cod::TensorSpace, dom::TensorSpace) where {T} - length(cod) == length(dom) == 1 && cod == dom || - throw(ArgumentError("requested space is not square")) - return similar_diagonal(t, T, cod) -end - # 3 arguments -function similar_diagonal(t::AbstractTensorMap, ::Type{T}, V::TensorMapSpace) where {T} - numout(V) == numin(V) == 1 && domain(V) == codomain(V) || - throw(ArgumentError("requested space is not square")) - return similar_diagonal(t, T, codomain(V)) -end -function similar_diagonal(t::AbstractTensorMap, ::Type{T}, V::ProductSpace) where {T} - length(V) == 1 || throw(ArgumentError()) - return similar_diagonal(t, T, only(V.spaces)) -end function similar_diagonal(t::AbstractTensorMap, ::Type{TorA}, V::ElementarySpace) where {TorA} if TorA <: Number T = TorA From 4a00c134bbc060dc9226aa319bdd63c96f744063 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 16 Dec 2025 20:14:52 -0500 Subject: [PATCH 8/9] Update src/tensors/abstracttensor.jl Co-authored-by: Jutho --- src/tensors/abstracttensor.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tensors/abstracttensor.jl b/src/tensors/abstracttensor.jl index 737a875d5..217429a4e 100644 --- a/src/tensors/abstracttensor.jl +++ b/src/tensors/abstracttensor.jl @@ -586,7 +586,7 @@ similar_diagonal(t::AbstractTensorMap, T::Type) = similar_diagonal(t, T, _diagsp function _diagspace(t) cod, dom = codomain(t), domain(t) length(cod) == 1 && cod == dom || - throw(ArgumentError("space does not support a diagonal TensorMap")) + throw(ArgumentError("space does not support a DiagonalTensorMap")) return only(cod) end From ef0d6ec83f845e485d513b344dfe1e843c48c51a Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 16 Dec 2025 20:16:54 -0500 Subject: [PATCH 9/9] small oopsie --- src/factorizations/truncation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 8e1899616..42d1e3248 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -132,7 +132,7 @@ for f! in (:eig_trunc!, :eigh_trunc!) ind = MAK.findtruncated(diagview(D), strategy) V_truncated = truncate_space(space(D, 1), ind) - D = similar_diagonal(D, V_truncated) + D̃ = similar_diagonal(D, V_truncated) truncate_diagonal!(D̃, D, ind) Ṽ = similar(V, codomain(V) ← V_truncated)