Segmentation fault when trying to read in site and wavefunctions

Question

Hello,

I'm attempting to read in a siteset and wavefunctions from a previous run, however I am running into some issues that I don't quite understand. From the documentation, I believe the code should look something like

int main()
{
    int N = 3;
    int dim = 10;

    auto sites = Exciton(N,{"on_site_dim", dim});
    readFromFile(format("sites_file_NB_%d",N),sites);
    auto psi = MPS(sites);
    readFromFile(format("psi_file_NB_%d",N),psi);

    auto ampo = AutoMPO(sites);

    ...

    auto H = toMPO(ampo);
}

where I define my Hamiltonian in the ellipses. However, the last line throws a segmentation fault. I have confirmed that if I don't try and read in a site file or wavefunction file, everything runs smoothly. Any ideas on what the problem might be? Thanks again for all of your help, Miles and Matt!

Comments

Hi, thanks for the question. It may have something to do with how you defined the Exciton class. (Unfortunately because of some fundamental issues it was hard for us to make reading/writing of custom site sets a complete no-brainer.) Could you post the code for that and/or tell me if you defined a read and write method at all? Or did you use one of the site set templates like BasicSiteSet or MixedSiteSet?

---

Miles,

No I definitely didn't define a read and write method. Here is the code for how I defined Exciton - I used BasicSiteSet:

namespace itensor {

class ExcitonSite;

using Exciton = BasicSiteSet<ExcitonSite>;

class ExcitonSite
    {
    Index s;
    int dim;
    public:

    ExcitonSite() { }

    ExcitonSite(Index I) : s(I) { }

    ExcitonSite(Args const& args = Args::global())
    {
        dim = args.getInt("on_site_dim", 5);
        auto ts = TagSet("Site,Exciton");
        int n = args.getInt("SiteNumber");
        ts.addTags("n="+str(n));
        s = Index(2 * dim + 1, ts);
    }
    
    ExcitonSite(int n, Args const& args = Args::global())
    {
        *this = ExcitonSite({args,"SiteNumber=",n});
    }

    Index
    index() const { return s; }

    IndexVal
    state(std::string const& state)
    {
        int j = -dim;
        while (j <= dim)
        {
            j++;
            if (state == str(j))
            {
                return s(j + 1 + dim);
                break;
            }
        }
        return IndexVal{};
    }

    // operators defined here

    };
} //namespace itensor

---

Ok thanks. As long as you used the BasicSiteSet template, which you did, that will define read and write for you.

So at this point I would need more information to diagnose the bug. Have you run the code in debug mode? Often that gives better error messages, and hopefully a more informative one than a segfault.

---

Here is the error I get in debug mode:

inds = 22 21
From line 495, file /home/dkweiss/danny/software/v3/itensor/itensor/itdata/task_types.h

Out of range: IndexVal at position 1 has val 22 > (21|id=787|n=1,Site,Exciton)

Out of range: IndexVal at position 1 has val 22 > (21|id=787|n=1,Site,Exciton)
Aborted (core dumped)

This would seem to imply that I've changed dim (right?), but I've kept it constant from each run to the next.

---

Hi, that's very helpful although I think it's still unclear what the real cause of the bug is. I don't think the readFromFile function would make an IndexVal at any point, so it could be code later down that's causing the error. I think a good thing to do would be to print out the site set just after you read it back from disk to make sure it looks ok (like the Index sizes are right).

Also it would be really helpful to run the code in a debugger to see what line of your driver code is the last one that runs before you get the error. Likely it's toMPO which is making operators in the process of MPO construction.

---

Indeed it looks ok after reading it in - here is the output of the code

    auto sites = Exciton(N,{"on_site_dim", dim});
    Print(sites);
    readFromFile(format("sites_file_NB_%d",N),sites);
    Print(sites);

sites =
SiteSet:
site 1 = (21|id=524|n=1,Site,Exciton)
site 2 = (21|id=576|n=2,Site,Exciton)
site 3 = (21|id=420|n=3,Site,Exciton)

sites =
SiteSet:
site 1 = (21|id=787|n=1,Site,Exciton)
site 2 = (21|id=131|n=2,Site,Exciton)
site 3 = (21|id=147|n=3,Site,Exciton)

And yes, it is the toMPO line that is the last line run.

---

Ok great. Then one culprit could be your “op” function which makes your custom operators. Is any of that code trying to access a tensor element past the size of the site indices?

---

Gotcha - I don't think so, I just checked it carefully. Though if that were the case, I would think I should've run into this issue before now? Here is the code though, for the sake of completeness. Thanks again Miles for helping with this!

ITensor
    op(std::string const& opname,
       Args const& args) const
        {
        auto sP = prime(s);
        auto Op = ITensor(dag(s),sP);

        if (opname == "n")
        {
            for (int j = -dim; j<= dim; ++j)
            {
                Op.set(s(j + 1 + dim), sP(j + 1 + dim), (float)j);
            }
        }
        else
        if (opname == "nsq")
        {
            for (int j = -dim; j<= dim; ++j)
            {
                Op.set(s(j + 1 + dim), sP(j + 1 + dim), (float)(j * j));
            }
        }
        else
        if (opname == "gp")
        {
            for (int j = -dim; j<= dim - 1; ++j)
            {
                Op.set(s(j + 1 + dim), sP(j + 2 + dim), +1.0);
            }
        }
        else
        if (opname == "gm")
        {
            for (int j = -dim + 1; j<= dim; ++j)
            {
                Op.set(s(j + 1 + dim), sP(j + dim), +1.0);
            }
        }
        else
        if (opname == "i_o")
        {
            for (int j=-dim; j<= dim; ++j)
            {
                Op.set(s(j + 1 + dim), sP(j + 1 + dim), +1.0);
            }
        }
        else
            {
            Error("Operator \"" + opname + "\" name not recognized");
            }

        return Op;
        }

---

One question I have is: how does the variable "dim" get defined in this function? I don't see it being defined or passed in. So is it a global variable?

---

Yes it is a global variable - see my first comment where it is defined underneath

Index s;

---

Hi, I see. There's a strong chance this is causing problems. I would recommend redesigning op to analyze the tensor index (calling dim(s) on it) to figure out what dim is, not by using a global variable.

Also you could put a line Print(dim); inside the op function to check whether dim is being set to 22 or 21 or the correct value or not.

---

This helped me better diagnose the problem, though I'm still not sure what the solution would be. The problem is that when I specify dim, typically dim=10, that creates an index with dimension 2*dim+1=21. Then, when that siteset is read back in, the code believes it's reading in an index with dim=21, and proceeds to try to access entries in s from 1 to 2*21+1.

So, long story short, my code might need some serious surgery in order to be able to use the read and write functions?

---

Yes, that sounds like the problem. Why does the code think that when it's reading the indices back in, that it should interpret the dimension of the index as the number "n" that goes into the 2*n+1 formula? I don't think that's something ITensor is doing itself, but it likely a logical flaw in your code somewhere.

I think a much better design that would avoid this is to not have dim be defined in your ExcitonSite class itself. Rather, it just holds an index of some size, and knows how to extract the number you're calling dim from the size of this index (as  (dim(i)-1)/2 ). This would be used in the op function to set the for loop bounds, for example.

How does that sound?

---

That solution works perfectly. As usual, thanks a million Miles.

---

Glad to hear it!

Answer 1

Question is answered in above comments