orth for non-AbstractMatrix

Author: mtfishman

src/implementations/orthnull.jl

@@ -107,13 +107,26 @@ function left_orth_polar!(A, VC, alg)
     return left_polar!(A, VC, alg′)
 end
 function left_orth_svd!(A, VC, alg, trunc::Nothing=nothing)
+    alg′ = select_algorithm(svd_compact!, A, alg)
+    U, S, Vᴴ = svd_compact!(A, alg′)
+    V, C = VC
+    return copy!(V, U), mul!(C, S, Vᴴ)
+end
+function left_orth_svd!(A::AbstractMatrix, VC, alg, trunc::Nothing=nothing)
     alg′ = select_algorithm(svd_compact!, A, alg)
     V, C = VC
     S = Diagonal(initialize_output(svd_vals!, A, alg′))
     U, S, Vᴴ = svd_compact!(A, (V, S, C), alg′)
     return U, lmul!(S, Vᴴ)
 end
 function left_orth_svd!(A, VC, alg, trunc)
+    alg′ = select_algorithm(svd_compact!, A, alg)
+    alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc)
+    U, S, Vᴴ = svd_trunc!(A, alg_trunc)
+    V, C = VC
+    return copy!(V, U), mul!(C, S, Vᴴ)
+end
+function left_orth_svd!(A::AbstractMatrix, VC, alg, trunc)
     alg′ = select_algorithm(svd_compact!, A, alg)
     alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc)
     V, C = VC
@@ -148,13 +161,26 @@ function right_orth_polar!(A, CVᴴ, alg)
     return right_polar!(A, CVᴴ, alg′)
 end
 function right_orth_svd!(A, CVᴴ, alg, trunc::Nothing=nothing)
+    alg′ = select_algorithm(svd_compact!, A, alg)
+    U, S, Vᴴ′ = svd_compact!(A, alg′)
+    C, Vᴴ = CVᴴ
+    return mul!(C, U, S), copy!(Vᴴ, Vᴴ′)
+end
+function right_orth_svd!(A::AbstractMatrix, CVᴴ, alg, trunc::Nothing=nothing)
     alg′ = select_algorithm(svd_compact!, A, alg)
     C, Vᴴ = CVᴴ
     S = Diagonal(initialize_output(svd_vals!, A, alg′))
     U, S, Vᴴ = svd_compact!(A, (C, S, Vᴴ), alg′)
     return rmul!(U, S), Vᴴ
 end
 function right_orth_svd!(A, CVᴴ, alg, trunc)
+    alg′ = select_algorithm(svd_compact!, A, alg)
+    alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc)
+    U, S, Vᴴ′ = svd_trunc!(A, alg_trunc)
+    C, Vᴴ = CVᴴ
+    return mul!(C, U, S), copy!(Vᴴ, Vᴴ′)
+end
+function right_orth_svd!(A::AbstractMatrix, CVᴴ, alg, trunc)
     alg′ = select_algorithm(svd_compact!, A, alg)
     alg_trunc = select_algorithm(svd_trunc!, A, alg′; trunc)
     C, Vᴴ = CVᴴ

src/interface/orthnull.jl

@@ -69,10 +69,10 @@ See also [`right_orth(!)`](@ref right_orth), [`left_null(!)`](@ref left_null), [
 """
 function left_orth end
 function left_orth! end
-function left_orth!(A::AbstractMatrix; kwargs...)
+function left_orth!(A; kwargs...)
     return left_orth!(A, initialize_output(left_orth!, A); kwargs...)
 end
-function left_orth(A::AbstractMatrix; kwargs...)
+function left_orth(A; kwargs...)
     return left_orth!(copy_input(left_orth, A); kwargs...)
 end
 
@@ -128,10 +128,10 @@ See also [`left_orth(!)`](@ref left_orth), [`left_null(!)`](@ref left_null), [`r
 """
 function right_orth end
 function right_orth! end
-function right_orth!(A::AbstractMatrix; kwargs...)
+function right_orth!(A; kwargs...)
     return right_orth!(A, initialize_output(right_orth!, A); kwargs...)
 end
-function right_orth(A::AbstractMatrix; kwargs...)
+function right_orth(A; kwargs...)
     return right_orth!(copy_input(right_orth, A); kwargs...)
 end
 
@@ -180,10 +180,10 @@ See also [`right_null(!)`](@ref right_null), [`left_orth(!)`](@ref left_orth), [
 """
 function left_null end
 function left_null! end
-function left_null!(A::AbstractMatrix; kwargs...)
+function left_null!(A; kwargs...)
     return left_null!(A, initialize_output(left_null!, A); kwargs...)
 end
-function left_null(A::AbstractMatrix; kwargs...)
+function left_null(A; kwargs...)
     return left_null!(copy_input(left_null, A); kwargs...)
 end
 
@@ -230,9 +230,9 @@ See also [`left_null(!)`](@ref left_null), [`left_orth(!)`](@ref left_orth), [`r
 """
 function right_null end
 function right_null! end
-function right_null!(A::AbstractMatrix; kwargs...)
+function right_null!(A; kwargs...)
     return right_null!(A, initialize_output(right_null!, A); kwargs...)
 end
-function right_null(A::AbstractMatrix; kwargs...)
+function right_null(A; kwargs...)
     return right_null!(copy_input(right_null, A); kwargs...)
 end

test/orthnull.jl

@@ -5,6 +5,41 @@ using StableRNGs
 using LinearAlgebra: LinearAlgebra, I
 using MatrixAlgebraKit: TruncationKeepAbove, TruncationKeepBelow
 
+# Used to test non-AbstractMatrix codepaths.
+struct LinearMap{P<:AbstractMatrix}
+    parent::P
+end
+using MatrixAlgebraKit: LAPACK_SVDAlgorithm, check_input, copy_input, default_svd_algorithm,
+                        initialize_output
+function MatrixAlgebraKit.copy_input(::typeof(qr_compact), A::LinearMap)
+    return LinearMap(copy_input(qr_compact, A.parent))
+end
+function MatrixAlgebraKit.copy_input(::typeof(lq_compact), A::LinearMap)
+    return LinearMap(copy_input(lq_compact, A.parent))
+end
+function MatrixAlgebraKit.initialize_output(::typeof(left_orth!), A::LinearMap)
+    return initialize_output(left_orth!, A.parent)
+end
+function MatrixAlgebraKit.initialize_output(::typeof(right_orth!), A::LinearMap)
+    return initialize_output(right_orth!, A.parent)
+end
+function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::LinearMap, VC)
+    return check_input(left_orth!, A.parent, VC)
+end
+function MatrixAlgebraKit.check_input(::typeof(right_orth!), A::LinearMap, VC)
+    return check_input(right_orth!, A.parent, VC)
+end
+function MatrixAlgebraKit.default_svd_algorithm(A::LinearMap)
+    return default_svd_algorithm(A.parent)
+end
+function MatrixAlgebraKit.initialize_output(::typeof(svd_compact!), A::LinearMap,
+                                            alg::LAPACK_SVDAlgorithm)
+    return initialize_output(svd_compact!, A.parent, alg)
+end
+function MatrixAlgebraKit.svd_compact!(A::LinearMap, USVᴴ, alg::LAPACK_SVDAlgorithm)
+    return svd_compact!(A.parent, USVᴴ, alg)
+end
+
 @testset "left_orth and left_null for T = $T" for T in (Float32, Float64, ComplexF32,
                                                         ComplexF64)
     rng = StableRNG(123)
@@ -23,6 +58,10 @@ using MatrixAlgebraKit: TruncationKeepAbove, TruncationKeepBelow
         @test N' * N ≈ I
         @test V * V' + N * N' ≈ I
 
+        M = LinearMap(A)
+        V, C = @constinferred left_orth(M; kind=:svd)
+        @test V * C ≈ A
+
         if m > n
             nullity = 5
             V, C = @constinferred left_orth(A)
@@ -162,6 +201,10 @@ end
         @test Nᴴ * Nᴴ' ≈ I
         @test Vᴴ' * Vᴴ + Nᴴ' * Nᴴ ≈ I
 
+        M = LinearMap(A)
+        C, Vᴴ = @constinferred right_orth(M; kind=:svd)
+        @test C * Vᴴ ≈ A
+
         Ac = similar(A)
         C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ))
         Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ)