Fix sqrt decomp (#1669)

Co-authored-by: Matt Fishman <[email protected]>

Author: JoeyT1994

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensors"
 uuid = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 authors = ["Matthew Fishman <[email protected]>", "Miles Stoudenmire <[email protected]>"]
-version = "0.9.9"
+version = "0.9.10"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/tensor_operations/matrix_decomposition.jl

@@ -581,8 +581,8 @@ using NDTensors: map_diag!
 function sqrt_decomp(D::ITensor, u::Index, v::Index)
   (storage(D) isa Union{Diag,DiagBlockSparse}) ||
     error("Must be a diagonal matrix ITensor.")
-  sqrtDL = diag_itensor(u, dag(u)')
-  sqrtDR = diag_itensor(v, dag(v)')
+  sqrtDL = adapt(datatype(D), diag_itensor(u, dag(u)'))
+  sqrtDR = adapt(datatype(D), diag_itensor(v, dag(v)'))
   map_diag!(sqrt ∘ abs, sqrtDL, D)
   map_diag!(sqrt ∘ abs, sqrtDR, D)
   δᵤᵥ = copy(D)

test/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
@@ -7,6 +8,7 @@ Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
+JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MappedArrays = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
 NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"

test/base/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"

test/base/test_svd.jl

@@ -1,4 +1,6 @@
-using ITensors
+using Adapt: adapt
+using ITensors: datatype
+using JLArrays: JLArray
 using Test
 using Suppressor
 
@@ -205,9 +207,14 @@ include(joinpath(@__DIR__, "utils", "util.jl"))
     end
 
     for dir in [ITensors.Out, ITensors.In]
-      L, R, spec = ITensors.factorize_svd(A, l1, l2; dir, ortho="none")
-      @test dir == ITensors.dir(commonind(L, R))
-      @test norm(L * R - A) <= 1e-14
+      for arrayt in (Array, JLArray)
+        A′ = adapt(arrayt, A)
+        L, R, spec = ITensors.factorize_svd(A, l1, l2; dir, ortho="none")
+        @test datatype(L) == datatype(A)
+        @test datatype(R) == datatype(A)
+        @test dir == ITensors.dir(commonind(L, R))
+        @test norm(L * R - A) <= 1e-14
+      end
     end
   end