Minimize the number of swaps when applying nonlocal gates to MPS/MPO

benchmark/bench_dmrg.jl

@@ -45,6 +45,12 @@ let N = 100
   dmrg(H, psi, sweeps; outputlevel=0)
 
   suite["1d_S=1_heisenberg_qn"] = @benchmarkable dmrg($H, $psi, $sweeps; outputlevel=0)
+
+  # Precompile
+  H_splitblocks = splitblocks(linkinds, H)
+  dmrg(H_splitblocks, psi, sweeps; outputlevel=0)
+
+  suite["1d_S=1_heisenberg_qn_splitblocks"] = @benchmarkable dmrg($H_splitblocks, $psi, $sweeps; outputlevel=0)
 end
 
 #let Ny = 6,

docs/src/MPSandMPO.md

@@ -15,10 +15,10 @@ MPS(::Type{<:Number}, ::Vector{<:Index})
 randomMPS(sites::Vector{<:Index})
 randomMPS(::Type{<:Number}, sites::Vector{<:Index})
 randomMPS(::Vector{<:Index}, ::Any)
-productMPS(::Vector{<:Index}, ::Any)
-productMPS(::Type{<:Number}, ::Vector{<:Index}, ::Any)
-productMPS(::Vector{<:IndexVal})
-productMPS(::Type{<:Number}, ::Vector{<:IndexVal})
+MPS(::Vector{<:Index}, ::Any)
+MPS(::Type{<:Number}, ::Vector{<:Index}, ::Any)
+MPS(::Vector{<:IndexVal})
+MPS(::Type{<:Number}, ::Vector{<:IndexVal})
 ```
 
 ## MPO Constructors
@@ -40,6 +40,7 @@ deepcopy(::ITensors.AbstractMPS)
 
 ```@docs
 eltype(::ITensors.AbstractMPS)
+ITensors.promote_itensor_eltype(::ITensors.AbstractMPS)
 flux(::ITensors.AbstractMPS)
 hasqns(::ITensors.AbstractMPS)
 length(::ITensors.AbstractMPS)

docs/src/examples/DMRG.md

@@ -52,7 +52,7 @@ psi0 = randomMPS(sites,2)
 Here we have made a random MPS of bond dimension 2. We could have used a random product
 state instead, but choosing a slightly larger bond dimension can help DMRG avoid getting
 stuck in local minima. We could also set psi to some specific initial state using the 
-function `productMPS`, which is actually required if we were conserving QNs.
+function `MPS`, which is actually required if we were conserving QNs.
 
 Finally, we are ready to call DMRG:
 

docs/src/getting_started/Tutorials.md

@@ -216,7 +216,7 @@ by the lines
 
 ```julia
 state = [isodd(n) ? "Up" : "Dn" for n=1:N]
-psi0 = productMPS(sites,state)
+psi0 = MPS(sites,state)
 ```
 
 The first line of the new code above makes an array of strings which 
@@ -251,21 +251,21 @@ computing the quantum-number flux of `psi0`
     use the initial state:
     ```julia
     state = ["Up" for n=1:N]
-    psi0 = productMPS(sites,state)
+    psi0 = MPS(sites,state)
     ```
     Or to initialize this 10-site system to have a total "Sz" of +16
     in ITensor units (``S^z=8`` in physical units):
     ```julia
     state = ["Dn","Up","Up","Up","Up","Up","Up","Up","Up","Up"]
-    psi0 = productMPS(sites,state)
+    psi0 = MPS(sites,state)
     ```
     would work (as would any `state` with one "Dn" and nine "Up"'s
     in any order).
     Or you could initialize to a total "Sz" of +18
     in ITensor units (``S^z=9`` in physical units) as
     ```julia
     state = ["Z0","Up","Up","Up","Up","Up","Up","Up","Up","Up"]
-    psi0 = productMPS(sites,state)
+    psi0 = MPS(sites,state)
     ```
     where "Z0" refers to the ``S^z=0`` state of a spin-one spin.
 
@@ -294,7 +294,7 @@ let
   H = MPO(ampo,sites)
 
   state = [isodd(n) ? "Up" : "Dn" for n=1:N]
-  psi0 = productMPS(sites,state)
+  psi0 = MPS(sites,state)
   @show flux(psi0)
 
   sweeps = Sweeps(5)
@@ -413,7 +413,7 @@ let
   end
 
   # Initialize psi to be a product state (alternating up and down)
-  psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")
+  psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
 
   c = div(N,2)
 
@@ -490,7 +490,7 @@ we will take to be near the center of the MPS. The details of this
 function are outside the scope of this tutorial, but are explained in the
 example code for measuring MPS.
 
-The line of code `psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")`
+The line of code `psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")`
 initializes our MPS `psi` as a product state of alternating 
 up and down spins. We call `measure_Sz` before starting the
 time evolution.

src/ITensors.jl

@@ -7,6 +7,10 @@ The ITensor library also includes composable and extensible algorithms for optim
 """
 module ITensors
 
+#if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@compiler_options"))
+#  @eval Base.Experimental.@compiler_options compile=min optimize=1 infer=false
+#end
+
 #####################################
 # External packages
 #

src/indexset.jl

@@ -456,28 +456,45 @@ function swaptags(is::Indices, tags1, tags2, args...; kwargs...)
   return swaptags(fmatch(args...; kwargs...), is, tags1, tags2)
 end
 
-function replaceinds(is::Indices, rep_inds::Pair{<:Index,<:Index}...)
-  return replaceinds(is, zip(rep_inds...)...)
-end
+# Check that the QNs are all the same
+hassameflux(i1::Index, i2::Index) = (dim(i1) == dim(i2))
 
-function replaceinds(is::Indices, rep_inds::Vector{<:Pair{<:Index,<:Index}})
-  return replaceinds(is, rep_inds...)
+# For notation:
+# `replaceinds(is, [i, j] => [k, l])`
+# `replaceinds(is, (i, j) => (k, l))`
+function replaceinds(is::Indices, rep_inds::Pair)
+  return replaceinds(is, first(rep_inds), last(rep_inds))
 end
 
-function replaceinds(is::Indices, rep_inds::Tuple{Vararg{Pair{<:Index,<:Index}}})
-  return replaceinds(is, rep_inds...)
+# For notation:
+# `replaceinds(is, i => k, j => l)`
+function replaceinds(is::Indices, rep_inds::Pair{<:Index,<:Index}...)
+  return replaceinds(is, first.(rep_inds), last.(rep_inds))
 end
 
-function replaceinds(is::Indices, rep_inds::Pair)
-  return replaceinds(is, Tuple(first(rep_inds)) .=> Tuple(last(rep_inds)))
+# For notation:
+# `replaceinds(is, i, k)`
+function replaceinds(is::Indices, inds1::Index, inds2::Index)
+  return replaceinds(is, (inds1,), (inds2,))
 end
 
-# Check that the QNs are all the same
-hassameflux(i1::Index, i2::Index) = (dim(i1) == dim(i2))
+# For notation:
+# `replaceinds(is, [i => k, j => l])`
+# `replaceinds(is, (i => k, j => l))`
+function replaceinds(is::Indices, inds1_inds2)
+  @assert all(x -> x isa Pair, inds1_inds2)
+  return replaceinds(is, first.(inds1_inds2), last.(inds1_inds2))
+end
 
+# For notation:
+# `replaceinds(is, [i, j], [k, l])`
+# `replaceinds(is, (i, j), (k, l))`
 function replaceinds(is::Indices, inds1, inds2)
+  @assert length(inds1) == length(inds2)
+  isempty(inds1) && return is
   is1 = (inds1)
   poss = indexin(is1, is)
+  # TODO: don't convert to Tuple in the case of a Vector
   is_tuple = Tuple(is)
   for (j, pos) in enumerate(poss)
     isnothing(pos) && continue