Implement left_orth/right_orth

Project.toml

@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.6.4"
+version = "0.6.5"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/BlockSparseArrays.jl

@@ -49,5 +49,6 @@ include("factorizations/truncation.jl")
 include("factorizations/qr.jl")
 include("factorizations/lq.jl")
 include("factorizations/polar.jl")
+include("factorizations/orthnull.jl")
 
 end

src/factorizations/orthnull.jl

@@ -0,0 +1,70 @@
+using MatrixAlgebraKit:
+  MatrixAlgebraKit,
+  left_orth!,
+  left_polar!,
+  lq_compact!,
+  qr_compact!,
+  right_orth!,
+  right_polar!,
+  select_algorithm,
+  svd_compact!
+
+function MatrixAlgebraKit.initialize_output(
+  ::typeof(left_orth!), A::AbstractBlockSparseMatrix
+)
+  return nothing
+end
+function MatrixAlgebraKit.check_input(
+  ::typeof(left_orth!), A::AbstractBlockSparseMatrix, F::Nothing
+)
+  return nothing
+end
+
+function MatrixAlgebraKit.left_orth_qr!(A::AbstractBlockSparseMatrix, F, alg)
+  alg′ = select_algorithm(qr_compact!, A, alg)
+  return qr_compact!(A, alg′)
+end
+function MatrixAlgebraKit.left_orth_polar!(A::AbstractBlockSparseMatrix, F, alg)
+  alg′ = select_algorithm(left_polar!, A, alg)
+  return left_polar!(A, alg′)
+end
+function MatrixAlgebraKit.left_orth_svd!(
+  A::AbstractBlockSparseMatrix, F, 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.initialize_output(
+  ::typeof(right_orth!), A::AbstractBlockSparseMatrix
+)
+  return nothing
+end
+function MatrixAlgebraKit.check_input(
+  ::typeof(right_orth!), A::AbstractBlockSparseMatrix, F::Nothing
+)
+  return nothing
+end
+
+function MatrixAlgebraKit.right_orth_lq!(A::AbstractBlockSparseMatrix, F, alg)
+  alg′ = select_algorithm(lq_compact!, A, alg)
+  return lq_compact!(A, alg′)
+end
+function MatrixAlgebraKit.right_orth_polar!(A::AbstractBlockSparseMatrix, F, alg)
+  alg′ = select_algorithm(right_polar!, A, alg)
+  return right_polar!(A, alg′)
+end
+function MatrixAlgebraKit.right_orth_svd!(
+  A::AbstractBlockSparseMatrix, F, 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_svd!(A::AbstractBlockSparseMatrix, F, 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

test/test_factorizations.jl

@@ -1,11 +1,18 @@
 using BlockArrays: Block, BlockedMatrix, BlockedVector, blocks, mortar
 using BlockSparseArrays: BlockSparseArray, BlockDiagonal, eachblockstoredindex
 using MatrixAlgebraKit:
+  left_orth,
+  left_polar,
+  qr_compact,
+  qr_full,
+  right_orth,
+  left_orth,
   left_polar,
   lq_compact,
   lq_full,
   qr_compact,
   qr_full,
+  right_orth,
   right_polar,
   svd_compact,
   svd_full,
@@ -166,7 +173,7 @@ end
   end
 end
 
-@testset "qr_compact" for T in (Float32, Float64, ComplexF32, ComplexF64)
+@testset "qr_compact (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]
     A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
     A[Block(1, 1)] = randn(T, i, k)
@@ -177,7 +184,7 @@ end
   end
 end
 
-@testset "qr_full" for T in (Float32, Float64, ComplexF32, ComplexF64)
+@testset "qr_full (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
   for i in [2, 3], j in [2, 3], k in [2, 3], l in [2, 3]
     A = BlockSparseArray{T}(undef, ([i, j], [k, l]))
     A[Block(1, 1)] = randn(T, i, k)
@@ -237,3 +244,27 @@ end
   @test C * U ≈ A
   @test Matrix(U * U') ≈ LinearAlgebra.I
 end
+
+@testset "left_orth (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+  A = BlockSparseArray{T}(undef, ([3, 4], [2, 3]))
+  A[Block(1, 1)] = randn(T, 3, 2)
+  A[Block(2, 2)] = randn(T, 4, 3)
+
+  for kind in (:polar, :qr, :svd)
+    U, C = left_orth(A; kind)
+    @test U * C ≈ A
+    @test Matrix(U'U) ≈ LinearAlgebra.I
+  end
+end
+
+@testset "right_orth (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+  A = BlockSparseArray{T}(undef, ([2, 3], [3, 4]))
+  A[Block(1, 1)] = randn(T, 2, 3)
+  A[Block(2, 2)] = randn(T, 3, 4)
+
+  for kind in (:lq, :polar, :svd)
+    C, U = right_orth(A; kind)
+    @test C * U ≈ A
+    @test Matrix(U * U') ≈ LinearAlgebra.I
+  end
+end