Better implementation of left_orth/right_orth

Author: mtfishman

src/BlockSparseArrays.jl

@@ -47,5 +47,6 @@ include("BlockArraysSparseArraysBaseExt/BlockArraysSparseArraysBaseExt.jl")
 include("factorizations/svd.jl")
 include("factorizations/truncation.jl")
 include("factorizations/qr.jl")
+include("factorizations/orthnull.jl")
 
 end

src/factorizations/orthnull.jl

@@ -0,0 +1,84 @@
+using MatrixAlgebraKit:
+  MatrixAlgebraKit, left_polar, qr_compact, select_algorithm, svd_compact
+
+function MatrixAlgebraKit.left_orth(
+  A::AbstractBlockSparseMatrix;
+  trunc=nothing,
+  kind=isnothing(trunc) ? :qr : :svd,
+  alg_qr=(; positive=true),
+  alg_polar=(;),
+  alg_svd=(;),
+)
+  if !isnothing(trunc) && kind != :svd
+    throw(ArgumentError("truncation not supported for `left_orth` with `kind=$kind`"))
+  end
+  if kind == :qr
+    return left_orth_qr(A, alg_qr)
+  elseif kind == :polar
+    # TODO: Implement this.
+    # return left_orth_polar(A, alg_polar)
+    return left_orth_svd(A, alg_svd)
+  elseif kind == :svd
+    return left_orth_svd(A, alg_svd, trunc)
+  else
+    throw(ArgumentError("`left_orth` received unknown value `kind = $kind`"))
+  end
+end
+function left_orth_qr(A, alg)
+  alg′ = select_algorithm(qr_compact, A, alg)
+  return qr_compact(A, alg′)
+end
+function left_orth_polar(A, alg)
+  alg′ = select_algorithm(left_polar, A, alg)
+  return left_polar(A, alg′)
+end
+function left_orth_svd(A, alg, trunc::Nothing=nothing)
+  alg′ = select_algorithm(svd_compact, A, alg)
+  U, S, Vᴴ = svd_compact(A, alg′)
+  return U, S * Vᴴ
+end
+
+function MatrixAlgebraKit.right_orth(
+  A;
+  trunc=nothing,
+  kind=isnothing(trunc) ? :lq : :svd,
+  alg_lq=(; positive=true),
+  alg_polar=(;),
+  alg_svd=(;),
+)
+  if !isnothing(trunc) && kind != :svd
+    throw(ArgumentError("truncation not supported for `right_orth` with `kind=$kind`"))
+  end
+  if kind == :qr
+    # TODO: Implement this.
+    # return right_orth_lq(A, alg_lq)
+    return right_orth_svd(A, alg_svd)
+  elseif kind == :polar
+    # TODO: Implement this.
+    # return right_orth_polar(A, alg_polar)
+    return right_orth_svd(A, alg_svd)
+  elseif kind == :svd
+    return right_orth_svd(A, alg_svd, trunc)
+  else
+    throw(ArgumentError("`right_orth` received unknown value `kind = $kind`"))
+  end
+end
+function right_orth_lq(A, alg)
+  alg′ = select_algorithm(lq_compact, A, alg)
+  return lq_compact(A, alg′)
+end
+function right_orth_polar(A, alg)
+  alg′ = select_algorithm(right_polar, A, alg)
+  return right_polar(A, alg′)
+end
+function right_orth_svd(A, alg, trunc::Nothing=nothing)
+  alg′ = select_algorithm(svd_compact, A, alg)
+  U, S, Vᴴ = svd_compact(A, alg′)
+  return U * S, Vᴴ
+end
+function right_orth_svd(A, 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)
+  return U * S, Vᴴ
+end

src/factorizations/qr.jl

@@ -1,4 +1,4 @@
-using MatrixAlgebraKit: MatrixAlgebraKit, qr_compact!, qr_full!
+using MatrixAlgebraKit: MatrixAlgebraKit, lq_compact!, lq_full!, qr_compact!, qr_full!
 
 # TODO: this is a hardcoded for now to get around this function not being defined in the
 # type domain
@@ -19,6 +19,23 @@ function MatrixAlgebraKit.default_algorithm(
   return default_blocksparse_qr_algorithm(A; kwargs...)
 end
 
+function default_blocksparse_lq_algorithm(A::AbstractMatrix; kwargs...)
+  blocktype(A) <: StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat} ||
+    error("unsupported type: $(blocktype(A))")
+  alg = MatrixAlgebraKit.LAPACK_HouseholderLQ(; kwargs...)
+  return BlockPermutedDiagonalAlgorithm(alg)
+end
+function MatrixAlgebraKit.default_algorithm(
+  ::typeof(lq_compact!), A::AbstractBlockSparseMatrix; kwargs...
+)
+  return default_blocksparse_lq_algorithm(A; kwargs...)
+end
+function MatrixAlgebraKit.default_algorithm(
+  ::typeof(lq_full!), A::AbstractBlockSparseMatrix; kwargs...
+)
+  return default_blocksparse_lq_algorithm(A; kwargs...)
+end
+
 function similar_output(
   ::typeof(qr_compact!), A, R_axis, alg::MatrixAlgebraKit.AbstractAlgorithm
 )

src/factorizations/svd.jl

@@ -265,7 +265,3 @@ function MatrixAlgebraKit.svd_full!(
 
   return USVᴴ
 end
-
-function MatrixAlgebraKit.left_orth_svd!(A::AbstractBlockSparseMatrix, VC, alg, trunc)
-  return @invoke MatrixAlgebraKit.left_orth_svd!(A::Any, VC::Any, alg::Any, trunc::Any)
-end