[ITensors] Port VUMPS functions

[ryanlevy]

This PR adds a few of the functions living in https://github.com/ITensor/ITensorInfiniteMPS.jl/blob/main/src/ITensors.jl

The one thing missing at the moment is the more efficient sqrt(T::DiagTensor) which I couldn't get working, but the eigen implementation seems sane at the moment so there may be no need

[ryanlevy]

Looks like there's some assumption being made that for a sliced mps e.g. psi[2:4] returns a vector of itensors instead of a MPS, not sure if we want to adjust that here or in ITensorInfiniteMPS

[ryanlevy]

Re the sliced MPS, seems the assumption is setting/getting a unit range

function deepcopy_ortho_center(M::AbstractMPS)
  M = copy(M)
  c = ortho_lims(M)
  # TODO: define `getindex(::AbstractMPS, I)` to return `AbstractMPS`
  M[c] = deepcopy(typeof(M)(M[c]))
  return M
end

What are your thoughts given the TODO Matt?

[mtfishman]

I think slicing an MPS/MPO with a unit range should return an MPS/MPO, if we make that change we can just update other code like the one you quote accordingly.

However, the definition in ITensorInfiniteMPS.jl: https://github.com/ITensor/ITensorInfiniteMPS.jl/blob/fc1078a60442081010447caf65c9a22c7278e4ce/src/ITensors.jl#L186 is too naive, as the comment states, since it doesn't preserve information about the orthogonality center.

[ryanlevy]

Just pushed a hopefully better performance version of sqrt. The original logic is a bit strange to me, it will fail for any matrix that isn't PSD without a large atol, how do you feel about allowing for conversion to complex like Base.sqrt does?

julia> eigvals(array(A,i,i'))
2-element Vector{Float64}:
 0.031100580805835387
 1.5475689615743486

julia> sqrt(Hermitian(array(A,i,i'))) ≈ array(sqrt(A),i,i')
true

julia> eigvals(array(B,i,i'))
2-element Vector{Float64}:
 -0.9688994191941646
  0.5475689615743488

jjulia> sqrt(array(B,i,i'))
2×2 Matrix{ComplexF64}:
 0.638411+0.135107im  0.254641-0.338726im
 0.254641-0.338726im  0.101568+0.84922im

julia> sqrt(B)
ERROR: DomainError with -0.9688994191941646:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] sqrt
    @ ./math.jl:608 [inlined]
...