diff --git a/Project.toml b/Project.toml index 8a88983a..973c7045 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BlockSparseArrays" uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" authors = ["ITensor developers and contributors"] -version = "0.6.4" +version = "0.6.5" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl index bb3bf9e0..0c8e7716 100644 --- a/src/BlockSparseArrays.jl +++ b/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 diff --git a/src/factorizations/orthnull.jl b/src/factorizations/orthnull.jl new file mode 100644 index 00000000..b06af8e3 --- /dev/null +++ b/src/factorizations/orthnull.jl @@ -0,0 +1,124 @@ +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) + !isnothing(F) && throw( + ArgumentError( + "`left_orth!` on block sparse matrices does not support specifying the output" + ), + ) + return nothing +end + +function MatrixAlgebraKit.left_orth_qr!(A::AbstractBlockSparseMatrix, F, alg) + !isnothing(F) && throw( + ArgumentError( + "`left_orth!` on block sparse matrices does not support specifying the output" + ), + ) + alg′ = select_algorithm(qr_compact!, A, alg) + return qr_compact!(A, alg′) +end +function MatrixAlgebraKit.left_orth_polar!(A::AbstractBlockSparseMatrix, F, alg) + !isnothing(F) && throw( + ArgumentError( + "`left_orth!` on block sparse matrices does not support specifying the output" + ), + ) + 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 +) + !isnothing(F) && throw( + ArgumentError( + "`left_orth!` on block sparse matrices does not support specifying the output" + ), + ) + alg′ = select_algorithm(svd_compact!, A, alg) + U, S, Vᴴ = svd_compact!(A, alg′) + return U, S * Vᴴ +end +function MatrixAlgebraKit.left_orth_svd!(A::AbstractBlockSparseMatrix, F, alg, trunc) + !isnothing(F) && throw( + ArgumentError( + "`left_orth!` on block sparse matrices does not support specifying the output" + ), + ) + 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 + +function MatrixAlgebraKit.initialize_output( + ::typeof(right_orth!), A::AbstractBlockSparseMatrix +) + return nothing +end +function MatrixAlgebraKit.check_input( + ::typeof(right_orth!), A::AbstractBlockSparseMatrix, F::Nothing +) + !isnothing(F) && throw( + ArgumentError( + "`right_orth!` on block sparse matrices does not support specifying the output" + ), + ) + return nothing +end + +function MatrixAlgebraKit.right_orth_lq!(A::AbstractBlockSparseMatrix, F, alg) + !isnothing(F) && throw( + ArgumentError( + "`right_orth!` on block sparse matrices does not support specifying the output" + ), + ) + alg′ = select_algorithm(lq_compact!, A, alg) + return lq_compact!(A, alg′) +end +function MatrixAlgebraKit.right_orth_polar!(A::AbstractBlockSparseMatrix, F, alg) + !isnothing(F) && throw( + ArgumentError( + "`right_orth!` on block sparse matrices does not support specifying the output" + ), + ) + 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 +) + !isnothing(F) && throw( + ArgumentError( + "`right_orth!` on block sparse matrices does not support specifying the output" + ), + ) + 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) + !isnothing(F) && throw( + ArgumentError( + "`right_orth!` on block sparse matrices does not support specifying the output" + ), + ) + 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 diff --git a/src/factorizations/qr.jl b/src/factorizations/qr.jl index c28da925..8ec1ccdf 100644 --- a/src/factorizations/qr.jl +++ b/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 diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index 70639203..2ae4ebce 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -1,11 +1,13 @@ using BlockArrays: Block, BlockedMatrix, BlockedVector, blocks, mortar using BlockSparseArrays: BlockSparseArray, BlockDiagonal, eachblockstoredindex using MatrixAlgebraKit: + left_orth, left_polar, lq_compact, lq_full, qr_compact, qr_full, + right_orth, right_polar, svd_compact, svd_full, @@ -166,7 +168,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 +179,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 +239,37 @@ 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 + + U, C = left_orth(A; trunc=(; maxrank=2)) + @test size(U, 2) ≤ 2 + @test size(C, 1) ≤ 2 + @test Matrix(U'U) ≈ LinearAlgebra.I +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 + + C, U = right_orth(A; trunc=(; maxrank=2)) + @test size(C, 2) ≤ 2 + @test size(U, 1) ≤ 2 + @test Matrix(U * U') ≈ LinearAlgebra.I +end