ITensor3 is slower than same code on ITensor2

Question

Hello,

I recently switched from ITensor 2 to version 3 and I noticed that the same code runs noticeably slower on the newer version.
A calculation that was taking ~7 hours in ITensor 2 now takes almost twice as much time.
I am not sure whether it is ITensor that is responsible to this or something in the cluster configuration where I am running the code, but I would appreciate any help or ideas about what might be causing this.

In particular, on the cluster I am using, I only have access to certain gcc compilers, and in order to use c++17 I had to switch from gcc 4.9.3 to gcc9.1.0 (the only one that supports c++17).
In both cases, I am using LAPACK 3.6.1.

Thanks,

Comments

I agree having a gperftools profile of your code would be very clarifying, Ron. One tip that Matt taught me is that you need to make sure to set OMP_NUM_THREADS=1 (shell environment variable) when using it to get meaningful results. (Otherwise it reports a ton of threading functions that are misleading.) But apart from that it's really great at revealing performance bottlenecks.

---

I ran the dmrgj1j2 code with the same modifications as Matt.
Both v2 and v3 seems approximately the same. I dont see the speedup that Matt sees.
 
I will look into gperftools.

---

V2 Code:

#include "itensor/all.h"

using namespace itensor;

int main(int argc, char* argv[])
    {
    int N = 100;

    println("Input J2 value:");
    Real J2 = 0;
    std::cin >> J2;

    //
    // Initialize the site degrees of freedom.
    //
    auto sites = SpinOne(N);

    //
    // Create the Hamiltonian matrix product operator.
    // Here we use the IQMPO class which is an MPO of
    // IQTensors, tensors whose indices are sorted
    // with respect to quantum numbers
    //
    auto ampo = AutoMPO(sites);
    for(int j = 1; j < N; ++j)
        {
        ampo += 0.5,"S+",j,"S-",j+1;
        ampo += 0.5,"S-",j,"S+",j+1;
        ampo +=     "Sz",j,"Sz",j+1;
        }
    for(int j = 1; j < N-1; ++j)
        {
        ampo += 0.5*J2,"S+",j,"S-",j+2;
        ampo += 0.5*J2,"S-",j,"S+",j+2;
        ampo +=     J2,"Sz",j,"Sz",j+2;
        }
    auto H = MPO(ampo);

    // Set the initial wavefunction matrix product state
    // to be a Neel state.
    //
    auto state = InitState(sites);
    for(int i = 1; i <= N; ++i)
        state.set(i,(i%2==1 ? "Up" : "Dn"));

    auto psi = MPS(state);

    //
    // overlap calculates matrix elements of MPO's with respect to MPS's
    // overlap(psi,H,psi) = <psi|H|psi>
    //
    printfln("Initial energy = %.5f",overlap(psi,H,psi));

    //
    // Set the parameters controlling the accuracy of the DMRG
    // calculation for each DMRG sweep. Here less than 10 maxm
    // values are provided, so all remaining sweeps will use the
    // last maxm (= 200).
    //
    auto sweeps = Sweeps(5);
    sweeps.maxm() = 50,50,100,100,200;
    sweeps.cutoff() = 1E-8;
    println(sweeps);

    //
    // Begin the DMRG calculation
    //
    auto energy = dmrg(psi,H,sweeps,"Quiet");

    //
    // Print the final energy reported by DMRG
    //
    printfln("\nGround State Energy = %.10f",energy);
    printfln("\nUsing overlap = %.10f\n", overlap(psi,H,psi) );

    //println("\nTotal QN of Ground State = ",totalQN(psi));

    //
    // Measure S.S on every bond
    //
    for(int b = 1; b < N; ++b)
        {
        psi.position(b);
        auto ketzz = psi.A(b)*psi.A(b+1)*sites.op("Sz",b)*sites.op("Sz",b+1);
        auto ketpm = psi.A(b)*psi.A(b+1)*sites.op("Sp",b)*sites.op("Sm",b+1)*0.5;
        auto ketmp = psi.A(b)*psi.A(b+1)*sites.op("Sm",b)*sites.op("Sp",b+1)*0.5;
        auto bra = dag(psi.A(b)*psi.A(b+1));
        bra.prime(Site);
        auto SdS = (bra*ketzz).real() + (bra*ketpm).real() + (bra*ketmp).real();
        printfln("S.S b %d = %.10f",b,SdS);
        }

    return 0;
    }

Output:

    vN Entropy at center bond b=50 = 1.281288148412
    Eigs at center bond b=50: 0.4629 0.3626 0.0686 0.0426 0.0402 0.0115 0.0094 0.0017
    Largest m during sweep 1/5 was 9
    Largest truncation error: 5.70055e-16
    Energy after sweep 1/5 is -133.047725652099
    Sweep 1/5 CPU time = 0.0630s (Wall time = 0.0677s)

    vN Entropy at center bond b=50 = 1.538710479300
    Eigs at center bond b=50: 0.4453 0.2765 0.1089 0.0847 0.0356 0.0177 0.0088 0.0059 0.0036 0.0027
    Largest m during sweep 2/5 was 50
    Largest truncation error: 2.64072e-05
    Energy after sweep 2/5 is -149.576391371753
    Sweep 2/5 CPU time = 1.439s (Wall time = 1.861s)

    vN Entropy at center bond b=50 = 1.679144043470
    Eigs at center bond b=50: 0.3261 0.2293 0.2017 0.1810 0.0098 0.0086 0.0081 0.0036 0.0033 0.0030
    Largest m during sweep 3/5 was 100
    Largest truncation error: 5.82099e-06
    Energy after sweep 3/5 is -151.121715673937
    Sweep 3/5 CPU time = 17.58s (Wall time = 21.05s)

    vN Entropy at center bond b=50 = 1.720886079247
    Eigs at center bond b=50: 0.3090 0.2085 0.2072 0.2060 0.0101 0.0100 0.0100 0.0033 0.0033 0.0033
    Largest m during sweep 4/5 was 100
    Largest truncation error: 1.34707e-05
    Energy after sweep 4/5 is -151.147226028195
    Sweep 4/5 CPU time = 21.91s (Wall time = 26.23s)

    vN Entropy at center bond b=50 = 1.731658618466
    Eigs at center bond b=50: 0.3080 0.2069 0.2068 0.2068 0.0103 0.0103 0.0103 0.0033 0.0033 0.0033
    Largest m during sweep 5/5 was 200
    Largest truncation error: 7.18139e-07
    Energy after sweep 5/5 is -151.154828974666
    Sweep 5/5 CPU time = 1m, 43.2s (Wall time = 1m, 52.6s)

---

V3 Code:

#include "itensor/all.h"

using namespace itensor;

int main(int argc, char* argv[])
    {
    int N = 100;

    println("Input J2 value:");
    Real J2 = 0;
    std::cin >> J2;

    //
    // Initialize the site degrees of freedom.
    //
    auto sites = SpinOne(N,{"ConserveQNs=",false});

    //
    // Create the Hamiltonian matrix product operator.
    // Here we use the MPO class which is an MPO of
    // IQTensors, tensors whose indices are sorted
    // with respect to quantum numbers
    //
    auto ampo = AutoMPO(sites);
    for(int j = 1; j < N; ++j)
        {
        ampo += 0.5,"S+",j,"S-",j+1;
        ampo += 0.5,"S-",j,"S+",j+1;
        ampo +=     "Sz",j,"Sz",j+1;
        }
    for(int j = 1; j < N-1; ++j)
        {
        ampo += 0.5*J2,"S+",j,"S-",j+2;
        ampo += 0.5*J2,"S-",j,"S+",j+2;
        ampo +=     J2,"Sz",j,"Sz",j+2;
        }
    auto H = toMPO(ampo);

    // Set the initial wavefunction matrix product state
    // to be a Neel state.
    //
    auto state = InitState(sites);
    for(int i = 1; i <= N; ++i)
        state.set(i,(i%2==1 ? "Up" : "Dn"));

    auto psi0 = MPS(state);

    //
    // inner calculates matrix elements of MPO's with respect to MPS's
    // inner(psi0,H,psi0) = <psi0|H|psi0>
    //
    printfln("Initial energy = %.5f",inner(psi0,H,psi0));

    //
    // Set the parameters controlling the accuracy of the DMRG
    // calculation for each DMRG sweep. Here less than 10 maxdim
    // values are provided, so all remaining sweeps will use the
    // last maxdim (= 200).
    //
    auto sweeps = Sweeps(5);
    sweeps.maxdim() = 50,50,100,100,200;
    sweeps.cutoff() = 1E-8;
    println(sweeps);

    //
    // Begin the DMRG calculation
    //
    auto [energy,psi] = dmrg(H,psi0,sweeps,"Quiet");

    //
    // Print the final energy reported by DMRG
    //
    printfln("\nGround State Energy = %.10f",energy);
    printfln("\nUsing inner = %.10f\n", inner(psi,H,psi) );

    //println("\nTotal QN of Ground State = ",totalQN(psi));

    //
    // Measure S.S on every bond
    //
    for(int b = 1; b < N; ++b)
        {
        psi.position(b);
        auto ketzz = psi(b)*psi(b+1)*op(sites,"Sz",b)*op(sites,"Sz",b+1);
        auto ketpm = psi(b)*psi(b+1)*op(sites,"Sp",b)*op(sites,"Sm",b+1)*0.5;
        auto ketmp = psi(b)*psi(b+1)*op(sites,"Sm",b)*op(sites,"Sp",b+1)*0.5;
        auto bra = dag(psi(b)*psi(b+1));
        bra.prime("Site");
        auto SdS = elt(bra*ketzz) + elt(bra*ketpm) + elt(bra*ketmp);
        printfln("S.S b %d = %.10f",b,SdS);
        }

    return 0;
    }


Output:


    vN Entropy at center bond b=50 = 1.281288148412
    Eigs at center bond b=50: 0.4629 0.3626 0.0686 0.0426 0.0402 0.0115 0.0094 0.0017
    Largest link dim during sweep 1/5 was 9
    Largest truncation error: 3.99609e-16
    Energy after sweep 1/5 is -133.047725652099
    Sweep 1/5 CPU time = 0.109s (Wall time = 0.112s)

    vN Entropy at center bond b=50 = 1.538717778431
    Eigs at center bond b=50: 0.4453 0.2765 0.1089 0.0847 0.0356 0.0177 0.0088 0.0059 0.0036 0.0027
    Largest link dim during sweep 2/5 was 50
    Largest truncation error: 1.19503e-05
    Energy after sweep 2/5 is -149.576253943409
    Sweep 2/5 CPU time = 1.534s (Wall time = 1.952s)

    vN Entropy at center bond b=50 = 1.679145692220
    Eigs at center bond b=50: 0.3261 0.2293 0.2017 0.1810 0.0098 0.0086 0.0081 0.0036 0.0033 0.0030
    Largest link dim during sweep 3/5 was 100
    Largest truncation error: 1.86552e-06
    Energy after sweep 3/5 is -151.121709478775
    Sweep 3/5 CPU time = 18.79s (Wall time = 22.47s)

    vN Entropy at center bond b=50 = 1.720891954614
    Eigs at center bond b=50: 0.3090 0.2085 0.2072 0.2060 0.0101 0.0100 0.0100 0.0033 0.0033 0.0033
    Largest link dim during sweep 4/5 was 100
    Largest truncation error: 4.20862e-06
    Energy after sweep 4/5 is -151.147225986774
    Sweep 4/5 CPU time = 22.59s (Wall time = 27.49s)

    vN Entropy at center bond b=50 = 1.731658867223
    Eigs at center bond b=50: 0.3080 0.2069 0.2068 0.2068 0.0103 0.0103 0.0103 0.0033 0.0033 0.0033
    Largest link dim during sweep 5/5 was 200
    Largest truncation error: 2.21218e-07
    Energy after sweep 5/5 is -151.154828916378
    Sweep 5/5 CPU time = 1m, 45.5s (Wall time = 1m, 54.5s)

---

It looks like you may be running MKL with a single thread (since your CPU time and Wall time are about equal). When I run with a single thread I also get a similar time for ITensor v2 and v3. For example:

ITensor v2 with a single thread (export OMP_NUM_THREADS=1):

    vN Entropy at center bond b=50 = 1.731717097752
    Eigs at center bond b=50: 0.3080 0.2068 0.2068 0.2068 0.0103 0.0103 0.0103 0.0033 0.0033 0.0033
    Largest m during sweep 5/5 was 200
    Largest truncation error: 7.20283e-07
    Energy after sweep 5/5 is -151.154829221677
    Sweep 5/5 CPU time = 42.71s (Wall time = 1m, 0.0s)

ITensor v3 with a single thread:

    vN Entropy at center bond b=50 = 1.731717627565
    Eigs at center bond b=50: 0.3080 0.2068 0.2068 0.2068 0.0103 0.0103 0.0103 0.0033 0.0033 0.0033
    Largest link dim during sweep 5/5 was 200
    Largest truncation error: 2.21877e-07
    Energy after sweep 5/5 is -151.154829255446
    Sweep 5/5 CPU time = 39.90s (Wall time = 54.32s)

In this case, it seems as though using only a single BLAS thread slows down the matrix multiplication and therefore the optimization from ITensor v2 to v3 is not as pronounced (since the optimization was in speeding up the allocation of ITensor data, so it is a subleading cost compared to the matrix multiplication performed in the contraction). You should make sure you are running MKL with multithreading turned on for optimal performance (and note the MKL performance may depend on the machine you are running on).

Likely the slowdown in your code is not due to standard tensor operations like svd and contraction (since you are seeing that dmrg is comparable between ITensor v2 and v3, and dmrg is dominated by tensor contractions and svd operations). I think the next step would be to do some timing and profiling to try to pinpoint what operations are dominating your calculation.

Cheers,
Matt