diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md index 6309d0c03..2b5d6ca87 100644 --- a/docs/src/Changelog.md +++ b/docs/src/Changelog.md @@ -24,13 +24,13 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - Extended support for selecting storage types in the `TensorMap` constructors ([#327](https://github.com/QuantumKitHub/TensorKit.jl/pull/327)) - `similar_diagonal` to handle storage types when constructing diagonals ([#330](https://github.com/QuantumKitHub/TensorKit.jl/pull/330)) -- `LinearAlgebra.svd` overloads ### Fixed - Issue with using relative tolerances in truncation schemes ([#314](https://github.com/QuantumKitHub/TensorKit.jl/issues/314)) - Using `scalartype` instead of `eltype` in BLAS contraction ([#326](https://github.com/QuantumKitHub/TensorKit.jl/pull/326)) - Divide by zero error in `show` for empty tensors ([#329](https://github.com/QuantumKitHub/TensorKit.jl/pull/329)) +- `svd_vals(::DiagonalTensorMap)` correctly outputs `SectorVector` and implementation fix. ([#333](https://github.com/QuantumKitHub/TensorKit.jl/pull/329)) ## [0.16.0](https://github.com/QuantumKitHub/TensorKit.jl/releases/tag/v0.16.0) - 2025-12-08 diff --git a/src/factorizations/diagonal.jl b/src/factorizations/diagonal.jl index a101b1f21..cec4a9061 100644 --- a/src/factorizations/diagonal.jl +++ b/src/factorizations/diagonal.jl @@ -1,6 +1,7 @@ # DiagonalTensorMap # ----------------- _repack_diagonal(d::DiagonalTensorMap) = Diagonal(d.data) +_repack_diagonal(d::SectorVector) = Diagonal(parent(d)) MAK.diagview(t::DiagonalTensorMap) = SectorVector(t.data, TensorKit.diagonalblockstructure(space(t))) @@ -93,17 +94,10 @@ function MAK.svd_compact!(t::AbstractTensorMap, USVᴴ, alg::DiagonalAlgorithm) return svd_full!(t, USVᴴ, alg) end -# f_vals -# ------ -for f! in (:eig_vals!, :eigh_vals!, :svd_vals!) - @eval function MAK.$f!(d::AbstractTensorMap, V, alg::DiagonalAlgorithm) - $f!(_repack_diagonal(d), diagview(_repack_diagonal(V)), alg) - return V - end - @eval function MAK.initialize_output( - ::typeof($f!), d::DiagonalTensorMap, alg::DiagonalAlgorithm - ) - data = MAK.initialize_output($f!, _repack_diagonal(d), alg) - return DiagonalTensorMap(data, d.domain) - end +# For diagonal inputs we don't have to promote the scalartype since we know they are symmetric +function MAK.initialize_output(::typeof(eig_vals!), t::AbstractTensorMap, alg::DiagonalAlgorithm) + V_D = fuse(domain(t)) + Tc = scalartype(t) + A = similarstoragetype(t, Tc) + return SectorVector{Tc, sectortype(t), A}(undef, V_D) end diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 6a23ef045..0c2272ac8 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -49,11 +49,6 @@ function LinearAlgebra.eigvals(t::AbstractTensorMap; kwargs...) end LinearAlgebra.eigvals!(t::AbstractTensorMap; kwargs...) = eig_vals!(t) -LinearAlgebra.svd(t::AbstractTensorMap; full::Bool = false) = - full ? svd_full(t) : svd_compact(t) -LinearAlgebra.svd!(t::AbstractTensorMap; full::Bool = false) = - full ? svd_full!(t) : svd_compact!(t) - function LinearAlgebra.svdvals(t::AbstractTensorMap) tcopy = copy_oftype(t, factorisation_scalartype(svd_vals!, t)) return LinearAlgebra.svdvals!(tcopy) diff --git a/test/tensors/factorizations.jl b/test/tensors/factorizations.jl index b017b746b..176d62657 100644 --- a/test/tensors/factorizations.jl +++ b/test/tensors/factorizations.jl @@ -1,6 +1,7 @@ using Test, TestExtras using TensorKit using LinearAlgebra: LinearAlgebra +using MatrixAlgebraKit: diagview @isdefined(TestSetup) || include("../setup.jl") using .TestSetup @@ -191,10 +192,9 @@ for V in spacelist @test isposdef(s) @test isisometric(vᴴ; side = :right) - s′ = LinearAlgebra.diag(s) - for (c, b) in pairs(LinearAlgebra.svdvals(t)) - @test b ≈ s′[c] - end + s′ = @constinferred svd_vals(t) + @test s′ ≈ diagview(s) + @test s′ isa TensorKit.SectorVector v, c = @constinferred left_orth(t; alg = :svd) @test v * c ≈ t @@ -261,14 +261,14 @@ for V in spacelist @test norm(t - U1 * S1 * Vᴴ1) ≈ ϵ1 atol = eps(real(T))^(4 / 5) @test dim(domain(S1)) <= nvals - λ = minimum(minimum, values(LinearAlgebra.diag(S1))) + λ = minimum(diagview(S1)) trunc = trunctol(; atol = λ - 10eps(λ)) U2, S2, Vᴴ2, ϵ2 = @constinferred svd_trunc(t; trunc) @test t * Vᴴ2' ≈ U2 * S2 @test isisometric(U2) @test isisometric(Vᴴ2; side = :right) @test norm(t - U2 * S2 * Vᴴ2) ≈ ϵ2 atol = eps(real(T))^(4 / 5) - @test minimum(minimum, values(LinearAlgebra.diag(S1))) >= λ + @test minimum(diagview(S1)) >= λ @test U2 ≈ U1 @test S2 ≈ S1 @test Vᴴ2 ≈ Vᴴ1 @@ -297,7 +297,7 @@ for V in spacelist @test isisometric(U5) @test isisometric(Vᴴ5; side = :right) @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5) - @test minimum(minimum, values(LinearAlgebra.diag(S5))) >= λ + @test minimum(diagview(S5)) >= λ @test dim(domain(S5)) ≤ nvals end end @@ -312,13 +312,11 @@ for V in spacelist d, v = @constinferred eig_full(t) @test t * v ≈ v * d - d′ = LinearAlgebra.diag(d) - for (c, b) in pairs(LinearAlgebra.eigvals(t)) - @test sort(b; by = abs) ≈ sort(d′[c]; by = abs) - end + d′ = @constinferred eig_vals(t) + @test d′ ≈ diagview(d) + @test d′ isa TensorKit.SectorVector - vdv = v' * v - vdv = (vdv + vdv') / 2 + vdv = project_hermitian!(v' * v) @test @constinferred isposdef(vdv) t isa DiagonalTensorMap || @test !isposdef(t) # unlikely for non-hermitian map @@ -327,35 +325,34 @@ for V in spacelist @test t * v ≈ v * d @test dim(domain(d)) ≤ nvals - t2 = (t + t') + t2 = @constinferred project_hermitian(t) D, V = eigen(t2) @test isisometric(V) D̃, Ṽ = @constinferred eigh_full(t2) @test D ≈ D̃ @test V ≈ Ṽ - λ = minimum( - minimum(real(LinearAlgebra.diag(b))) - for (c, b) in blocks(D) - ) + λ = minimum(real, diagview(D)) @test cond(Ṽ) ≈ one(real(T)) @test isposdef(t2) == isposdef(λ) @test isposdef(t2 - λ * one(t2) + 0.1 * one(t2)) @test !isposdef(t2 - λ * one(t2) - 0.1 * one(t2)) - add!(t, t') - - d, v = @constinferred eigh_full(t) - @test t * v ≈ v * d + d, v = @constinferred eigh_full(t2) + @test t2 * v ≈ v * d @test isunitary(v) - λ = minimum(minimum(real(LinearAlgebra.diag(b))) for (c, b) in blocks(d)) + d′ = @constinferred eigh_vals(t2) + @test d′ ≈ diagview(d) + @test d′ isa TensorKit.SectorVector + + λ = minimum(real, diagview(d)) @test cond(v) ≈ one(real(T)) - @test isposdef(t) == isposdef(λ) - @test isposdef(t - λ * one(t) + 0.1 * one(t)) - @test !isposdef(t - λ * one(t) - 0.1 * one(t)) + @test isposdef(t2) == isposdef(λ) + @test isposdef(t2 - λ * one(t) + 0.1 * one(t2)) + @test !isposdef(t2 - λ * one(t) - 0.1 * one(t2)) - d, v = @constinferred eigh_trunc(t; trunc = truncrank(nvals)) - @test t * v ≈ v * d + d, v = @constinferred eigh_trunc(t2; trunc = truncrank(nvals)) + @test t2 * v ≈ v * d @test dim(domain(d)) ≤ nvals end end @@ -390,7 +387,7 @@ for V in spacelist @test cond(t2) == 0.0 end for T in eltypes, t in (rand(T, W, W), rand(T, W, W)') - add!(t, t') + project_hermitian!(t) vals = @constinferred LinearAlgebra.eigvals(t) λmax = maximum(s -> maximum(abs, s), values(vals)) λmin = minimum(s -> minimum(abs, s), values(vals))