Truncation composition

Project.toml

@@ -1,7 +1,7 @@
 name = "MatrixAlgebraKit"
 uuid = "6c742aac-3347-4629-af66-fc926824e5e4"
 authors = ["Jutho <[email protected]> and contributors"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/implementations/truncation.jl

@@ -11,10 +11,17 @@
     if isnothing(maxrank) && isnothing(atol) && isnothing(rtol)
         return NoTruncation()
     elseif isnothing(maxrank)
-        @assert isnothing(rtol) "TODO: rtol"
-        return trunctol(atol)
+        atol = @something atol 0
+        rtol = @something rtol 0
+        return TruncationKeepAbove(atol, rtol)
     else
-        return truncrank(maxrank)
+        if isnothing(atol) && isnothing(rtol)
+            return truncrank(maxrank)
+        else
+            atol = @something atol 0
+            rtol = @something rtol 0
+            return truncrank(maxrank) & TruncationKeepAbove(atol, rtol)
+        end
     end
 end
 
@@ -82,6 +89,27 @@
 """
 truncabove(atol) = TruncationKeepFiltered(≤(atol) ∘ abs)
 
+"""
+    TruncationComposition(trunc1::TruncationStrategy, trunc2::TruncationStrategy)
+Compose two truncation strategies, keeping values common between the two strategies.
+"""
+struct TruncationComposition{T<:Tuple{Vararg{TruncationStrategy}}} <:
+       TruncationStrategy
+    components::T
+end
+function Base.:&(trunc1::TruncationStrategy, trunc2::TruncationStrategy)
+    return TruncationComposition((trunc1, trunc2))
+end
+function Base.:&(trunc1::TruncationComposition, trunc2::TruncationComposition)
+    return TruncationComposition((trunc1.components..., trunc2.components...))
+end
+function Base.:&(trunc1::TruncationComposition, trunc2::TruncationStrategy)
+    return TruncationComposition((trunc1.components..., trunc2))
+end
+function Base.:&(trunc1::TruncationStrategy, trunc2::TruncationComposition)
+    return TruncationComposition((trunc1, trunc2.components...))
+end
+
 # truncate!
 # ---------
 # Generic implementation: `findtruncated` followed by indexing
@@ -147,6 +175,11 @@
     return 1:i
 end
 
+function findtruncated(values::AbstractVector, strategy::TruncationComposition)
+    inds = map(Base.Fix1(findtruncated, values), strategy.components)
+    return intersect(inds...)
+end
+
 """
     TruncatedAlgorithm(alg::AbstractAlgorithm, trunc::TruncationAlgorithm)
 

test/svd.jl

@@ -115,3 +115,30 @@ end
         end
     end
 end
+
+@testset "svd_trunc! mix maxrank and tol for T = $T" for T in
+                                                         (Float32, Float64, ComplexF32,
+                                                          ComplexF64)
+    rng = StableRNG(123)
+    if LinearAlgebra.LAPACK.version() < v"3.12.0"
+        algs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection())
+    else
+        algs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection(),
+                LAPACK_Jacobi())
+    end
+    m = 4
+    @testset "algorithm $alg" for alg in algs
+        U = qr_compact(randn(rng, T, m, m))[1]
+        S = Diagonal([0.9, 0.3, 0.1, 0.01])
+        Vᴴ = qr_compact(randn(rng, T, m, m))[1]
+        A = U * S * Vᴴ
+
+        U1, S1, V1ᴴ = svd_trunc(A; alg, trunc=(; rtol=0.2, maxrank=1))
+        @test length(S1.diag) == 1
+        @test S1.diag ≈ S.diag[1:1] rtol = sqrt(eps(real(T)))
+
+        U2, S2, V2ᴴ = svd_trunc(A; alg, trunc=(; rtol=0.2, maxrank=3))
+        @test length(S2.diag) == 2
+        @test S2.diag ≈ S.diag[1:2] rtol = sqrt(eps(real(T)))
+    end
+end