From b444a44414b41851433a603862647687b48e55dd Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 5 May 2025 17:17:49 -0400 Subject: [PATCH 01/12] Add initial truncation implementation --- src/BlockSparseArrays.jl | 1 + src/factorizations/truncation.jl | 62 ++++++++++++++++++++++++++++++++ test/test_factorizations.jl | 56 ++++++++++++++++++++++++++++- 3 files changed, 118 insertions(+), 1 deletion(-) create mode 100644 src/factorizations/truncation.jl diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl index 78ce8237..ffbfe937 100644 --- a/src/BlockSparseArrays.jl +++ b/src/BlockSparseArrays.jl @@ -44,5 +44,6 @@ include("BlockArraysSparseArraysBaseExt/BlockArraysSparseArraysBaseExt.jl") # factorizations include("factorizations/svd.jl") +include("factorizations/truncation.jl") end diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl new file mode 100644 index 00000000..37b12ca0 --- /dev/null +++ b/src/factorizations/truncation.jl @@ -0,0 +1,62 @@ +using MatrixAlgebraKit: TruncationStrategy, diagview + +const TBlockUSVᴴ = Tuple{ + <:AbstractBlockSparseMatrix,<:AbstractBlockSparseMatrix,<:AbstractBlockSparseMatrix +} + +function MatrixAlgebraKit.truncate!( + ::typeof(svd_trunc!), (U, S, Vᴴ)::TBlockUSVᴴ, strategy::TruncationStrategy +) + ind = MatrixAlgebraKit.findtruncated(diagview(S), strategy) + # cannot use regular slicing here: I want to slice without altering blockstructure + # solution: use boolean indexing and slice the mask, effectively cheaply inverting the map + indexmask = falses(size(S, 1)) + indexmask[ind] .= true + + # first determine the block structure of the output to avoid having assumptions on the + # data structures + ax = axes(S, 1) + counter = Base.Fix1(count, Base.Fix1(getindex, indexmask)) + Slengths = filter!(>(0), map(counter, blocks(ax))) + Sax = blockedrange(Slengths) + Ũ = similar(U, axes(U, 1), Sax) + S̃ = similar(S, Sax, Sax) + Ṽᴴ = similar(Vᴴ, Sax, axes(Vᴴ, 2)) + + # then loop over the blocks and assign the data + # TODO: figure out if we can presort and loop over the blocks - + # for now this has issues with missing blocks + bI_Us = collect(eachblockstoredindex(U)) + bI_Ss = collect(eachblockstoredindex(S)) + bI_Vᴴs = collect(eachblockstoredindex(Vᴴ)) + + I′ = 0 # number of skipped blocks that got fully truncated + for (I, b) in enumerate(blocks(ax)) + mask = indexmask[b] + + if !any(mask) + I′ += 1 + continue + end + + bU_id = @something findfirst(x -> last(Tuple(x)) == Block(I), bI_Us) error( + "No U-block found for $I" + ) + bU = Tuple(bI_Us[bU_id]) + Ũ[bU[1], bU[2] - Block(I′)] = view(U, bU...)[:, mask] + + bVᴴ_id = @something findfirst(x -> first(Tuple(x)) == Block(I), bI_Vᴴs) error( + "No Vᴴ-block found for $I" + ) + bVᴴ = Tuple(bI_Vᴴs[bVᴴ_id]) + Ṽᴴ[bVᴴ[1] - Block(I′), bVᴴ[2]] = view(Vᴴ, bVᴴ...)[mask, :] + + bS_id = @something findfirst(x -> last(Tuple(x)) == Block(I), bI_Ss) error( + "No S-block found for $I" + ) + bS = Tuple(bI_Ss[bS_id]) + S̃[(bS .- Block(I′))...] = Diagonal(diagview(view(S, bS...))[mask]) + end + + return Ũ, S̃, Ṽᴴ +end diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 7a28fdb7..394ab674 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,6 +1,6 @@ using BlockArrays: Block, BlockedMatrix, BlockedVector, blocks, mortar using BlockSparseArrays: BlockSparseArray, BlockDiagonal, eachblockstoredindex -using MatrixAlgebraKit: svd_compact, svd_full +using MatrixAlgebraKit: svd_compact, svd_full, svd_trunc, truncrank using LinearAlgebra: LinearAlgebra using Random: Random using Test: @inferred, @testset, @test @@ -83,3 +83,57 @@ end usv = svd_full(c) @test test_svd(c, usv; full=true) end + +# svd_trunc! +# ---------- + +@testset "svd_trunc ($m, $n) BlockSparseMatri{$T}" for ((m, n), T) in test_params + (m, n), T = first(test_params) + a = BlockSparseArray{T}(undef, m, n) + + # test blockdiagonal + for i in LinearAlgebra.diagind(blocks(a)) + I = CartesianIndices(blocks(a))[i] + a[Block(I.I...)] = rand(T, size(blocks(a)[i])) + end + + minmn = min(size(a)...) + r = max(1, minmn - 2) + + U1, S1, V1ᴴ = svd_trunc(a; trunc=truncrank(r)) + U2, S2, V2ᴴ = svd_trunc(Matrix(a); trunc=truncrank(r)) + @test size(U1) == size(U2) + @test size(S1) == size(S2) + @test size(V1ᴴ) == size(V2ᴴ) + @test Matrix(U1 * S1 * V1ᴴ) ≈ U2 * S2 * V2ᴴ + + @test (U1' * U1 ≈ LinearAlgebra.I) + @test (V1ᴴ * V1ᴴ' ≈ LinearAlgebra.I) + + # test permuted blockdiagonal + perm = Random.randperm(length(m)) + b = a[Block.(perm), Block.(1:length(n))] + U1, S1, V1ᴴ = svd_trunc(b; trunc=truncrank(r)) + U2, S2, V2ᴴ = svd_trunc(Matrix(b); trunc=truncrank(r)) + @test size(U1) == size(U2) + @test size(S1) == size(S2) + @test size(V1ᴴ) == size(V2ᴴ) + @test Matrix(U1 * S1 * V1ᴴ) ≈ U2 * S2 * V2ᴴ + + @test (U1' * U1 ≈ LinearAlgebra.I) + @test (V1ᴴ * V1ᴴ' ≈ LinearAlgebra.I) + + # test permuted blockdiagonal with missing row/col + I_removed = rand(eachblockstoredindex(b)) + c = copy(b) + delete!(blocks(c).storage, CartesianIndex(Int.(Tuple(I_removed)))) + U1, S1, V1ᴴ = svd_trunc(c; trunc=truncrank(r)) + U2, S2, V2ᴴ = svd_trunc(Matrix(c); trunc=truncrank(r)) + @test size(U1) == size(U2) + @test size(S1) == size(S2) + @test size(V1ᴴ) == size(V2ᴴ) + @test Matrix(U1 * S1 * V1ᴴ) ≈ U2 * S2 * V2ᴴ + + @test (U1' * U1 ≈ LinearAlgebra.I) + @test (V1ᴴ * V1ᴴ' ≈ LinearAlgebra.I) +end From b2f448c53caf7f66af53b330a8c850b65bba6355 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:00:53 -0400 Subject: [PATCH 02/12] Add dedicated truncation type --- src/factorizations/truncation.jl | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 37b12ca0..4d81676f 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -1,13 +1,27 @@ using MatrixAlgebraKit: TruncationStrategy, diagview +""" + BlockPermutedDiagonalTruncationStrategy(strategy::MatrixAlgebraKit.TruncationStrategy) + +A wrapper for `MatrixAlgebraKit.TruncationStrategy` that implements the wrapped strategy on +a block-by-block basis, which is possible if the input matrix is a block-diagonal matrix or +a block permuted block-diagonal matrix. +""" +struct BlockPermutedDiagonalTruncationStrategy{T<:MatrixAlgebraKit.TruncationStrategy} <: + MatrixAlgebraKit.TruncationStrategy + strategy::T +end + const TBlockUSVᴴ = Tuple{ <:AbstractBlockSparseMatrix,<:AbstractBlockSparseMatrix,<:AbstractBlockSparseMatrix } function MatrixAlgebraKit.truncate!( - ::typeof(svd_trunc!), (U, S, Vᴴ)::TBlockUSVᴴ, strategy::TruncationStrategy + ::typeof(svd_trunc!), + (U, S, Vᴴ)::TBlockUSVᴴ, + strategy::BlockPermutedDiagonalTruncationStrategy, ) - ind = MatrixAlgebraKit.findtruncated(diagview(S), strategy) + ind = MatrixAlgebraKit.findtruncated(diagview(S), strategy.strategy) # cannot use regular slicing here: I want to slice without altering blockstructure # solution: use boolean indexing and slice the mask, effectively cheaply inverting the map indexmask = falses(size(S, 1)) @@ -60,3 +74,4 @@ function MatrixAlgebraKit.truncate!( return Ũ, S̃, Ṽᴴ end + From 792fbd325d2ba59de525265d9553b80ba3061f81 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:01:01 -0400 Subject: [PATCH 03/12] format docstring --- src/factorizations/svd.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl index 4f64ff57..f83c11b3 100644 --- a/src/factorizations/svd.jl +++ b/src/factorizations/svd.jl @@ -4,7 +4,8 @@ using MatrixAlgebraKit: MatrixAlgebraKit, svd_compact!, svd_full! BlockPermutedDiagonalAlgorithm(A::MatrixAlgebraKit.AbstractAlgorithm) A wrapper for `MatrixAlgebraKit.AbstractAlgorithm` that implements the wrapped algorithm on -a block-by-block basis, which is possible if the input matrix is a block-diagonal matrix or a block permuted block-diagonal matrix. +a block-by-block basis, which is possible if the input matrix is a block-diagonal matrix or +a block permuted block-diagonal matrix. """ struct BlockPermutedDiagonalAlgorithm{A<:MatrixAlgebraKit.AbstractAlgorithm} <: MatrixAlgebraKit.AbstractAlgorithm From 92b706014693ab2762839ed83a93a79791d3e6b1 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:02:41 -0400 Subject: [PATCH 04/12] fix imports --- src/factorizations/truncation.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 4d81676f..cc55035f 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -1,14 +1,13 @@ -using MatrixAlgebraKit: TruncationStrategy, diagview +using MatrixAlgebraKit: TruncationStrategy, diagview, svd_trunc! """ - BlockPermutedDiagonalTruncationStrategy(strategy::MatrixAlgebraKit.TruncationStrategy) + BlockPermutedDiagonalTruncationStrategy(strategy::TruncationStrategy) -A wrapper for `MatrixAlgebraKit.TruncationStrategy` that implements the wrapped strategy on -a block-by-block basis, which is possible if the input matrix is a block-diagonal matrix or -a block permuted block-diagonal matrix. +A wrapper for `TruncationStrategy` that implements the wrapped strategy on a block-by-block +basis, which is possible if the input matrix is a block-diagonal matrix or a block permuted +block-diagonal matrix. """ -struct BlockPermutedDiagonalTruncationStrategy{T<:MatrixAlgebraKit.TruncationStrategy} <: - MatrixAlgebraKit.TruncationStrategy +struct BlockPermutedDiagonalTruncationStrategy{T<:TruncationStrategy} <: TruncationStrategy strategy::T end From b1ae30ffed132b8621763d8f935ec77b0dc65824 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:05:08 -0400 Subject: [PATCH 05/12] pass correct truncation strategy --- src/factorizations/truncation.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index cc55035f..edb88c77 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -15,6 +15,15 @@ const TBlockUSVᴴ = Tuple{ <:AbstractBlockSparseMatrix,<:AbstractBlockSparseMatrix,<:AbstractBlockSparseMatrix } +function MatrixAlgebraKit.truncate!( + ::typeof(svd_trunc!), (U, S, Vᴴ)::TBlockUSVᴴ, strategy::TruncationStrategy +) + # TODO assert blockdiagonal + return MatrixAlgebraKit.truncate!( + svd_trunc!, (U, S, Vᴴ), BlockPermutedDiagonalTruncationStrategy(strategy) + ) +end + function MatrixAlgebraKit.truncate!( ::typeof(svd_trunc!), (U, S, Vᴴ)::TBlockUSVᴴ, From 14596402469f6a26498a14c389be914783324c82 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 8 May 2025 18:07:17 +0000 Subject: [PATCH 06/12] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/factorizations/truncation.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index edb88c77..3f4e0572 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -82,4 +82,3 @@ function MatrixAlgebraKit.truncate!( return Ũ, S̃, Ṽᴴ end - From 2ea41a3577f37920b922e520906d460489172373 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:08:03 -0400 Subject: [PATCH 07/12] avoid `enumerate` --- src/factorizations/truncation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index edb88c77..99382178 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -53,7 +53,8 @@ function MatrixAlgebraKit.truncate!( bI_Vᴴs = collect(eachblockstoredindex(Vᴴ)) I′ = 0 # number of skipped blocks that got fully truncated - for (I, b) in enumerate(blocks(ax)) + for I in 1:blocksize(ax, 1) + b = ax[Block(I)] mask = indexmask[b] if !any(mask) @@ -82,4 +83,3 @@ function MatrixAlgebraKit.truncate!( return Ũ, S̃, Ṽᴴ end - From 3718571e04913cc6325d34bd5592efb24940089d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:10:07 -0400 Subject: [PATCH 08/12] Refactor truncation slightly --- src/factorizations/truncation.jl | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 99382178..44bcda71 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -24,16 +24,23 @@ function MatrixAlgebraKit.truncate!( ) end +# cannot use regular slicing here: I want to slice without altering blockstructure +# solution: use boolean indexing and slice the mask, effectively cheaply inverting the map +function MatrixAlgebraKit.findtruncated( + values::AbstractVector, strategy::BlockPermutedDiagonalTruncationStrategy +) + ind = MatrixAlgebraKit.findtruncated(values, strategy.strategy) + indexmask = falses(length(values)) + indexmask[ind] .= true + return indexmask +end + function MatrixAlgebraKit.truncate!( ::typeof(svd_trunc!), (U, S, Vᴴ)::TBlockUSVᴴ, strategy::BlockPermutedDiagonalTruncationStrategy, ) - ind = MatrixAlgebraKit.findtruncated(diagview(S), strategy.strategy) - # cannot use regular slicing here: I want to slice without altering blockstructure - # solution: use boolean indexing and slice the mask, effectively cheaply inverting the map - indexmask = falses(size(S, 1)) - indexmask[ind] .= true + indexmask = MatrixAlgebraKit.findtruncated(diagview(S), strategy) # first determine the block structure of the output to avoid having assumptions on the # data structures From dfa1c46a65b23761d12f9ea4ef5badda2362eb4f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:10:36 -0400 Subject: [PATCH 09/12] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index b13be4aa..9e61a2e4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BlockSparseArrays" uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" authors = ["ITensor developers and contributors"] -version = "0.5.0" +version = "0.5.1" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From 0c9f42d7e59b3156e45b966e8548e69806d079d8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:32:32 -0400 Subject: [PATCH 10/12] Create nicer objects for diagview --- src/factorizations/truncation.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index 44bcda71..ec7bad31 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -1,5 +1,15 @@ using MatrixAlgebraKit: TruncationStrategy, diagview, svd_trunc! +function MatrixAlgebraKit.diagview(A::BlockSparseMatrix{T,Diagonal{T,Vector{T}}}) where {T} + D = BlockSparseVector{T}(undef, axes(A, 1)) + for I in eachblockstoredindex(A) + if ==(Int.(Tuple(I))...) + D[Tuple(I)[1]] = diagview(A[I]) + end + end + return D +end + """ BlockPermutedDiagonalTruncationStrategy(strategy::TruncationStrategy) From 945c660a366670cd114b8f5b9a72171e1f1d5135 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:46:49 -0400 Subject: [PATCH 11/12] relax error criterion --- src/factorizations/truncation.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl index ec7bad31..0029c7a1 100644 --- a/src/factorizations/truncation.jl +++ b/src/factorizations/truncation.jl @@ -91,11 +91,11 @@ function MatrixAlgebraKit.truncate!( bVᴴ = Tuple(bI_Vᴴs[bVᴴ_id]) Ṽᴴ[bVᴴ[1] - Block(I′), bVᴴ[2]] = view(Vᴴ, bVᴴ...)[mask, :] - bS_id = @something findfirst(x -> last(Tuple(x)) == Block(I), bI_Ss) error( - "No S-block found for $I" - ) - bS = Tuple(bI_Ss[bS_id]) - S̃[(bS .- Block(I′))...] = Diagonal(diagview(view(S, bS...))[mask]) + bS_id = findfirst(x -> last(Tuple(x)) == Block(I), bI_Ss) + if !isnothing(bS_id) + bS = Tuple(bI_Ss[bS_id]) + S̃[(bS .- Block(I′))...] = Diagonal(diagview(view(S, bS...))[mask]) + end end return Ũ, S̃, Ṽᴴ From 17ccf8e5828f7c6334470f49ae8c239bc38dbe6b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 8 May 2025 14:47:03 -0400 Subject: [PATCH 12/12] Add tests --- test/test_factorizations.jl | 53 ++++++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 18 deletions(-) diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 394ab674..5152237b 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,6 +1,6 @@ using BlockArrays: Block, BlockedMatrix, BlockedVector, blocks, mortar using BlockSparseArrays: BlockSparseArray, BlockDiagonal, eachblockstoredindex -using MatrixAlgebraKit: svd_compact, svd_full, svd_trunc, truncrank +using MatrixAlgebraKit: svd_compact, svd_full, svd_trunc, truncrank, trunctol using LinearAlgebra: LinearAlgebra using Random: Random using Test: @inferred, @testset, @test @@ -88,7 +88,6 @@ end # ---------- @testset "svd_trunc ($m, $n) BlockSparseMatri{$T}" for ((m, n), T) in test_params - (m, n), T = first(test_params) a = BlockSparseArray{T}(undef, m, n) # test blockdiagonal @@ -99,9 +98,10 @@ end minmn = min(size(a)...) r = max(1, minmn - 2) + trunc = truncrank(r) - U1, S1, V1ᴴ = svd_trunc(a; trunc=truncrank(r)) - U2, S2, V2ᴴ = svd_trunc(Matrix(a); trunc=truncrank(r)) + U1, S1, V1ᴴ = svd_trunc(a; trunc) + U2, S2, V2ᴴ = svd_trunc(Matrix(a); trunc) @test size(U1) == size(U2) @test size(S1) == size(S2) @test size(V1ᴴ) == size(V2ᴴ) @@ -110,11 +110,11 @@ end @test (U1' * U1 ≈ LinearAlgebra.I) @test (V1ᴴ * V1ᴴ' ≈ LinearAlgebra.I) - # test permuted blockdiagonal - perm = Random.randperm(length(m)) - b = a[Block.(perm), Block.(1:length(n))] - U1, S1, V1ᴴ = svd_trunc(b; trunc=truncrank(r)) - U2, S2, V2ᴴ = svd_trunc(Matrix(b); trunc=truncrank(r)) + atol = minimum(LinearAlgebra.diag(S1)) + 10 * eps(real(T)) + trunc = trunctol(atol) + + U1, S1, V1ᴴ = svd_trunc(a; trunc) + U2, S2, V2ᴴ = svd_trunc(Matrix(a); trunc) @test size(U1) == size(U2) @test size(S1) == size(S2) @test size(V1ᴴ) == size(V2ᴴ) @@ -123,17 +123,34 @@ end @test (U1' * U1 ≈ LinearAlgebra.I) @test (V1ᴴ * V1ᴴ' ≈ LinearAlgebra.I) + # test permuted blockdiagonal + perm = Random.randperm(length(m)) + b = a[Block.(perm), Block.(1:length(n))] + for trunc in (truncrank(r), trunctol(atol)) + U1, S1, V1ᴴ = svd_trunc(b; trunc) + U2, S2, V2ᴴ = svd_trunc(Matrix(b); trunc) + @test size(U1) == size(U2) + @test size(S1) == size(S2) + @test size(V1ᴴ) == size(V2ᴴ) + @test Matrix(U1 * S1 * V1ᴴ) ≈ U2 * S2 * V2ᴴ + + @test (U1' * U1 ≈ LinearAlgebra.I) + @test (V1ᴴ * V1ᴴ' ≈ LinearAlgebra.I) + end + # test permuted blockdiagonal with missing row/col I_removed = rand(eachblockstoredindex(b)) c = copy(b) delete!(blocks(c).storage, CartesianIndex(Int.(Tuple(I_removed)))) - U1, S1, V1ᴴ = svd_trunc(c; trunc=truncrank(r)) - U2, S2, V2ᴴ = svd_trunc(Matrix(c); trunc=truncrank(r)) - @test size(U1) == size(U2) - @test size(S1) == size(S2) - @test size(V1ᴴ) == size(V2ᴴ) - @test Matrix(U1 * S1 * V1ᴴ) ≈ U2 * S2 * V2ᴴ - - @test (U1' * U1 ≈ LinearAlgebra.I) - @test (V1ᴴ * V1ᴴ' ≈ LinearAlgebra.I) + for trunc in (truncrank(r), trunctol(atol)) + U1, S1, V1ᴴ = svd_trunc(c; trunc) + U2, S2, V2ᴴ = svd_trunc(Matrix(c); trunc) + @test size(U1) == size(U2) + @test size(S1) == size(S2) + @test size(V1ᴴ) == size(V2ᴴ) + @test Matrix(U1 * S1 * V1ᴴ) ≈ U2 * S2 * V2ᴴ + + @test (U1' * U1 ≈ LinearAlgebra.I) + @test (V1ᴴ * V1ᴴ' ≈ LinearAlgebra.I) + end end