diff --git a/Project.toml b/Project.toml index bff14a8..6d636b7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BlockSparseArrays" uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" +version = "0.10.11" authors = ["ITensor developers and contributors"] -version = "0.10.10" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" @@ -41,7 +41,7 @@ GPUArraysCore = "0.1, 0.2" LinearAlgebra = "1.10" MacroTools = "0.5.13" MapBroadcast = "0.1.5" -MatrixAlgebraKit = "0.5" +MatrixAlgebraKit = "0.6" SparseArraysBase = "0.7.1" SplitApplyCombine = "1.2.3" TensorAlgebra = "0.3, 0.4" diff --git a/src/factorizations/orthnull.jl b/src/factorizations/orthnull.jl index f1adf27..c023f8f 100644 --- a/src/factorizations/orthnull.jl +++ b/src/factorizations/orthnull.jl @@ -1,26 +1,21 @@ -using MatrixAlgebraKit: - MatrixAlgebraKit, - default_svd_algorithm, - left_null!, - left_null_svd!, - left_orth!, - left_polar!, - lq_compact!, - null_truncation_strategy, - qr_compact!, - right_null!, - right_null_svd!, - right_orth!, - right_polar!, - select_algorithm, - svd_compact! +using MatrixAlgebraKit: MatrixAlgebraKit, AbstractAlgorithm, default_svd_algorithm, + left_null!, left_orth!, left_polar!, lq_compact!, null_truncation_strategy, qr_compact!, + right_null!, right_orth!, right_polar!, select_algorithm, svd_compact! -function MatrixAlgebraKit.initialize_output( - ::typeof(left_orth!), A::AbstractBlockSparseMatrix - ) - return nothing +using MatrixAlgebraKit: LeftOrthAlgorithm +for kind in ("polar", "qr", "svd") + @eval begin + function MatrixAlgebraKit.initialize_output( + ::typeof(left_orth!), A::AbstractBlockSparseMatrix, + alg::LeftOrthAlgorithm{Symbol($kind)}, + ) + return nothing + end + end end -function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::AbstractBlockSparseMatrix, F) +function MatrixAlgebraKit.check_input( + ::typeof(left_orth!), A::AbstractBlockSparseMatrix, F, alg::AbstractAlgorithm + ) !isnothing(F) && throw( ArgumentError( "`left_orth!` on block sparse matrices does not support specifying the output" @@ -29,55 +24,60 @@ function MatrixAlgebraKit.check_input(::typeof(left_orth!), A::AbstractBlockSpar return nothing end -function MatrixAlgebraKit.left_orth_qr!(A::AbstractBlockSparseMatrix, F, alg) +using MatrixAlgebraKit: LeftOrthViaQR +function MatrixAlgebraKit.left_orth!(A::AbstractBlockSparseMatrix, F, alg::LeftOrthViaQR) !isnothing(F) && throw( ArgumentError( "`left_orth!` on block sparse matrices does not support specifying the output" ), ) - alg′ = select_algorithm(qr_compact!, A, alg) + alg′ = select_algorithm(qr_compact!, A, alg.alg) return qr_compact!(A, alg′) end -function MatrixAlgebraKit.left_orth_polar!(A::AbstractBlockSparseMatrix, F, alg) +using MatrixAlgebraKit: LeftOrthViaPolar +function MatrixAlgebraKit.left_orth!(A::AbstractBlockSparseMatrix, F, alg::LeftOrthViaPolar) !isnothing(F) && throw( ArgumentError( "`left_orth!` on block sparse matrices does not support specifying the output" ), ) - alg′ = select_algorithm(left_polar!, A, alg) + alg′ = select_algorithm(left_polar!, A, alg.alg) return left_polar!(A, alg′) end -function MatrixAlgebraKit.left_orth_svd!( - A::AbstractBlockSparseMatrix, F, alg, trunc::Nothing = nothing +using MatrixAlgebraKit: LeftOrthViaSVD, does_truncate +function MatrixAlgebraKit.left_orth!( + A::AbstractBlockSparseMatrix, F, alg::LeftOrthViaSVD ) !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) + U, S, Vᴴ = if !does_truncate(alg.alg) + alg′ = select_algorithm(svd_compact!, A, alg.alg) + svd_compact!(A, alg′) + else + alg′ = select_algorithm(svd_compact!, A, alg.alg) + alg_trunc = select_algorithm(svd_trunc!, A, alg′) + svd_trunc!(A, alg_trunc) + end return U, S * Vᴴ end -function MatrixAlgebraKit.initialize_output( - ::typeof(right_orth!), A::AbstractBlockSparseMatrix - ) - return nothing +using MatrixAlgebraKit: RightOrthAlgorithm +for kind in ("lq", "polar", "svd") + @eval begin + function MatrixAlgebraKit.initialize_output( + ::typeof(right_orth!), A::AbstractBlockSparseMatrix, + alg::RightOrthAlgorithm{Symbol($kind)}, + ) + return nothing + end + end end function MatrixAlgebraKit.check_input( - ::typeof(right_orth!), A::AbstractBlockSparseMatrix, F::Nothing + ::typeof(right_orth!), A::AbstractBlockSparseMatrix, F::Nothing, + alg::AbstractAlgorithm, ) !isnothing(F) && throw( ArgumentError( @@ -87,65 +87,63 @@ function MatrixAlgebraKit.check_input( return nothing end -function MatrixAlgebraKit.right_orth_lq!(A::AbstractBlockSparseMatrix, F, alg) +using MatrixAlgebraKit: RightOrthViaLQ +function MatrixAlgebraKit.right_orth!(A::AbstractBlockSparseMatrix, F, alg::RightOrthViaLQ) !isnothing(F) && throw( ArgumentError( "`right_orth!` on block sparse matrices does not support specifying the output" ), ) - alg′ = select_algorithm(lq_compact!, A, alg) + alg′ = select_algorithm(lq_compact!, A, alg.alg) return lq_compact!(A, alg′) end -function MatrixAlgebraKit.right_orth_polar!(A::AbstractBlockSparseMatrix, F, alg) +using MatrixAlgebraKit: RightOrthViaPolar +function MatrixAlgebraKit.right_orth!( + A::AbstractBlockSparseMatrix, F, alg::RightOrthViaPolar + ) !isnothing(F) && throw( ArgumentError( "`right_orth!` on block sparse matrices does not support specifying the output" ), ) - alg′ = select_algorithm(right_polar!, A, alg) + alg′ = select_algorithm(right_polar!, A, alg.alg) return right_polar!(A, alg′) end -function MatrixAlgebraKit.right_orth_svd!( - A::AbstractBlockSparseMatrix, F, alg, trunc::Nothing = nothing +using MatrixAlgebraKit: RightOrthViaSVD +function MatrixAlgebraKit.right_orth!( + A::AbstractBlockSparseMatrix, F, alg::RightOrthViaSVD ) !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) + U, S, Vᴴ = if !does_truncate(alg.alg) + alg′ = select_algorithm(svd_compact!, A, alg.alg) + svd_compact!(A, alg′) + else + alg′ = select_algorithm(svd_compact!, A, alg.alg) + alg_trunc = select_algorithm(svd_trunc!, A, alg′) + svd_trunc!(A, alg_trunc) + end return U * S, Vᴴ end function MatrixAlgebraKit.initialize_output( - ::typeof(left_null!), A::AbstractBlockSparseMatrix + ::typeof(left_null!), A::AbstractBlockSparseMatrix, alg::AbstractAlgorithm ) return nothing end function MatrixAlgebraKit.check_input( - ::typeof(left_null!), A::AbstractBlockSparseMatrix, N::Nothing + ::typeof(left_null!), A::AbstractBlockSparseMatrix, N::Nothing, + alg::AbstractAlgorithm, ) return nothing end -function MatrixAlgebraKit.left_null_qr!(A::AbstractBlockSparseMatrix, N, alg) - return left_null_svd!(A, N, default_svd_algorithm(A)) -end -function MatrixAlgebraKit.left_null_svd!( - A::AbstractBlockSparseMatrix, N, alg, trunc::Nothing +function MatrixAlgebraKit.left_null!( + A::AbstractBlockSparseMatrix, N, alg::AbstractAlgorithm ) - return left_null_svd!(A, N, alg, null_truncation_strategy(; atol = 0, rtol = 0)) + return error("Not implemented.") end function MatrixAlgebraKit.truncate( ::typeof(left_null!), @@ -156,16 +154,19 @@ function MatrixAlgebraKit.truncate( end function MatrixAlgebraKit.initialize_output( - ::typeof(right_null!), A::AbstractBlockSparseMatrix + ::typeof(right_null!), A::AbstractBlockSparseMatrix, alg::AbstractAlgorithm ) return nothing end function MatrixAlgebraKit.check_input( - ::typeof(right_null!), A::AbstractBlockSparseMatrix, N::Nothing + ::typeof(right_null!), A::AbstractBlockSparseMatrix, N::Nothing, + alg::AbstractAlgorithm, ) return nothing end -function MatrixAlgebraKit.right_null_lq!(A::AbstractBlockSparseMatrix, N, alg) +function MatrixAlgebraKit.right_null!( + A::AbstractBlockSparseMatrix, N, alg::AbstractAlgorithm + ) return error("Not implement.") end function MatrixAlgebraKit.truncate( diff --git a/src/factorizations/polar.jl b/src/factorizations/polar.jl index c217600..2dd607e 100644 --- a/src/factorizations/polar.jl +++ b/src/factorizations/polar.jl @@ -27,7 +27,7 @@ end function MatrixAlgebraKit.left_polar!(A::AbstractBlockSparseMatrix, alg::PolarViaSVD) check_input(left_polar!, A) # TODO: Use more in-place operations here, avoid `copy`. - U, S, Vᴴ = svd_compact!(A, alg.svdalg) + U, S, Vᴴ = svd_compact!(A, alg.svd_alg) W = U * Vᴴ # TODO: `copy` is required for now because of: # https://github.com/ITensor/BlockSparseArrays.jl/issues/24 @@ -38,7 +38,7 @@ end function MatrixAlgebraKit.right_polar!(A::AbstractBlockSparseMatrix, alg::PolarViaSVD) check_input(right_polar!, A) # TODO: Use more in-place operations here, avoid `copy`. - U, S, Vᴴ = svd_compact!(A, alg.svdalg) + U, S, Vᴴ = svd_compact!(A, alg.svd_alg) Wᴴ = U * Vᴴ # TODO: `copy` is required for now because of: # https://github.com/ITensor/BlockSparseArrays.jl/issues/24 diff --git a/src/factorizations/svd.jl b/src/factorizations/svd.jl index d3c2fa2..373276e 100644 --- a/src/factorizations/svd.jl +++ b/src/factorizations/svd.jl @@ -1,6 +1,6 @@ using DiagonalArrays: diagonaltype -using MatrixAlgebraKit: - MatrixAlgebraKit, check_input, default_svd_algorithm, svd_compact!, svd_full!, svd_vals! +using MatrixAlgebraKit: MatrixAlgebraKit, check_input, default_svd_algorithm, svd_compact!, + svd_full!, svd_vals! using TypeParameterAccessors: realtype function MatrixAlgebraKit.default_svd_algorithm( @@ -239,3 +239,14 @@ function MatrixAlgebraKit.svd_vals!( end return S end + +# Computing truncation errors. Currently not supported for block sparse matrices, +# needs to be implemented. +struct BlockSparseDiagView{T, A <: AbstractBlockSparseMatrix{T}} <: + AbstractBlockSparseVector{T} + a::A +end +MatrixAlgebraKit.diagview(a::AbstractBlockSparseMatrix) = BlockSparseDiagView(a) +# TODO: Implement this, and/or develop `BlockSparseDiagView` so that the generic +# version in MatrixAlgebraKit works. +MatrixAlgebraKit.truncation_error!(a::BlockSparseDiagView, ind) = nothing diff --git a/src/factorizations/utility.jl b/src/factorizations/utility.jl index f0e7323..fbbf9d4 100644 --- a/src/factorizations/utility.jl +++ b/src/factorizations/utility.jl @@ -1,5 +1,22 @@ +using LinearAlgebra: LinearAlgebra using MatrixAlgebraKit: MatrixAlgebraKit +# Friendlier definitions for block sparse matrices. MatrixAlgebraKit.jl v0.6 +# defines them in terms of `diagview`: +# https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/blob/v0.6.0/src/common/matrixproperties.jl#L40-L64 +# which isn't well supported right now for block sparse matrices. +# TODO: Improve `diagview` for block sparse matrices and remove these definitions. +function MatrixAlgebraKit.is_left_isometric( + A::AnyAbstractBlockSparseMatrix; isapprox_kwargs... + ) + return isapprox(A' * A, LinearAlgebra.I; isapprox_kwargs...) +end +function MatrixAlgebraKit.is_right_isometric( + A::AnyAbstractBlockSparseMatrix; isapprox_kwargs... + ) + return MatrixAlgebraKit.is_left_isometric(A'; isapprox_kwargs...) +end + function infimum(r1::AbstractUnitRange, r2::AbstractUnitRange) (isone(first(r1)) && isone(first(r2))) || throw(ArgumentError("infimum only defined for ranges starting at 1")) diff --git a/test/Project.toml b/test/Project.toml index 391c540..65a615d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -31,7 +31,7 @@ DiagonalArrays = "0.3" GPUArraysCore = "0.2" JLArrays = "0.2, 0.3" LinearAlgebra = "1" -MatrixAlgebraKit = "0.5" +MatrixAlgebraKit = "0.6" Random = "1" SafeTestsets = "0.1" SparseArraysBase = "0.7" diff --git a/test/test_abstract_blocktype.jl b/test/test_abstract_blocktype.jl index 297b54f..94d5750 100644 --- a/test/test_abstract_blocktype.jl +++ b/test/test_abstract_blocktype.jl @@ -3,25 +3,9 @@ using BlockArrays: Block using BlockSparseArrays: BlockSparseMatrix, blockstoredlength using JLArrays: JLArray using LinearAlgebra: hermitianpart, norm -using MatrixAlgebraKit: - eig_full, - eig_trunc, - eig_vals, - eigh_full, - eigh_trunc, - eigh_vals, - isisometry, - left_orth, - left_polar, - lq_compact, - lq_full, - qr_compact, - qr_full, - right_orth, - right_polar, - svd_compact, - svd_full, - svd_trunc +using MatrixAlgebraKit: eig_full, eig_trunc, eig_vals, eigh_full, eigh_trunc, eigh_vals, + isisometric, left_orth, left_polar, lq_compact, lq_full, qr_compact, qr_full, + right_orth, right_polar, svd_compact, svd_full, svd_trunc using SparseArraysBase: storedlength using Test: @test, @test_broken, @testset @@ -106,31 +90,23 @@ arrayts = (Array, JLArray) a = BlockSparseMatrix{elt, AbstractMatrix{elt}}(undef, [2, 3], [2, 3]) a[Block(1, 1)] = dev(randn(elt, 2, 2)) for f in (left_orth, left_polar, qr_compact, qr_full) - if arrayt ≢ Array && f ≡ left_orth - @test_broken f(a) + u, c = f(a) + @test u * c ≈ a + if arrayt ≡ Array + @test isisometric(u; side = :left) else - u, c = f(a) - @test u * c ≈ a - if arrayt ≡ Array - @test isisometry(u; side = :left) - else - # TODO: Fix comparison with UniformScaling on GPU. - @test_broken isisometry(u; side = :left) - end + # TODO: Fix comparison with UniformScaling on GPU. + @test_broken isisometric(u; side = :left) end end for f in (right_orth, right_polar, lq_compact, lq_full) - if arrayt ≢ Array && f ≡ right_orth - @test_broken f(a) + c, u = f(a) + @test c * u ≈ a + if arrayt ≡ Array + @test isisometric(u; side = :right) else - c, u = f(a) - @test c * u ≈ a - if arrayt ≡ Array - @test isisometry(u; side = :right) - else - # TODO: Fix comparison with UniformScaling on GPU. - @test_broken isisometry(u; side = :right) - end + # TODO: Fix comparison with UniformScaling on GPU. + @test_broken isisometric(u; side = :right) end end for f in (svd_compact, svd_full, svd_trunc) diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index a9188c2..c29ce48 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -351,8 +351,8 @@ end A[Block(1, 1)] = randn(rng, T, 3, 2) A[Block(2, 2)] = randn(rng, T, 4, 3) - for kind in (:polar, :qr, :svd) - U, C = left_orth(A; kind) + for alg in (:polar, :qr, :svd) + U, C = left_orth(A; alg) @test U * C ≈ A @test Matrix(U'U) ≈ LinearAlgebra.I end @@ -369,8 +369,8 @@ end A[Block(1, 1)] = randn(rng, T, 2, 3) A[Block(2, 2)] = randn(rng, T, 3, 4) - for kind in (:lq, :polar, :svd) - C, U = right_orth(A; kind) + for alg in (:lq, :polar, :svd) + C, U = right_orth(A; alg) @test C * U ≈ A @test Matrix(U * U') ≈ LinearAlgebra.I end