Matrix functions with abstract block type (#155)

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.7.18"
+version = "0.7.19"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

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

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