Merge 835c523 into 13d0332

Author: mtfishman

Project.toml

@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 authors = ["ITensor developers <support@itensor.org> 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,96 @@
+using MatrixAlgebraKit:
+  MatrixAlgebraKit,
+  left_orth_polar!,
+  left_orth_qr!,
+  left_orth_svd!,
+  left_polar!,
+  lq_compact!,
+  qr_compact!,
+  right_orth_lq!,
+  right_orth_polar!,
+  right_orth_svd!,
+  right_polar!,
+  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
+    return left_orth_polar!(A, alg_polar)
+  elseif kind == :svd
+    return left_orth_svd!(A, alg_svd, trunc)
+  else
+    throw(ArgumentError("`left_orth` received unknown value `kind = $kind`"))
+  end
+end
+function MatrixAlgebraKit.left_orth_qr!(A::AbstractBlockSparseMatrix, alg)
+  alg′ = select_algorithm(qr_compact!, A, alg)
+  return qr_compact!(A, alg′)
+end
+function MatrixAlgebraKit.left_orth_polar!(A::AbstractBlockSparseMatrix, alg)
+  alg′ = select_algorithm(left_polar!, A, alg)
+  return left_polar!(A, alg′)
+end
+function MatrixAlgebraKit.left_orth_svd!(
+  A::AbstractBlockSparseMatrix, 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::AbstractBlockSparseMatrix;
+  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
+    return right_orth_polar!(A, alg_polar)
+  elseif kind == :svd
+    return right_orth_svd!(A, alg_svd, trunc)
+  else
+    throw(ArgumentError("`right_orth` received unknown value `kind = $kind`"))
+  end
+end
+function MatrixAlgebraKit.right_orth_lq!(A::AbstractBlockSparseMatrix, alg)
+  alg′ = select_algorithm(lq_compact, A, alg)
+  return lq_compact!(A, alg′)
+end
+function MatrixAlgebraKit.right_orth_polar!(A::AbstractBlockSparseMatrix, alg)
+  alg′ = select_algorithm(right_polar!, A, alg)
+  return right_polar!(A, alg′)
+end
+function MatrixAlgebraKit.right_orth_svd!(
+  A::AbstractBlockSparseMatrix, 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, 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
 )

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 (:qr, :polar, :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 (:qr, :polar, :svd)
+    C, U = right_orth(A; kind)
+    @test C * U ≈ A
+    @test Matrix(U * U') ≈ LinearAlgebra.I
+  end
+end