Support for Diagonal in orthnull

src/interface/orthnull.jl

@@ -443,7 +443,9 @@ left_orth_alg(alg::LeftOrthAlgorithm) = alg
 left_orth_alg(alg::QRAlgorithms) = LeftOrthViaQR(alg)
 left_orth_alg(alg::PolarAlgorithms) = LeftOrthViaPolar(alg)
 left_orth_alg(alg::SVDAlgorithms) = LeftOrthViaSVD(alg)
+left_orth_alg(alg::DiagonalAlgorithm) = LeftOrthViaQR(alg)
 left_orth_alg(alg::TruncatedAlgorithm{<:SVDAlgorithms}) = LeftOrthViaSVD(alg)
+left_orth_alg(alg::TruncatedAlgorithm{DiagonalAlgorithm}) = LeftOrthViaSVD(alg)
 
 """
     right_orth_alg(alg::AbstractAlgorithm) -> RightOrthAlgorithm
@@ -478,7 +480,9 @@ right_orth_alg(alg::RightOrthAlgorithm) = alg
 right_orth_alg(alg::LQAlgorithms) = RightOrthViaLQ(alg)
 right_orth_alg(alg::PolarAlgorithms) = RightOrthViaPolar(alg)
 right_orth_alg(alg::SVDAlgorithms) = RightOrthViaSVD(alg)
+right_orth_alg(alg::DiagonalAlgorithm) = RightOrthViaLQ(alg)
 right_orth_alg(alg::TruncatedAlgorithm{<:SVDAlgorithms}) = RightOrthViaSVD(alg)
+right_orth_alg(alg::TruncatedAlgorithm{DiagonalAlgorithm}) = RightOrthViaSVD(alg)
 
 """
     left_null_alg(alg::AbstractAlgorithm) -> LeftNullAlgorithm

test/orthnull.jl

@@ -19,16 +19,18 @@ for T in (BLASFloats..., GenericFloats...), n in (37, m, 63)
     if T ∈ BLASFloats
         if CUDA.functional()
             TestSuite.test_orthnull(CuMatrix{T}, (m, n); test_nullity = false)
-            n == m && TestSuite.test_orthnull(Diagonal{T, CuVector{T}}, m; test_orthnull = false)
+            n == m && TestSuite.test_orthnull(Diagonal{T, CuVector{T}}, m)
         end
         if AMDGPU.functional()
             TestSuite.test_orthnull(ROCMatrix{T}, (m, n); test_nullity = false)
-            n == m && TestSuite.test_orthnull(Diagonal{T, ROCVector{T}}, m; test_orthnull = false)
+            n == m && TestSuite.test_orthnull(Diagonal{T, ROCVector{T}}, m)
+        end
+        if AMDGPU.functional()
         end
     end
     if !is_buildkite
         TestSuite.test_orthnull(T, (m, n))
         AT = Diagonal{T, Vector{T}}
-        TestSuite.test_orthnull(AT, m; test_orthnull = false)
+        TestSuite.test_orthnull(AT, m)
     end
 end

test/testsuite/TestSuite.jl

@@ -73,9 +73,11 @@ is_pivoted(alg::MatrixAlgebraKit.LQViaTransposedQR) = is_pivoted(alg.qr_alg)
 isleftcomplete(V, N) = V * V' + N * N' ≈ I
 isleftcomplete(V::AnyCuMatrix, N::AnyCuMatrix) = isleftcomplete(collect(V), collect(N))
 isleftcomplete(V::AnyROCMatrix, N::AnyROCMatrix) = isleftcomplete(collect(V), collect(N))
+isleftcomplete(V::Diagonal{TV, <:AnyROCVector}, N::Diagonal{TN, <:AnyROCVector}) where {TV, TN} = isleftcomplete(Diagonal(collect(V.diag)), Diagonal(collect(N.diag)))
 isrightcomplete(Vᴴ, Nᴴ) = Vᴴ' * Vᴴ + Nᴴ' * Nᴴ ≈ I
 isrightcomplete(V::AnyCuMatrix, N::AnyCuMatrix) = isrightcomplete(collect(V), collect(N))
 isrightcomplete(V::AnyROCMatrix, N::AnyROCMatrix) = isrightcomplete(collect(V), collect(N))
+isrightcomplete(V::Diagonal{TV, <:AnyROCVector}, N::Diagonal{TN, <:AnyROCVector}) where {TV, TN} = isrightcomplete(Diagonal(collect(V.diag)), Diagonal(collect(N.diag)))
 
 instantiate_unitary(T, A, sz) = qr_compact(randn!(similar(A, eltype(T), sz, sz)))[1]
 # AMDGPU can't generate ComplexF32 random numbers

test/testsuite/orthnull.jl

@@ -0,0 +1,324 @@
+using TestExtras
+using LinearAlgebra
+
+include("../linearmap.jl")
+
+_left_orth_svd(x; kwargs...) = left_orth(x; alg = :svd, kwargs...)
+_left_orth_svd!(x, VC; kwargs...) = left_orth!(x, VC; alg = :svd, kwargs...)
+_left_orth_qr(x; kwargs...) = left_orth(x; alg = :qr, kwargs...)
+_left_orth_qr!(x, VC; kwargs...) = left_orth!(x, VC; alg = :qr, kwargs...)
+_left_orth_polar(x; kwargs...) = left_orth(x; alg = :polar, kwargs...)
+_left_orth_polar!(x, VC; kwargs...) = left_orth!(x, VC; alg = :polar, kwargs...)
+
+_right_orth_svd(x; kwargs...) = right_orth(x; alg = :svd, kwargs...)
+_right_orth_svd!(x, CVᴴ; kwargs...) = right_orth!(x, CVᴴ; alg = :svd, kwargs...)
+_right_orth_lq(x; kwargs...) = right_orth(x; alg = :lq, kwargs...)
+_right_orth_lq!(x, CVᴴ; kwargs...) = right_orth!(x, CVᴴ; alg = :lq, kwargs...)
+_right_orth_polar(x; kwargs...) = right_orth(x; alg = :polar, kwargs...)
+_right_orth_polar!(x, CVᴴ; kwargs...) = right_orth!(x, CVᴴ; alg = :polar, kwargs...)
+
+function test_orthnull(T::Type, sz; test_nullity = true, kwargs...)
+    summary_str = testargs_summary(T, sz)
+    return @testset "orthnull $summary_str" begin
+        test_left_orthnull(T, sz; kwargs...)
+        test_nullity  && test_left_nullity(T, sz; kwargs...)
+        test_right_orthnull(T, sz; kwargs...)
+        test_nullity  && test_right_nullity(T, sz; kwargs...)
+    end
+end
+
+function test_left_nullity(
+        T::Type, sz;
+        atol::Real = 0, rtol::Real = precision(T),
+        kwargs...
+    )
+    summary_str = testargs_summary(T, sz)
+    return @testset "left_nullity $summary_str" begin
+        A = instantiate_matrix(T, sz)
+        Ac = deepcopy(A)
+        m, n = size(A)
+        minmn = min(m, n)
+        if m > n
+            nullity = 5
+            V, C = @testinferred left_orth(A)
+            N = @testinferred left_null(A; trunc = (; maxnullity = nullity))
+            @test V isa typeof(A) && size(V) == (m, minmn)
+            @test C isa typeof(A) && size(C) == (minmn, n)
+            @test eltype(N) == eltype(A) && size(N) == (m, nullity)
+            @test V * C ≈ A
+            @test isisometric(V)
+            @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+            @test isisometric(N)
+        end
+
+        rtol = eps(real(eltype(T)))
+        for (trunc_orth, trunc_null) in (
+                ((; rtol = rtol), (; rtol = rtol)),
+                (trunctol(; rtol), trunctol(; rtol, keep_below = true)),
+            )
+            V, C = left_orth(A)
+            N = left_null(A)
+            V2, C2 = @testinferred left_orth!(copy!(Ac, A), (V, C); trunc = trunc_orth)
+            N2 = @testinferred left_null!(copy!(Ac, A), N; trunc = trunc_null)
+            @test V2 * C2 ≈ A
+            @test isisometric(V2)
+            @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+            @test isisometric(N2)
+            @test isleftcomplete(V2, N2)
+        end
+
+        alg = :svd
+        V2, C2 = @testinferred _left_orth_svd!(copy!(Ac, A), (V, C); trunc = (; atol))
+        N2 = @testinferred left_null!(copy!(Ac, A), N; alg, trunc = (; atol))
+        @test V2 * C2 ≈ A
+        @test isisometric(V2)
+        @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+        @test isisometric(N2)
+        @test isleftcomplete(V2, N2)
+
+        V2, C2 = @testinferred _left_orth_svd!(copy!(Ac, A), (V, C); trunc = (; rtol))
+        N2 = @testinferred left_null!(copy!(Ac, A), N; alg, trunc = (; rtol))
+        @test V2 * C2 ≈ A
+        @test isisometric(V2)
+        @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+        @test isisometric(N2)
+        @test isleftcomplete(V2, N2)
+
+        # doesn't work on AMD...
+        atol = eps(real(eltype(T)))
+        V2, C2 = @testinferred left_orth!(copy!(Ac, A), (V, C); trunc = (; atol = atol))
+        N2 = @testinferred left_null!(copy!(Ac, A), N; trunc = (; atol = atol))
+        @test V2 * C2 ≈ A
+        @test isisometric(V2)
+        @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
+        @test isisometric(N2)
+        @test isleftcomplete(V2, N2)
+
+    end
+end
+
+function test_left_orthnull(
+        T::Type, sz;
+        atol::Real = 0, rtol::Real = precision(T),
+        kwargs...
+    )
+    summary_str = testargs_summary(T, sz)
+    return @testset "left_orth! and left_null! $summary_str" begin
+        A = instantiate_matrix(T, sz)
+        Ac = deepcopy(A)
+        V, C = @testinferred left_orth(A)
+        N = @testinferred left_null(A)
+        m, n = size(A)
+        minmn = min(m, n)
+        @test V isa typeof(A) && size(V) == (m, minmn)
+        @test C isa typeof(A) && size(C) == (minmn, n)
+        @test eltype(N) == eltype(A) && size(N) == (m, m - minmn)
+        @test V * C ≈ A
+        @test isisometric(V)
+        @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
+        @test isisometric(N)
+        @test isleftcomplete(V, N)
+
+        M = LinearMap(A)
+        VM, CM = @testinferred _left_orth_svd(M)
+        @test parent(VM) * parent(CM) ≈ A
+
+        # passing a kind and some kwargs
+        V, C = @testinferred _left_orth_qr(A; positive = true)
+        N = @testinferred left_null(A; alg = :qr, positive = true)
+        @test V isa typeof(A) && size(V) == (m, minmn)
+        @test C isa typeof(A) && size(C) == (minmn, n)
+        @test eltype(N) == eltype(A) && size(N) == (m, m - minmn)
+        @test V * C ≈ A
+        @test isisometric(V)
+        @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
+        @test isisometric(N)
+        @test isleftcomplete(V, N)
+
+        # passing an algorithm
+        V, C = @testinferred left_orth(A; alg = MatrixAlgebraKit.default_qr_algorithm(A))
+        N = @testinferred left_null(A; alg = :qr, positive = true)
+        @test V isa typeof(A) && size(V) == (m, minmn)
+        @test C isa typeof(A) && size(C) == (minmn, n)
+        @test eltype(N) == eltype(A) && size(N) == (m, m - minmn)
+        @test V * C ≈ A
+        @test isisometric(V)
+        @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
+        @test isisometric(N)
+        @test isleftcomplete(V, N)
+
+        V2, C2 = @testinferred left_orth!(copy!(Ac, A), (V, C))
+        N2 = @testinferred left_null!(copy!(Ac, A), N)
+        @test V2 * C2 ≈ A
+        @test isisometric(V2)
+        @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
+        @test isisometric(N2)
+        @test isleftcomplete(V2, N2)
+
+        for alg in (:qr, :polar, :svd) # explicit kind kwarg
+            m < n && alg === :polar && continue
+            if alg == :svd
+                V2, C2 = @testinferred _left_orth_svd!(copy!(Ac, A), (V, C))
+            elseif alg == :qr
+                V2, C2 = @testinferred _left_orth_qr!(copy!(Ac, A), (V, C))
+            elseif alg == :polar
+                V2, C2 = @testinferred _left_orth_polar!(copy!(Ac, A), (V, C))
+            end
+            @test V2 * C2 ≈ A
+            @test isisometric(V2)
+            if alg != :polar
+                N2 = @testinferred left_null!(copy!(Ac, A), N; alg)
+                @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
+                @test isisometric(N2)
+                @test isleftcomplete(V2, N2)
+            end
+
+            # with kind and tol kwargs
+            if alg != :svd
+                @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; atol))
+                @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; rtol))
+                alg == :polar && continue
+                @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; atol))
+                @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; rtol))
+            end
+        end
+    end
+end
+
+function test_right_nullity(
+        T::Type, sz;
+        atol::Real = 0, rtol::Real = precision(T),
+        kwargs...
+    )
+    summary_str = testargs_summary(T, sz)
+    return @testset "right_nullity $summary_str" begin
+        A = instantiate_matrix(T, sz)
+        Ac = deepcopy(A)
+        m, n = size(A)
+        minmn = min(m, n)
+
+        C, Vᴴ = @testinferred right_orth(A)
+        Nᴴ = @testinferred right_null(A)
+        atol = eps(real(eltype(T)))
+        C2, Vᴴ2 = @testinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; atol))
+        Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; atol))
+        @test C2 * Vᴴ2 ≈ A
+        @test isisometric(Vᴴ2; side = :right)
+        @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+        @test isisometric(Nᴴ; side = :right)
+        @test isrightcomplete(Vᴴ2, Nᴴ2)
+
+        rtol = eps(real(eltype(T)))
+        C2, Vᴴ2 = @testinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; rtol))
+        Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; rtol))
+        @test C2 * Vᴴ2 ≈ A
+        @test isisometric(Vᴴ2; side = :right)
+        @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+        @test isisometric(Nᴴ2; side = :right)
+        @test isrightcomplete(Vᴴ2, Nᴴ2)
+
+        alg = :svd
+        C2, Vᴴ2 = @testinferred _right_orth_svd!(copy!(Ac, A), (C, Vᴴ); trunc = (; atol))
+        Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; alg = alg, trunc = (; atol))
+        @test C2 * Vᴴ2 ≈ A
+        @test isisometric(Vᴴ2; side = :right)
+        @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
+        @test isisometric(Nᴴ2; side = :right)
+        @test isrightcomplete(Vᴴ2, Nᴴ2)
+
+        C2, Vᴴ2 = @testinferred _right_orth_svd!(copy!(Ac, A), (C, Vᴴ); trunc = (; rtol))
+        Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; alg = alg, trunc = (; rtol))
+        @test C2 * Vᴴ2 ≈ A
+        @test isisometric(Vᴴ2; side = :right)
+        @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
+        @test isisometric(Nᴴ2; side = :right)
+        @test isrightcomplete(Vᴴ2, Nᴴ2)
+    end
+end
+
+function test_right_orthnull(
+        T::Type, sz;
+        atol::Real = 0, rtol::Real = precision(T),
+        kwargs...
+    )
+    summary_str = testargs_summary(T, sz)
+    return @testset "right_orth! and right_null! $summary_str" begin
+        A = instantiate_matrix(T, sz)
+        m, n = size(A)
+        minmn = min(m, n)
+        Ac = deepcopy(A)
+        C, Vᴴ = @testinferred right_orth(A)
+        Nᴴ = @testinferred right_null(A)
+        @test C isa typeof(A) && size(C) == (m, minmn)
+        @test Vᴴ isa typeof(A) && size(Vᴴ) == (minmn, n)
+        @test eltype(Nᴴ) == eltype(A) && size(Nᴴ) == (n - minmn, n)
+        @test C * Vᴴ ≈ A
+        @test isisometric(Vᴴ; side = :right)
+        @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+        @test isisometric(Nᴴ; side = :right)
+        @test isrightcomplete(Vᴴ, Nᴴ)
+
+        M = LinearMap(A)
+        CM, VMᴴ = @testinferred _right_orth_svd(M)
+        @test parent(CM) * parent(VMᴴ) ≈ A
+
+        # passing a kind and some kwargs
+        C, Vᴴ = @testinferred _right_orth_lq(A; positive = true)
+        Nᴴ = @testinferred right_null(A; alg = :lq, positive = true)
+        @test C isa typeof(A) && size(C) == (m, minmn)
+        @test Vᴴ isa typeof(A) && size(Vᴴ) == (minmn, n)
+        @test eltype(Nᴴ) == eltype(A) && size(Nᴴ) == (n - minmn, n)
+        @test C * Vᴴ ≈ A
+        @test isisometric(Vᴴ; side = :right)
+        @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+        @test isisometric(Nᴴ; side = :right)
+        @test isrightcomplete(Vᴴ, Nᴴ)
+
+        # passing an algorithm
+        C, Vᴴ = @testinferred right_orth(A; alg = MatrixAlgebraKit.default_lq_algorithm(A))
+        @test C isa typeof(A) && size(C) == (m, minmn)
+        @test Vᴴ isa typeof(A) && size(Vᴴ) == (minmn, n)
+        Nᴴ = @testinferred right_null(A; alg = :lq, positive = true)
+        @test eltype(Nᴴ) == eltype(A) && size(Nᴴ) == (n - minmn, n)
+        @test C * Vᴴ ≈ A
+        @test isisometric(Vᴴ; side = :right)
+        @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+        @test isisometric(Nᴴ; side = :right)
+        @test isrightcomplete(Vᴴ, Nᴴ)
+
+        C2, Vᴴ2 = @testinferred right_orth!(copy!(Ac, A), (C, Vᴴ))
+        Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ)
+        @test C2 * Vᴴ2 ≈ A
+        @test isisometric(Vᴴ2; side = :right)
+        @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+        @test isisometric(Nᴴ; side = :right)
+        @test isrightcomplete(Vᴴ2, Nᴴ2)
+
+        for alg in (:lq, :polar, :svd)
+            n < m && alg == :polar && continue
+            if alg == :lq
+                C2, Vᴴ2 = @testinferred _right_orth_lq!(copy!(Ac, A), (C, Vᴴ))
+            elseif alg == :polar
+                C2, Vᴴ2 = @testinferred _right_orth_polar!(copy!(Ac, A), (C, Vᴴ))
+            elseif alg == :svd
+                C2, Vᴴ2 = @testinferred _right_orth_svd!(copy!(Ac, A), (C, Vᴴ))
+            end
+            @test C2 * Vᴴ2 ≈ A
+            @test isisometric(Vᴴ2; side = :right)
+            if alg != :polar
+                Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; alg = alg)
+                @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T))
+                @test isisometric(Nᴴ2; side = :right)
+                @test isrightcomplete(Vᴴ2, Nᴴ2)
+            end
+
+            if alg != :svd
+                @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; atol))
+                @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; rtol))
+                alg == :polar && continue
+                @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; atol))
+                @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; rtol))
+            end
+        end
+    end
+end