Add tests for svd_trunc (#32)

Author: mtfishman

Project.toml

@@ -1,7 +1,7 @@
 name = "KroneckerArrays"
 uuid = "05d0b138-81bc-4ff7-84be-08becefb1ccc"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.1.25"
+version = "0.1.26"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -23,7 +23,7 @@ KroneckerArraysBlockSparseArraysExt = ["BlockArrays", "BlockSparseArrays"]
 [compat]
 Adapt = "4.3"
 BlockArrays = "1.6"
-BlockSparseArrays = "0.8"
+BlockSparseArrays = "0.8.1"
 DerivableInterfaces = "0.5"
 DiagonalArrays = "0.3.5"
 FillArrays = "1.13"

src/kroneckerarray.jl

@@ -2,6 +2,17 @@
 function _convert(A::Type{<:AbstractArray}, a::AbstractArray)
   return convert(A, a)
 end
+# Custom `_convert` works around the issue that
+# `convert(::Type{<:Diagonal}, ::AbstractMatrix)` isnt' defined
+# in Julia v1.10 (https://github.com/JuliaLang/julia/pull/48895,
+# https://github.com/JuliaLang/julia/pull/52487).
+# TODO: Delete once we drop support for Julia v1.10.
+using LinearAlgebra: LinearAlgebra, Diagonal, diag, isdiag
+_construct(A::Type{<:Diagonal}, a::AbstractMatrix) = A(diag(a))
+function _convert(A::Type{<:Diagonal}, a::AbstractMatrix)
+  LinearAlgebra.checksquare(a)
+  return isdiag(a) ? _construct(A, a) : throw(InexactError(:convert, A, a))
+end
 
 struct KroneckerArray{T,N,A<:AbstractArray{T,N},B<:AbstractArray{T,N}} <: AbstractArray{T,N}
   a::A
@@ -39,7 +50,7 @@ function Base.copyto!(dest::KroneckerArray, src::KroneckerArray)
 end
 
 function Base.convert(::Type{KroneckerArray{T,N,A,B}}, a::KroneckerArray) where {T,N,A,B}
-  return KroneckerArray(convert(A, arg1(a)), convert(B, arg2(a)))
+  return KroneckerArray(_convert(A, arg1(a)), _convert(B, arg2(a)))
 end
 
 # Like `similar` but allows some custom behavior, such as for `FillArrays.Eye`.

test/Project.toml

@@ -21,7 +21,7 @@ TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 Adapt = "4"
 Aqua = "0.8"
 BlockArrays = "1.6"
-BlockSparseArrays = "0.8"
+BlockSparseArrays = "0.8.1"
 DerivableInterfaces = "0.5"
 DiagonalArrays = "0.3.7"
 FillArrays = "1"

test/test_blocksparsearrays.jl

@@ -6,7 +6,8 @@ using FillArrays: Eye, SquareEye
 using JLArrays: JLArray
 using KroneckerArrays: KroneckerArray, ⊗, ×, arg1, arg2
 using LinearAlgebra: norm
-using MatrixAlgebraKit: svd_compact
+using MatrixAlgebraKit: svd_compact, svd_trunc
+using StableRNGs: StableRNG
 using Test: @test, @test_broken, @testset
 using TestExtras: @constinferred
 
@@ -273,6 +274,23 @@ end
   # Broken operations
   @test_broken a[Block.(1:2), Block(2)]
 
+  # svd_trunc
+  dev = adapt(arrayt)
+  r = @constinferred blockrange([2 × 2, 3 × 3])
+  rng = StableRNG(1234)
+  d = Dict(
+    Block(1, 1) => Eye{elt}(2, 2) ⊗ randn(rng, elt, 2, 2),
+    Block(2, 2) => Eye{elt}(3, 3) ⊗ randn(rng, elt, 3, 3),
+  )
+  a = @constinferred dev(blocksparse(d, r, r))
+  if arrayt === Array
+    u, s, v = svd_trunc(a; trunc=(; maxrank=6))
+    u′, s′, v′ = svd_trunc(Matrix(a); trunc=(; maxrank=5))
+    @test Matrix(u * s * v) ≈ u′ * s′ * v′
+  else
+    @test_broken svd_trunc(a; trunc=(; maxrank=6))
+  end
+
   @testset "Block deficient" begin
     da = Dict(Block(1, 1) => Eye{elt}(2, 2) ⊗ dev(randn(elt, 2, 2)))
     a = @constinferred dev(blocksparse(da, r, r))