Hi all,

I found that the built in function correlation_matrix() can lead to a Out Of Memory Error easyily even for a small size system with a moderate bond dimensions.

The minimal code can be like

```
ITensors.Strided.set_num_threads(1)
BLAS.set_num_threads(1)
ITensors.enable_threaded_blocksparse()
Nx = 4
Ny = 4
N = Nx * Ny
Npart = Int(Nx * Ny / 2)
t = 1.0
D = 600
Nsweep = 4
sweeps = Sweeps(Nsweep)
maxdim!(sweeps, D)
cutoff!(sweeps, 1E-8)
sites = siteinds("Electron", N; conserve_qns=true)
lattice = square_lattice(Nx, Ny; yperiodic=true)
ampo = AutoMPO()
for b in lattice
ampo .+= -t, "Cdagup", b.s1, "Cup", b.s2
ampo .+= -t, "Cdagup", b.s2, "Cup", b.s1
ampo .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
ampo .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
end
H = MPO(ampo, sites)
H = splitblocks(linkinds, H)
state = ["Emp" for n = 1:N]
p = Npart
for i = N:-1:1
if p > i
println("Doubly occupying site $i")
state[i] = "UpDn"
p -= 2
elseif p > 0
println("Singly occupying site $i")
state[i] = (isodd(i) ? "Up" : "Dn")
p -= 1
end
end
psi0_ini = randomMPS(sites, state, 10)
energy0, psi0 = dmrg(H, psi0_ini, sweeps)
Cupij = correlation_matrix(psi0,"Cdagup","Cup")
```

So for a Hubbard model on a 4*4 cylinder with bond dimension D = 600, I can readily obtain the ground state wavefunction, but get a out of memory error report.

I am not sure whthere I miss something when using the built-in function correlation_matrix().

Any comments will be helpful :)

Xin-Chi