From 0ca39b0abf7879ca41c5f288a351a153e7f92357 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 20 Jun 2025 14:04:01 -0400 Subject: [PATCH 1/3] Fix matrix functions of block sparse arrays with abstract block type --- src/abstractblocksparsearray/linearalgebra.jl | 20 ++- test/test_factorizations.jl | 128 +++++++++--------- 2 files changed, 83 insertions(+), 65 deletions(-) diff --git a/src/abstractblocksparsearray/linearalgebra.jl b/src/abstractblocksparsearray/linearalgebra.jl index 70702f0c..21694cb9 100644 --- a/src/abstractblocksparsearray/linearalgebra.jl +++ b/src/abstractblocksparsearray/linearalgebra.jl @@ -93,8 +93,14 @@ const MATRIX_FUNCTIONS_UNSTABLE = [ ] function initialize_output_blocksparse(f::F, a::AbstractMatrix) where {F} - B = Base.promote_op(f, blocktype(a)) - return similar(a, BlockType(B)) + blockt = Base.promote_op(f, blocktype(a)) + elt′ = Base.promote_op(f, eltype(a)) + blockt′ = if blockt >: AbstractMatrix{elt′} || blockt === Union{} + AbstractMatrix{elt′} + else + blockt + end + return similar(a, BlockType(blockt′)) end function matrix_function_blocksparse(f::F, a::AbstractMatrix; kwargs...) where {F} @@ -117,8 +123,14 @@ end for f in MATRIX_FUNCTIONS_UNSTABLE @eval begin function initialize_output_blocksparse(::typeof($f), a::AbstractMatrix) - B = similartype(blocktype(a), complex(eltype(a))) - return similar(a, BlockType(B)) + elt′ = complex(eltype(a)) + blockt = Base.promote_op(similar, blocktype(a), elt′) + blockt′ = if blockt >: AbstractMatrix{elt′} || blockt === Union{} + AbstractMatrix{elt′} + else + blockt + end + return similar(a, BlockType(blockt′)) end end end diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl index d9fb6294..5113497b 100644 --- a/test/test_factorizations.jl +++ b/test/test_factorizations.jl @@ -34,79 +34,85 @@ using StableRNGs: StableRNG using Test: @inferred, @test, @test_broken, @test_throws, @testset @testset "Matrix functions (T=$elt)" for elt in (Float32, Float64, ComplexF64) - rng = StableRNG(123) - a = BlockSparseMatrix{elt}(undef, [2, 3], [2, 3]) - a[Block(1, 1)] = randn(rng, elt, 2, 2) - a[Block(2, 2)] = randn(rng, elt, 3, 3) - MATRIX_FUNCTIONS = BlockSparseArrays.MATRIX_FUNCTIONS - MATRIX_FUNCTIONS = [MATRIX_FUNCTIONS; [:inv, :pinv]] - # Only works when real, also isn't defined in Julia 1.10. - MATRIX_FUNCTIONS = setdiff(MATRIX_FUNCTIONS, [:cbrt]) - MATRIX_FUNCTIONS_LOW_ACCURACY = [:acoth] - for f in setdiff(MATRIX_FUNCTIONS, MATRIX_FUNCTIONS_LOW_ACCURACY) - @eval begin - fa = $f($a) - @test Matrix(fa) ≈ $f(Matrix($a)) rtol = √(eps(real($elt))) - @test fa isa BlockSparseMatrix - @test issetequal(eachblockstoredindex(fa), [Block(1, 1), Block(2, 2)]) + for matrixt in (Matrix, AbstractMatrix) + a = BlockSparseMatrix{elt,matrixt{elt}}(undef, [2, 3], [2, 3]) + rng = StableRNG(123) + a[Block(1, 1)] = randn(rng, elt, 2, 2) + a[Block(2, 2)] = randn(rng, elt, 3, 3) + MATRIX_FUNCTIONS = BlockSparseArrays.MATRIX_FUNCTIONS + MATRIX_FUNCTIONS = [MATRIX_FUNCTIONS; [:inv, :pinv]] + # Only works when real, also isn't defined in Julia 1.10. + MATRIX_FUNCTIONS = setdiff(MATRIX_FUNCTIONS, [:cbrt]) + MATRIX_FUNCTIONS_LOW_ACCURACY = [:acoth] + for f in setdiff(MATRIX_FUNCTIONS, MATRIX_FUNCTIONS_LOW_ACCURACY) + @eval begin + fa = $f($a) + @test Matrix(fa) ≈ $f(Matrix($a)) rtol = √(eps(real($elt))) + @test fa isa BlockSparseMatrix + @test issetequal(eachblockstoredindex(fa), [Block(1, 1), Block(2, 2)]) + end end - end - for f in MATRIX_FUNCTIONS_LOW_ACCURACY - @eval begin - fa = $f($a) - if !Sys.isapple() && ($elt <: Real) - # `acoth` appears to be broken on this matrix on Windows and Ubuntu - # for real matrices. - @test_broken Matrix(fa) ≈ $f(Matrix($a)) rtol = √eps(real($elt)) - else - @test Matrix(fa) ≈ $f(Matrix($a)) rtol = √eps(real($elt)) + for f in MATRIX_FUNCTIONS_LOW_ACCURACY + @eval begin + fa = $f($a) + if !Sys.isapple() && ($elt <: Real) + # `acoth` appears to be broken on this matrix on Windows and Ubuntu + # for real matrices. + @test_broken Matrix(fa) ≈ $f(Matrix($a)) rtol = √eps(real($elt)) + else + @test Matrix(fa) ≈ $f(Matrix($a)) rtol = √eps(real($elt)) + end + @test fa isa BlockSparseMatrix + @test issetequal(eachblockstoredindex(fa), [Block(1, 1), Block(2, 2)]) end - @test fa isa BlockSparseMatrix - @test issetequal(eachblockstoredindex(fa), [Block(1, 1), Block(2, 2)]) end end # Catch case of off-diagonal blocks. - rng = StableRNG(123) - a = BlockSparseMatrix{elt}(undef, [2, 3], [2, 3]) - a[Block(1, 1)] = randn(rng, elt, 2, 2) - a[Block(1, 2)] = randn(rng, elt, 2, 3) - for f in MATRIX_FUNCTIONS - @eval begin - @test_throws ArgumentError $f($a) + for matrixt in (Matrix, AbstractMatrix) + a = BlockSparseMatrix{elt,matrixt{elt}}(undef, [2, 3], [2, 3]) + rng = StableRNG(123) + a[Block(1, 1)] = randn(rng, elt, 2, 2) + a[Block(1, 2)] = randn(rng, elt, 2, 3) + for f in BlockSparseArrays.MATRIX_FUNCTIONS + @eval begin + @test_throws ArgumentError $f($a) + end end end # Missing diagonal blocks. - rng = StableRNG(123) - a = BlockSparseMatrix{elt}(undef, [2, 3], [2, 3]) - a[Block(2, 2)] = randn(rng, elt, 3, 3) - MATRIX_FUNCTIONS = BlockSparseArrays.MATRIX_FUNCTIONS - # These functions involve inverses so they break when there are zeros on the diagonal. - MATRIX_FUNCTIONS_SINGULAR = [ - :log, :acsc, :asec, :acot, :acsch, :asech, :acoth, :csc, :cot, :csch, :coth - ] - MATRIX_FUNCTIONS = setdiff(MATRIX_FUNCTIONS, MATRIX_FUNCTIONS_SINGULAR) - # Dense version is broken for some reason, investigate. - MATRIX_FUNCTIONS = setdiff(MATRIX_FUNCTIONS, [:cbrt]) - for f in MATRIX_FUNCTIONS - @eval begin - fa = $f($a) - @test Matrix(fa) ≈ $f(Matrix($a)) rtol = √(eps(real($elt))) - @test fa isa BlockSparseMatrix - @test issetequal(eachblockstoredindex(fa), [Block(1, 1), Block(2, 2)]) + for matrixt in (Matrix, AbstractMatrix) + a = BlockSparseMatrix{elt,matrixt{elt}}(undef, [2, 3], [2, 3]) + rng = StableRNG(123) + a[Block(2, 2)] = randn(rng, elt, 3, 3) + MATRIX_FUNCTIONS = BlockSparseArrays.MATRIX_FUNCTIONS + # These functions involve inverses so they break when there are zeros on the diagonal. + MATRIX_FUNCTIONS_SINGULAR = [ + :log, :acsc, :asec, :acot, :acsch, :asech, :acoth, :csc, :cot, :csch, :coth + ] + MATRIX_FUNCTIONS = setdiff(MATRIX_FUNCTIONS, MATRIX_FUNCTIONS_SINGULAR) + # Dense version is broken for some reason, investigate. + MATRIX_FUNCTIONS = setdiff(MATRIX_FUNCTIONS, [:cbrt]) + for f in MATRIX_FUNCTIONS + @eval begin + fa = $f($a) + @test Matrix(fa) ≈ $f(Matrix($a)) rtol = √(eps(real($elt))) + @test fa isa BlockSparseMatrix + @test issetequal(eachblockstoredindex(fa), [Block(1, 1), Block(2, 2)]) + end end - end - SINGULAR_EXCEPTION = if VERSION < v"1.11-" - # A different exception is thrown in older versions of Julia. - LinearAlgebra.LAPACKException - else - LinearAlgebra.SingularException - end - for f in setdiff(MATRIX_FUNCTIONS_SINGULAR, [:log]) - @eval begin - @test_throws $SINGULAR_EXCEPTION $f($a) + SINGULAR_EXCEPTION = if VERSION < v"1.11-" + # A different exception is thrown in older versions of Julia. + LinearAlgebra.LAPACKException + else + LinearAlgebra.SingularException + end + for f in setdiff(MATRIX_FUNCTIONS_SINGULAR, [:log]) + @eval begin + @test_throws $SINGULAR_EXCEPTION $f($a) + end end end end From da4f3e3e1e85c22b6d015013d4733af0da294cfd Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 20 Jun 2025 14:04:21 -0400 Subject: [PATCH 2/3] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index de1c132e..efc4a6c6 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.7.18" +version = "0.7.19" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From bf06bfc3ecb309f952d9edb263ded5d8bf9eb9d0 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 20 Jun 2025 14:33:25 -0400 Subject: [PATCH 3/3] Fix type comparison --- src/abstractblocksparsearray/linearalgebra.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/abstractblocksparsearray/linearalgebra.jl b/src/abstractblocksparsearray/linearalgebra.jl index 21694cb9..0c9483fb 100644 --- a/src/abstractblocksparsearray/linearalgebra.jl +++ b/src/abstractblocksparsearray/linearalgebra.jl @@ -95,7 +95,7 @@ const MATRIX_FUNCTIONS_UNSTABLE = [ function initialize_output_blocksparse(f::F, a::AbstractMatrix) where {F} blockt = Base.promote_op(f, blocktype(a)) elt′ = Base.promote_op(f, eltype(a)) - blockt′ = if blockt >: AbstractMatrix{elt′} || blockt === Union{} + blockt′ = if !(blockt <: AbstractMatrix{elt′}) || blockt === Union{} AbstractMatrix{elt′} else blockt @@ -125,7 +125,7 @@ for f in MATRIX_FUNCTIONS_UNSTABLE function initialize_output_blocksparse(::typeof($f), a::AbstractMatrix) elt′ = complex(eltype(a)) blockt = Base.promote_op(similar, blocktype(a), elt′) - blockt′ = if blockt >: AbstractMatrix{elt′} || blockt === Union{} + blockt′ = if !(blockt <: AbstractMatrix{elt′}) || blockt === Union{} AbstractMatrix{elt′} else blockt