Memory error (bad_alloc) with MPS/density matrix

ITensor C++ Questions
,
1

Hello all,

I’m working on a Guassian Boson Sampling setup from scratch, where I apply beam splitter 2-mode gates to an MPS state, and calculate the density matrix at the end to make local measurements. I am following this page from the official documentation for the details.

Below is the code in question:

    MPS psi = mps; // Make a local copy
    psi.position(site); // Set the ortho center
    auto psidag = dag(psi); // Create psi bra
    println(psi, psidag);

    psidag.prime("Link"); // Prime bra link indices to contract between them

    // ... Create the reduced density matrix ...
    auto rho = psi(1) * prime(psidag(1), "Site"); // Initialize the tensor
    for (int j = 2; j <= len; j++) { // Contract physical sites priming the site index
        rho *= psi(j);
        rho *= prime(psidag(j), "Site");
    }

When the number of physical sites (len in the code) is 2, everything works fine. However, as soon as I increase it just to 3, I get the following error:

terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc
Aborted (core dumped)

I managed to pinpoint the faulty lines to the two inside the for cycle, but I don’t understand the reason and cannot find a workaround. The local dimension for each site is 9, and the typical bond dimension is around 5, so it should definitely be manageable.

Thank you very much in advance