diff --git a/Project.toml b/Project.toml index a759bdaef..e8d17e10f 100644 --- a/Project.toml +++ b/Project.toml @@ -57,7 +57,7 @@ FiniteDifferences = "0.12" GPUArrays = "11.4.1" LRUCache = "1.6" LinearAlgebra = "1" -MatrixAlgebraKit = "0.6.9" +MatrixAlgebraKit = "0.6.13" Mooncake = "0.5.27" OhMyThreads = "0.8.0" PrecompileTools = "1" diff --git a/ext/TensorKitGPUArraysExt.jl b/ext/TensorKitGPUArraysExt.jl index 5e013865d..45f6df640 100644 --- a/ext/TensorKitGPUArraysExt.jl +++ b/ext/TensorKitGPUArraysExt.jl @@ -85,7 +85,7 @@ function MatrixAlgebraKit.findtruncated_svd(values::GPUSectorVector, strategy::S return Adapt.adapt(Vector, MatrixAlgebraKit.findtruncated(values, strategy)) end -for strat in (:(MatrixAlgebraKit.TruncationByOrder), :(MatrixAlgebraKit.TruncationByError), :(MatrixAlgebraKit.TruncationIntersection), :(TensorKit.Factorizations.TruncationSpace)) +for strat in (:(MatrixAlgebraKit.TruncationByOrder), :(MatrixAlgebraKit.TruncationByError), :(MatrixAlgebraKit.TruncationIntersection), :(MatrixAlgebraKit.TruncationMultiplets), :(TensorKit.Factorizations.TruncationSpace)) @eval function MatrixAlgebraKit.findtruncated_svd(values::GPUSectorVector, strategy::$strat) # returning a GPUSectorVector wrecks things in truncate_{co}domain # because of scalar indexing diff --git a/src/TensorKit.jl b/src/TensorKit.jl index 6bdbebfb5..2d1d806db 100644 --- a/src/TensorKit.jl +++ b/src/TensorKit.jl @@ -99,7 +99,7 @@ export @tensor, @tensoropt, @ncon, ncon, @planar, @plansor export scalar, add!, contract! # truncation schemes -export notrunc, truncrank, trunctol, truncfilter, truncspace, truncerror +export notrunc, truncrank, trunctol, truncfilter, truncspace, truncerror, truncmultiplets # cache management export empty_globalcaches! diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl index 7bada491c..a1753fa6a 100644 --- a/src/factorizations/factorizations.jl +++ b/src/factorizations/factorizations.jl @@ -23,7 +23,8 @@ using MatrixAlgebraKit import MatrixAlgebraKit as MAK using MatrixAlgebraKit: AbstractAlgorithm, TruncatedAlgorithm, DiagonalAlgorithm using MatrixAlgebraKit: TruncationStrategy, NoTruncation, TruncationByValue, TruncationByError, - TruncationIntersection, TruncationUnion, TruncationByFilter, TruncationByOrder + TruncationIntersection, TruncationUnion, TruncationByFilter, TruncationByOrder, + TruncationMultiplets using MatrixAlgebraKit: diagview include("utility.jl") diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index bd6956440..426e14927 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -337,6 +337,37 @@ function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationUnion) return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals) end +# Degenerate values can occur in different blocks, for example when the sectors are related by +# a symmetry, so groups of degenerate values are determined across all blocks simultaneously. +function MAK.findtruncated(values::SectorVector, strategy::TruncationMultiplets) + return _findtruncated_multiplets(values, MAK.findtruncated(values, strategy.trunc), strategy) +end +function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationMultiplets) + return _findtruncated_multiplets(values, MAK.findtruncated_svd(values, strategy.trunc), strategy) +end + +function _findtruncated_multiplets(values::SectorVector, inds, strategy::TruncationMultiplets) + keep = similar(values, Bool) + to_keep = collect(parent(_ind_to_mask!(keep, inds))) + kept_values = collect(parent(values)) + perm = sortperm(kept_values; by = abs, rev = true) + MAK._keep_multiplets!(to_keep, kept_values, perm, strategy) + copyto!(parent(keep), to_keep) + return keep +end + +function _ind_to_mask!(keep::SectorVector{Bool}, inds::SectorVector{Bool}) + copyto!(parent(keep), parent(inds)) + return keep +end +function _ind_to_mask!(keep::SectorVector{Bool}, inds::SectorDict) + fill!(parent(keep), false) + for (c, ind) in pairs(inds) + keep[c][ind] .= true + end + return keep +end + # Truncation error # ---------------- MAK.truncation_error(values::SectorVector, ind) = MAK.truncation_error!(copy(values), ind) diff --git a/test/amd/factorizations.jl b/test/amd/factorizations.jl index b65823964..d45d165d1 100644 --- a/test/amd/factorizations.jl +++ b/test/amd/factorizations.jl @@ -434,6 +434,31 @@ for V in spacelist @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5) @test minimum(diagview(S5)) >= λ @test abs(dim(domain(S5)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1) + + U0, S0, Vᴴ0 = svd_compact(t, alg) + Sd = similar(S0) + copyto!(Sd.data, real(T)[3 - mod(i, 3) for i in 0:(length(Sd.data) - 1)]) + td = U0 * Sd * Vᴴ0 + rtol = sqrt(eps(real(T))) + _, Sr, _, _ = svd_trunc(td; trunc = truncrank(nvals), alg) + for mode in (:down, :up) + trunc = truncmultiplets(truncrank(nvals); rtol, mode) + U6, S6, Vᴴ6, ϵ6 = @constinferred svd_trunc(td; trunc, alg) + @test td * Vᴴ6' ≈ U6 * S6 + @test isisometric(U6) + @test isisometric(Vᴴ6; side = :right) + @test norm(td - U6 * S6 * Vᴴ6) ≈ ϵ6 atol = eps(real(T))^(4 / 5) + if mode == :down + @test space(S6, 1) ≾ space(Sr, 1) + else + @test space(Sr, 1) ≾ space(S6, 1) + end + if dim(domain(S6)) > 0 + λ6 = minimum(diagview(S6)) + _, S7, _, _ = svd_trunc(td; trunc = trunctol(; atol = λ6 - 1 // 2), alg) + @test space(S6, 1) == space(S7, 1) + end + end end end diff --git a/test/cuda/factorizations.jl b/test/cuda/factorizations.jl index 3c482a04c..d560de296 100644 --- a/test/cuda/factorizations.jl +++ b/test/cuda/factorizations.jl @@ -453,6 +453,31 @@ for V in spacelist @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5) @test minimum(diagview(S5)) >= λ @test abs(dim(domain(S5)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1) + + U0, S0, Vᴴ0 = svd_compact(t) + Sd = similar(S0) + copyto!(Sd.data, real(T)[3 - mod(i, 3) for i in 0:(length(Sd.data) - 1)]) + td = U0 * Sd * Vᴴ0 + rtol = sqrt(eps(real(T))) + _, Sr, _, _ = svd_trunc(td; trunc = truncrank(nvals)) + for mode in (:down, :up) + trunc = truncmultiplets(truncrank(nvals); rtol, mode) + U6, S6, Vᴴ6, ϵ6 = @constinferred svd_trunc(td; trunc) + @test td * Vᴴ6' ≈ U6 * S6 + @test isisometric(U6) + @test isisometric(Vᴴ6; side = :right) + @test norm(td - U6 * S6 * Vᴴ6) ≈ ϵ6 atol = eps(real(T))^(4 / 5) + if mode == :down + @test space(S6, 1) ≾ space(Sr, 1) + else + @test space(Sr, 1) ≾ space(S6, 1) + end + if dim(domain(S6)) > 0 + λ6 = minimum(diagview(S6)) + _, S7, _, _ = svd_trunc(td; trunc = trunctol(; atol = λ6 - 1 // 2)) + @test space(S6, 1) == space(S7, 1) + end + end end end diff --git a/test/factorizations/svd.jl b/test/factorizations/svd.jl index e2eb19076..e6e8e9cf2 100644 --- a/test/factorizations/svd.jl +++ b/test/factorizations/svd.jl @@ -290,6 +290,31 @@ for V in spacelist @test norm(t - U5 * S5 * Vᴴ5) ≈ ϵ5 atol = eps(real(T))^(4 / 5) @test minimum(diagview(S5)) >= λ @test abs(dim(domain(S5)) - nvals) ≤ maximum(c -> blockdim(domain(t), c), blocksectors(t); init = 1) + + U0, S0, Vᴴ0 = svd_compact(t) + Sd = similar(S0) + copyto!(Sd.data, real(T)[3 - mod(i, 3) for i in 0:(length(Sd.data) - 1)]) + td = U0 * Sd * Vᴴ0 + rtol = sqrt(eps(real(T))) + _, Sr, _, _ = svd_trunc(td; trunc = truncrank(nvals)) + for mode in (:down, :up) + trunc = truncmultiplets(truncrank(nvals); rtol, mode) + U6, S6, Vᴴ6, ϵ6 = @constinferred svd_trunc(td; trunc) + @test td * Vᴴ6' ≈ U6 * S6 + @test isisometric(U6) + @test isisometric(Vᴴ6; side = :right) + @test norm(td - U6 * S6 * Vᴴ6) ≈ ϵ6 atol = eps(real(T))^(4 / 5) + if mode == :down + @test space(S6, 1) ≾ space(Sr, 1) + else + @test space(Sr, 1) ≾ space(S6, 1) + end + if dim(domain(S6)) > 0 + λ6 = minimum(diagview(S6)) + _, S7, _, _ = svd_trunc(td; trunc = trunctol(; atol = λ6 - 1 // 2)) + @test space(S6, 1) == space(S7, 1) + end + end end end end diff --git a/test/setup.jl b/test/setup.jl index ef8f1896f..dca3a6c4e 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -362,8 +362,8 @@ VIBMRepA4 = ( for V in (Vtr, VRepℤ₂, VRepℤ₃, VRepU₁, VfRepℤ₂, VTwistedVecℤ₄, VTwistedVecℤ₃, VRepCU₁, VRepSU₂, VfRepCU₁, VfRepSU₂, VfHubbard, Vfib, VRepA4, VfRepA4, VRepA4Twistedℤ₄, VIBM, VIBMRepA4) V1, V2, V3, V4, V5 = V for Vi in (V1, V2, V3, V4, V5) - @test allequal(leftunit, sectors(Vi)) - @test allequal(rightunit, sectors(Vi)) + @test TensorKit._allequal(leftunit, sectors(Vi)) + @test TensorKit._allequal(rightunit, sectors(Vi)) end @test dim(V1 ⊗ V2 ⊗ V3 ← (V4 ⊗ V5)') > 0 @test V1 ⊗ V2 ≾ (V3 ⊗ V4 ⊗ V5)' # wide blocks for QR