Sure! Below, I paste simple codes (Julia and C++) for subspace expansion a la [Hubig 2015] that uses directsum.
Julia version:
using ITensors
### DMRG just to obtain a psi with large χ
function dmrg()
N = 100
sites = siteinds("S=1/2",N; conserve_qns=true)
os = OpSum()
for j=1:N-1
os += 1, "Sz",j,"Sz",j+1
os += 0.5, "S+",j,"S-",j+1
os += 0.5, "S-",j,"S+",j+1
end
H = MPO(os,sites)
states = [isodd(n) ? "Up" : "Dn" for n in 1:N]
psi0 = MPS(sites, states)
nsweeps = 10
maxdim = [200]
cutoff = [1E-14]
noise = [0]
energy,psi = ITensors.dmrg(H, psi0; nsweeps, maxdim, cutoff, noise = noise);
return H, psi
end
let
H, psi = dmrg()
orthogonalize!(psi, 50) # OC, we will try to enlarge bond(50 <-> 51)
## Left environment at site 50
L = ITensor(true)
for b = 1 : 49
L *= psi[b]
L *= H[b]
L *= dag(prime(psi[b]))
end
## Enlargement term as per [Hubig et. al.]
pad = L * psi[50]
pad *= H[50]
noprime!(pad)
## Get the links right
nextlink = commonind(psi[50], psi[51]; tags = "Link")
nextlinkH = commonind(H[50], H[51]; tags = "Link")
comb = combiner(nextlink, nextlinkH; tags=tags(nextlink))
padind = combinedind(comb)
pad *= comb
## THIS IS THE PROBLEMATIC PART
## more than 20 times slower than the SVD (see below)
@time begin
phienlarged, sumind = directsum(psi[50] => nextlink, pad => padind)
end
## Padding psi[51] with zeros is omitted
## SVD
ij = commoninds(phienlarged, psi[50])
@time begin
U, S, V = svd(phienlarged, ij; maxdim=200);
end
end
C++ version:
#include <itensor/all.h>
using namespace std;
using namespace itensor;
int main() {
//**********************************************************
int N = 100;
auto sites = SpinHalf(N, {"ConserveSz", true});
auto ampo = AutoMPO(sites);
for (int ii = 1; ii < N; ++ii) {
ampo += 0.5, "Sp", ii, "Sm", ii+1;
ampo += 0.5, "Sm", ii, "Sp", ii+1;
ampo += 1.0, "Sz", ii, "Sz", ii+1;
}
auto H = toMPO(ampo);
auto state = InitState(sites);
for(int b = 1; b <= N; ++b) {
if (b % 2 == 1)
state.set(b, "Up");
else
state.set(b, "Dn");
}
auto psi0 = MPS(state);
auto sweeps = Sweeps(10);
sweeps.maxdim() = 200;
sweeps.cutoff() = 1e-14;
auto [energy,psi] = dmrg(H, psi0, sweeps, {"Silent", true});
//**********************************************************
psi.position(50);
auto L = ITensor(1);
for (int b = 1; b <= 49; ++b) {
L *= psi(b);
L *= H(b);
L *= dag(prime(psi(b)));
}
auto pad = L * psi(50);
pad *= H(50);
pad.noPrime();
auto nextlink = commonIndex(psi(50), psi(51), "Link");
auto nextlinkH = commonIndex(H(50), H(51), "Link");
auto [comb, padind] = combiner({nextlink, nextlinkH}, {"Tags", tags(nextlink)});
pad *= comb;
auto [phienlarged, sumind] = directSum(psi(50), pad, nextlink, padind);
// Then SVD and so on...
}
The times I get are as follows (with QN):
Julia directsum => 0.240269 seconds (~ 13% gc time)
C++ directSum => 0.005552 seconds (measured by std::chrono::steady_clock)
Julia SVD => 0.009786 seconds.
And without QN:|
Julia directsum => 0.033996 seconds.
C++ directSum => 0.054598 seconds (surprize !!)
Julia SVD => 0.080785 seconds.
On another note. In C++ I have sparingly used delta() to replace indices. What I discovered that this is super inefficient in Julia (even compared to QN directsum), so switched to replaceind!.