Problem in Lanczos with Hilbert space greater than 2^32 (Sz=0)

[llviteritti]

Hi, using the xcode version #49f8fe7, @riccardo-rende and I tried to diagonalize an Heisenberg chain using this code :

#include <xdiag/all.hpp>

using namespace xdiag;

int main() try {

  int n_sites = 36;
  int nup = n_sites / 2;
  Spinhalf block(n_sites, nup);

  // Define the nearest-neighbor Heisenberg model
  BondList bonds;
  for (int i = 0; i < n_sites; ++i) {
    bonds += Bond("HB", "J", {i, (i + 1) % n_sites});
  }
  bonds["J"] = 1.0;

  set_verbosity(2);                  // set verbosity for monitoring progress
  double e0 = eigval0(bonds, block); // compute ground state energy

  Log("Ground state energy: {:.12f}", e0);

} catch (Error e) {
  error_trace(e);
}

We checked that up to L=32 we get the correct results. However, for L=34 we obtain the following error:

Warning: Tmatrix zero dimensional in eig0_cplx
Ground state energy: nan

and for L=36 the code runs but the energy is wrong, namely it goes below the exact one :

Exact energy per site : -0.44378432839067183
Lanczos energy per site after 64 iterations : -0.4737482323668045

The same happens for the square lattice, the energy for the 4x4 is correct, instead for the 6x6 goes below the exact one.

[awietek]

Dear @llviteritti, thank you very much. This is actually the first issue someone is raising so it is good to know first people are using the code. This is a bug, that will need to be fixed. I will look into it and keep you updated.

May I ask, which linear algebra backend you are using, e.g. LAPACK, MKL?

There exist LAPACK interfaces that only allow for sizes smaller than 2^32. If you are using one of these, then you would need to link to a 64-bit interface.

[awietek]

I just tried to reproduce the error without success. On my side, the correct results are computed. I was using the Intel MKL linear algebra backend, which has a 64 bit interface.

It would be good to know if you are using a 32bit Lapack implementation. I think XDiag should warn or even throw an error if this is indeed the case.

[llviteritti]

Dear @awietek, thank you for your quick reply. We have now installed XDiag with the Intel MKL linear algebra backend and the simulations are running for L=36. We will update you as soon as the Lanczos iterations are complete.

[gcarleo]

thanks Alex ! this was a subtle issue 😄

[awietek]

Ah, cool, that's good to hear. Do you get the right energy now? Still, something like this should not happen. An error should be thrown at compile time if one tries to link to 32 bit, or an error at runtime if it's linked to 32 bit but the vector is too big. I'd like to understand this better.

[llviteritti]

Yes, now we get the right energies! Thank you again.

[awietek]

Good to know. I am keeping this issue open, since there should be a solution to when one is trying to link to a 32 bit Lapack implementation.

[awietek]

@llviteritti I now added a security check which checks that the block dimension is not too big if a 32 bit Lapack backend is used. Could you let me know if with your initial setup this error gets thrown correctly on your machine? If so I will close the issue.

[llviteritti]

Using the latest version of the commit #2640a8d and BLAS with 32-bit interface, now gives the error:

Trying to create a block whose dimension is too large for the backend BLAS routines. Your backend BLAS routines have a 32 bit interface, and a vector beyond a dimension of 2^31 cannot be represented. If you want to still perform this calculation, you will need to compile XDiag with a 64 bit backend, for example the ILP64 interface of IntelMKL.

So it should be fine now. Thank you very much for the support!

[awietek]

Thanks for checking. Yes, this is exactly what should be shown. I'll be closing this then