Supported contracting MPO-MPO by fitting algorithm

src/contract_mpo_mps.jl

@@ -30,16 +30,20 @@ function ITensors.contract(
   init_mps = deepcopy(init_mps)
   init_mps = sim(linkinds, init_mps)
   Ai = siteinds(A)
-  ti = Vector{Index}(undef, n)
+  init_mpsi = siteinds(init_mps)
   for j in 1:n
+    ti = nothing
     for i in Ai[j]
       if !hasind(psi0[j], i)
-        ti[j] = i
+        ti = i
         break
       end
     end
+    if ti !== nothing
+      ci = commoninds(init_mpsi[j], A[j])[1]
+      replaceind!(init_mps[j], ci=>ti)
+    end
   end
-  replace_siteinds!(init_mps, ti)
 
   t = Inf
   reverse_step = false

test/test_contract_mpo.jl

@@ -51,4 +51,23 @@ using Test
   @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1E-4
 end
 
+@testset "Contract MPO-MPO" begin
+  nbit = 5
+  sites = siteinds("Qubit", nbit)
+  M1 = randomMPO(sites) + randomMPO(sites)
+  M2 = randomMPO(sites) + randomMPO(sites)
+
+  # The function `apply` does not work correctly with the mapping-MPO-to-MPS trick.
+  M1 = replaceprime(M1, 1=>2, 0=>1)
+
+  M2_ = MPS(length(sites))
+  for n in eachindex(sites)
+    M2_[n] = M2[n]
+  end
+
+  M12_ref = contract(M1, M2; alg="naive")
+  M12 = contract(M1, M2_; alg="fit")
+  @test M12_ref ≈ M12
+end
+
 nothing