# nmultMPO very slow

+1 vote

I'm attempting to multiply two MPOs. I do, in full,

#include "itensor/itensor.h"
#include "itensor/mps/autompo.h"
#include "itensor/mps/sites/spinhalf.h"
#include "itensor/util/print_macro.h"

using namespace itensor;

int main(int argc, char* argv[])
{
int L      = 8;
float hz   = 2.0;

auto sites = SpinHalf(L);
auto H_ampo = AutoMPO(sites);

for(int j = 1; j <=  L; j++) { H_ampo += hz, "Sz",j; }

MPOt<ITensor> H  = MPO(H_ampo);
MPOt<ITensor> Hn = MPO(H_ampo);
MPOt<ITensor> res = MPO(H_ampo);

Print(H.A(2).inds());

nmultMPO<MPOt<ITensor>>(H, Hn, res);

return(0);
}


I compile with a Makefile only slightly modified from the sample---I'd be happy to post it, but it doesn't seem terribly interesting.

After printing something reasonable for the indices of H.A(2), the thing allocates hundreds of mb of memory (per top) and then runs for a long time (tens of minutes) before returning. This is not what I expect for an eight-site, bond-dimension-1 MPO!

I'm sure I'm misunderstanding something, here. What is it?

Hi, sorry about the slow reply! This is a good question: so nmultMPO expects additional named arguments saying what truncation parameters to use. Without them it was keeping exponentially many states as it went from left to right "compressing" the two MPOs.

I just pushed a fix for this in the latest version of ITensor. Now if you don't state any truncation parameters it defaults to a small cutoff of 1E-14.

But the best way to call this algorithm is as:

nmultMPO(A,B,C,{"Cutoff",1E-12});


where you may want to use a larger or smaller cutoff than 1E-12; it depends on the operators and what you need them for.

By the way, I would strongly recommend against using the float type in C++; use double instead or Real which is defined in ITensor to be double.

Also you don't need to write MPOt. The type MPO is an alias for this. Even better you can write auto H = MPO(H_ampo); and the compiler will deduce that the object on the left has the same type as on the right. Similarly you don't need to provide template parameters to call nmultMPO; the compiler will deduce the template parameters of functions for you whenever they are unambiguous.

Miles

commented by (49k points)
Here is an updated version of your code with the recommended C++ style:

int main()
{
int  L   = 8;
Real hz  = 2.0;
Real cutoff = 1E-12;

auto sites = SpinHalf(L);
auto H_ampo = AutoMPO(sites);

for(int j = 1; j <=  L; j++)
{
H_ampo += hz, "Sz",j;
}

auto H  = MPO(H_ampo);
auto Hn = MPO(H_ampo);
auto res = MPO{};

Print(H.A(2).inds());

nmultMPO(H, Hn, res,{"Cutoff",cutoff});

return 0;
}