Merge pull request #257 from LHerviou/Temp_ImprovingInfiniteMPOMatrix

Temp improving infinite mpo matrix

src/ITensorInfiniteMPS.jl

@@ -35,7 +35,7 @@ include("celledvectors.jl")
 include("abstractinfinitemps.jl")
 include("infinitemps.jl")
 include("infinitecanonicalmps.jl")
-include("infinitempomatrix.jl")
+include("infinitempoblock.jl")
 include("infinitempo.jl")
 include("transfermatrix.jl")
 include("models/models.jl")
@@ -60,7 +60,7 @@ export Cell,
   InfMPS,
   InfiniteSum,
   InfiniteMPO,
-  InfiniteMPOMatrix,
+  InfiniteBlockMPO,
   InfiniteSumLocalOps,
   ITensorMap,
   ITensorNetwork,

src/infinitempo.jl

@@ -60,11 +60,11 @@ function cat_to_itensor(Hm::Matrix{ITensor})
   return H, new_l, new_rs[1]
 end
 #
-#  Hm is the InfiniteMPOMatrix
+#  Hm is the InfiniteBlockMPO
 #  Hlrs is an array of {ITensor,Index,Index}s, one for each site in the unit cell.
 #  Hi is a CelledVector of ITensors.
 #
-function InfiniteMPO(Hm::InfiniteMPOMatrix)
+function InfiniteMPO(Hm::InfiniteBlockMPO)
   Hlrs = cat_to_itensor.(Hm) #return an array of {ITensor,Index,Index}
   #
   #  Unpack the array of tuples into three arrays.  And also get an array site indices.
@@ -115,8 +115,8 @@ function fuse_legs!(Hcl::InfiniteMPO, L, R)
       bottom => uniqueinds(bottom, dag(right_link));
       tags=[tags(right_link)],
     )
-    Hcl.data.data[j] = Hcl[j] * comb
-    Hcl.data.data[j + 1] = dag(comb) * Hcl[j + 1]
+    Hcl[j] = Hcl[j] * comb
+    Hcl[j + 1] = dag(comb) * Hcl[j + 1]
   end
   right_link = only(commoninds(Hcl[N], Hcl[N + 1]))
   right_link2 = translatecell(translator(Hcl), right_link, -1)
@@ -135,8 +135,8 @@ function fuse_legs!(Hcl::InfiniteMPO, L, R)
   )
   comb2 = translatecell(translator(Hcl), comb, -1)
   comb2_ind = translatecell(translator(Hcl), comb2, -1)
-  Hcl.data.data[N] = Hcl[N] * comb
-  Hcl.data.data[1] = dag(comb2) * Hcl[1]
+  Hcl[N] = Hcl[N] * comb
+  Hcl[1] = dag(comb2) * Hcl[1]
   L = L * comb2
   R = dag(comb) * R
   return L, R

src/infinitempomatrix.jl → src/infinitempoblock.jl

@@ -1,16 +1,16 @@
 
-mutable struct InfiniteMPOMatrix <: AbstractInfiniteMPS
+mutable struct InfiniteBlockMPO <: AbstractInfiniteMPS
   data::CelledVector{Matrix{ITensor}}
   llim::Int #RealInfinity
   rlim::Int #RealInfinity
   reverse::Bool
 end
 
-translator(mpo::InfiniteMPOMatrix) = mpo.data.translator
-data(mpo::InfiniteMPOMatrix) = mpo.data
+translator(mpo::InfiniteBlockMPO) = mpo.data.translator
+data(mpo::InfiniteBlockMPO) = mpo.data
 
 # TODO better printing?
-function Base.show(io::IO, M::InfiniteMPOMatrix)
+function Base.show(io::IO, M::InfiniteBlockMPO)
   print(io, "$(typeof(M))")
   (length(M) > 0) && print(io, "\n")
   for i in eachindex(M)
@@ -19,8 +19,8 @@ function Base.show(io::IO, M::InfiniteMPOMatrix)
     else
       A = M[i]
       println(io, "Matrix tensor of size $(size(A))")
-      for k in 1:size(A)[1], l in 1:size(A)[2]
-        if !isassigned(A, k + (size(A)[1] - 1) * l)
+      for k in 1:size(A, 1), l in 1:size(A, 2)
+        if !isassigned(A, k + (size(A, 1) - 1) * l)
           println(io, "[($k, $l)] #undef")
         elseif isempty(A[k, l])
           println(io, "[($k, $l)] empty")
@@ -32,27 +32,27 @@ function Base.show(io::IO, M::InfiniteMPOMatrix)
   end
 end
 
-function getindex(ψ::InfiniteMPOMatrix, n::Integer)
+function getindex(ψ::InfiniteBlockMPO, n::Integer)
   return ψ.data[n]
 end
 
-function InfiniteMPOMatrix(arrMat::Vector{Matrix{ITensor}})
-  return InfiniteMPOMatrix(arrMat, 0, size(arrMat)[1], false)
+function InfiniteBlockMPO(arrMat::Vector{Matrix{ITensor}})
+  return InfiniteBlockMPO(arrMat, 0, size(arrMat, 1), false)
 end
 
-function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function)
-  return InfiniteMPOMatrix(CelledVector(data, translator), 0, size(data)[1], false)
+function InfiniteBlockMPO(data::Vector{Matrix{ITensor}}, translator::Function)
+  return InfiniteBlockMPO(CelledVector(data, translator), 0, size(data, 1), false)
 end
 
-function InfiniteMPOMatrix(data::CelledVector{Matrix{ITensor}}, m::Int64, n::Int64)
-  return InfiniteMPOMatrix(data, m, n, false)
+function InfiniteBlockMPO(data::CelledVector{Matrix{ITensor}}, m::Int64, n::Int64)
+  return InfiniteBlockMPO(data, m, n, false)
 end
 
-function InfiniteMPOMatrix(data::CelledVector{Matrix{ITensor}})
-  return InfiniteMPOMatrix(data, 0, size(data)[1], false)
+function InfiniteBlockMPO(data::CelledVector{Matrix{ITensor}})
+  return InfiniteBlockMPO(data, 0, size(data, 1), false)
 end
 
-function ITensors.siteinds(A::InfiniteMPOMatrix)
+function ITensors.siteinds(A::InfiniteBlockMPO)
   data = [dag(only(filterinds(uniqueinds(A[1][1, 1], A[2][1, 1]); plev=0)))]
   for x in 2:(nsites(A) - 1)
     append!(
@@ -71,7 +71,7 @@ function ITensors.siteinds(A::InfiniteMPOMatrix)
   return CelledVector(data, translator(A))
 end
 
-function ITensors.splitblocks(H::InfiniteMPOMatrix)
+function ITensors.splitblocks(H::InfiniteBlockMPO)
   H = copy(H)
   N = nsites(H)
   for j in 1:N
@@ -115,22 +115,10 @@ function find_all_links(Hm::Matrix{ITensor})
   return left_links, right_links
 end
 
-# function convert_itensor_to_itensormatrix(tensor; kwargs...)
-#   if order(tensor) == 3
-#     return convert_itensor_3vector(tensor; kwargs...)
-#   elseif order(tensor) == 4
-#     return convert_itensor_33matrix(tensor; kwargs...)
-#   else
-#     error(
-#       "Conversion of ITensor into matrix of ITensor not planned for this type of tensors"
-#     )
-#   end
-# end
-
 """
     local_mpo_block_projectors(is::Index; new_tags = tags(is))
 
-  Build the projectors on the three parts of the itensor used to split a MPO into an InfiniteMPOMatrix
+  Build the projectors on the three parts of the itensor used to split a MPO into an InfiniteBlockMPO
   More precisely, create projectors on the first dimension, the 2:end-1 and the last dimension of the index
   Input: is the Index to split
   Output: the triplet of projectors (first, middle, last)

src/models/models.jl

@@ -87,14 +87,14 @@ function InfiniteMPO(model::Model, s::CelledVector; kwargs...)
 end
 
 function InfiniteMPO(model::Model, s::CelledVector, translator::Function; kwargs...)
-  return InfiniteMPO(InfiniteMPOMatrix(model, s, translator; kwargs...))
+  return InfiniteMPO(InfiniteBlockMPO(model, s, translator; kwargs...))
 end
 
-function InfiniteMPOMatrix(model::Model, s::CelledVector; kwargs...)
-  return InfiniteMPOMatrix(model, s, translator(s); kwargs...)
+function InfiniteBlockMPO(model::Model, s::CelledVector; kwargs...)
+  return InfiniteBlockMPO(model, s, translator(s); kwargs...)
 end
 
-function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function; kwargs...)
+function InfiniteBlockMPO(model::Model, s::CelledVector, translator::Function; kwargs...)
   N = length(s)
   temp_H = InfiniteSum{MPO}(model, s; kwargs...)
   range_H = maximum(nrange(temp_H)) #Should be improved
@@ -168,7 +168,7 @@ function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function;
     #mpos[j] += dense(Hmat) * setelt(ls[j-1] => total_dim) * setelt(ls[j] => total_dim)
   end
   #unify_indices and add virtual indices to the empty tensors
-  mpos = InfiniteMPOMatrix(mpos, translator)
+  mpos = InfiniteBlockMPO(mpos, translator)
   for x in 1:N
     sp = prime(s[x])
     sd = dag(s[x])
@@ -195,25 +195,25 @@ function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function;
     for j in 2:(size(mpos[x], 1) - 1)
       for k in 2:(size(mpos[x], 2) - 1)
         if isempty(mpos[x][j, k])
-          mpos.data.data[x][j, k] = ITensor(left_inds[j - 1], sd, sp, new_right_inds[k - 1])
+          mpos[x][j, k] = ITensor(left_inds[j - 1], sd, sp, new_right_inds[k - 1])
         else
-          replaceinds!(mpos.data.data[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
+          replaceinds!(mpos[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
         end
       end
     end
     for j in [1, size(mpos[x], 1)]
       for k in 2:(size(mpos[x], 2) - 1)
         if isempty(mpos[x][j, k])
-          mpos.data.data[x][j, k] = ITensor(sd, sp, new_right_inds[k - 1])
+          mpos[x][j, k] = ITensor(sd, sp, new_right_inds[k - 1])
         else
-          replaceinds!(mpos.data.data[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
+          replaceinds!(mpos[x][j, k], right_inds[k - 1] => new_right_inds[k - 1])
         end
       end
     end
     for j in 2:(size(mpos[x], 1) - 1)
       for k in [1, size(mpos[x], 2)]
         if isempty(mpos[x][j, k])
-          mpos.data.data[x][j, k] = ITensor(left_inds[j - 1], sd, sp)
+          mpos[x][j, k] = ITensor(left_inds[j - 1], sd, sp)
         end
       end
     end

src/subspace_expansion.jl

@@ -118,7 +118,7 @@ function generate_twobody_nullspace(
 end
 
 function generate_twobody_nullspace(
-  ψ::InfiniteCanonicalMPS, H::InfiniteMPOMatrix, b::Tuple{Int,Int}; atol=1e-2
+  ψ::InfiniteCanonicalMPS, H::InfiniteBlockMPO, b::Tuple{Int,Int}; atol=1e-2
 )
   n_1, n_2 = b
   L, _ = left_environment(H, ψ)

src/vumps_mpo.jl

@@ -1,7 +1,7 @@
 # Struct for use in linear system solver
 struct AOᴸ
   ψ::InfiniteCanonicalMPS
-  H::InfiniteMPOMatrix
+  H::InfiniteBlockMPO
   n::Int
 end
 
@@ -29,7 +29,7 @@ function (A::AOᴸ)(x)
 end
 
 function initialize_left_environment(
-  H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n::Int64; init_last=false
+  H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n::Int64; init_last=false
 )
   dₕ = size(H[n + 1], 1)
   sit = inds(H[n + 1][1, 1])
@@ -49,7 +49,7 @@ function initialize_left_environment(
 end
 
 function apply_local_left_transfer_matrix(
-  Lstart::Vector{ITensor}, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64
+  Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64
 )
   dₕ = length(Lstart)
   ψ′ = dag(ψ)'
@@ -67,7 +67,7 @@ end
 
 # apply the left transfer matrix at position n1 to the vector Lstart considering it at position m, adding to Ltarget
 function apply_local_left_transfer_matrix(
-  Lstart::ITensor, m::Int64, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64;
+  Lstart::ITensor, m::Int64, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64;
 )
   Ltarget = initialize_left_environment(H, ψ, n_1; init_last=false)
   for j in 1:m
@@ -78,7 +78,7 @@ end
 
 #apply the left transfer matrix n1:n1+nsites(ψ)-1
 function apply_left_transfer_matrix(
-  Lstart::ITensor, m::Int64, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64
+  Lstart::ITensor, m::Int64, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64
 )
   Ltarget = apply_local_left_transfer_matrix(Lstart, m, H, ψ, n_1)
   for j in 1:(nsites(ψ) - 1)
@@ -89,7 +89,7 @@ end
 
 # Also input C bond matrices to help compute the right fixed points
 # of ψ (R ≈ C * dag(C))
-function left_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1e-10)
+function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e-10)
   N = nsites(H)
   @assert N == nsites(ψ)
 
@@ -151,7 +151,7 @@ end
 # Struct for use in linear system solver
 struct AOᴿ
   ψ::InfiniteCanonicalMPS
-  H::InfiniteMPOMatrix
+  H::InfiniteBlockMPO
   n::Int
 end
 
@@ -179,7 +179,7 @@ function (A::AOᴿ)(x)
 end
 
 function initialize_right_environment(
-  H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n::Int64; init_first=false
+  H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n::Int64; init_first=false
 )
   dₕ = size(H[n - 1], 2)
   sit = inds(H[n - 1][1, 1])
@@ -199,7 +199,7 @@ function initialize_right_environment(
 end
 
 function apply_local_right_transfer_matrix(
-  Lstart::Vector{ITensor}, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64
+  Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64
 )
   dₕ = length(Lstart)
   ψ′ = dag(ψ)'
@@ -218,7 +218,7 @@ end
 function apply_local_right_transfer_matrix(
   Lstart::ITensor,
   m::Int64,
-  H::InfiniteMPOMatrix,
+  H::InfiniteBlockMPO,
   ψ::InfiniteCanonicalMPS,
   n_1::Int64;
   reset=true,
@@ -236,7 +236,7 @@ end
 
 #apply the right transfer matrix n1:n1+nsites(ψ)-1
 function apply_right_transfer_matrix(
-  Lstart::ITensor, m::Int64, H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS, n_1::Int64
+  Lstart::ITensor, m::Int64, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64
 )
   Ltarget = apply_local_right_transfer_matrix(Lstart, m, H, ψ, n_1)
   for j in 1:(nsites(ψ) - 1)
@@ -245,7 +245,7 @@ function apply_right_transfer_matrix(
   return Ltarget
 end
 
-function right_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1e-10)
+function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e-10)
   N = nsites(H)
   @assert N == nsites(ψ)
 
@@ -303,7 +303,7 @@ function right_environment(H::InfiniteMPOMatrix, ψ::InfiniteCanonicalMPS; tol=1
   return CelledVector(Rs), eᵣ[1]
 end
 
-function vumps(H::InfiniteMPOMatrix, ψ::InfiniteMPS; kwargs...)
+function vumps(H::InfiniteBlockMPO, ψ::InfiniteMPS; kwargs...)
   return vumps(H, orthogonalize(ψ, :); kwargs...)
 end
 
@@ -347,7 +347,7 @@ end
 
 function tdvp_iteration_sequential(
   solver::Function,
-  H::InfiniteMPOMatrix,
+  H::InfiniteBlockMPO,
   ψ::InfiniteCanonicalMPS;
   (ϵᴸ!)=fill(1e-15, nsites(ψ)),
   (ϵᴿ!)=fill(1e-15, nsites(ψ)),
@@ -409,7 +409,7 @@ end
 
 function tdvp_iteration_parallel(
   solver::Function,
-  H::InfiniteMPOMatrix,
+  H::InfiniteBlockMPO,
   ψ::InfiniteCanonicalMPS;
   (ϵᴸ!)=fill(1e-15, nsites(ψ)),
   (ϵᴿ!)=fill(1e-15, nsites(ψ)),