Refactor sqrt_inv_sqrt (#158)

src/ITensorNetworks.jl

@@ -36,7 +36,7 @@ include("caches/beliefpropagationcache.jl")
 include("contraction_tree_to_graph.jl")
 include("gauging.jl")
 include("utils.jl")
-include("ITensorsExt/itensorutils.jl")
+include("ITensorsExtensions/ITensorsExtensions.jl")
 include("solvers/local_solvers/eigsolve.jl")
 include("solvers/local_solvers/exponentiate.jl")
 include("solvers/local_solvers/dmrg_x.jl")

src/ITensorsExtensions/ITensorsExtensions.jl

@@ -12,7 +12,9 @@ using ITensors:
   map_diag,
   noncommonind,
   noprime,
+  replaceind,
   replaceinds,
+  sim,
   space,
   sqrt_decomp
 using ITensors.NDTensors:
@@ -52,16 +54,24 @@ invsqrt_diag(it::ITensor) = map_diag(inv ∘ sqrt, it)
 pinv_diag(it::ITensor) = map_diag(pinv, it)
 pinvsqrt_diag(it::ITensor) = map_diag(pinv ∘ sqrt, it)
 
-function map_itensor(
-  f::Function, A::ITensor, lind=first(inds(A)); regularization=nothing, kwargs...
-)
-  USV = svd(A, lind; kwargs...)
-  U, S, V, spec, u, v = USV
-  S = map_diag(s -> f(s + regularization), S)
-  sqrtDL, δᵤᵥ, sqrtDR = sqrt_decomp(S, u, v)
-  sqrtDR = denseblocks(sqrtDR) * denseblocks(δᵤᵥ)
-  L, R = U * sqrtDL, V * sqrtDR
-  return L * R
+#TODO: Make this work for non-hermitian A
+function eigendecomp(A::ITensor, linds, rinds; ishermitian=false, kwargs...)
+  @assert ishermitian
+  D, U = eigen(A, linds, rinds; ishermitian, kwargs...)
+  ul, ur = noncommonind(D, U), commonind(D, U)
+  Ul = replaceinds(U, vcat(rinds, ur), vcat(linds, ul))
+
+  return Ul, D, dag(U)
+end
+
+function map_eigvals(f::Function, A::ITensor, inds...; ishermitian=false, kwargs...)
+  if isdiag(A)
+    return map_diag(f, A)
+  end
+
+  Ul, D, Ur = eigendecomp(A, inds...; ishermitian, kwargs...)
+
+  return Ul * map_diag(f, D) * Ur
 end
 
 # Analagous to `denseblocks`.