Proper indexing for multiple (>1) redundant sites.

diff --git a/ham_ss.cpp b/ham_ss.cpp
index b0793fb..573b355 100644 (file)
--- a/ham_ss.cpp
+++ b/ham_ss.cpp
@@ -310,7 +310,8 @@ int PottsHamiltonian::insertReferenceStates()
     //    return 0;
 
     Coupling newJ;
-    newJ.reserve(J.size() + numX*size);
+    newJ.reserve(J.size() + numX*size - numX*(numX-1)/2);
+    //newJ.reserve(size*(size+1)/2);
 
     Coupling::iterator Jbegin = J.begin();
     //Coupling::iterator end = J.begin();
@@ -325,6 +326,8 @@ int PottsHamiltonian::insertReferenceStates()
     unsigned srcBeg = 0;
     unsigned srcEnd = 0;
     unsigned Xdone = 0;
+    unsigned sectionBeg = 0;
+    bool copyToEnd = true;
     //unsigned dstBeg = 0;
     //unsigned dstEnd = 0;
     // iterate over regions in the destination vector
@@ -332,12 +335,15 @@ int PottsHamiltonian::insertReferenceStates()
         unsigned secLen = size - sec;   // dst
         //dstBeg = index(sec,sec+1,size) - size + sec;
         //equivalently, srcBeg = srcEnd; ... but this may prevent loop unrolling?
-        srcBeg = index(sec,sec+1,sizeOld) - sizeOld + sec;
-        srcEnd = srcBeg;
+        //srcBeg = index(sec,sec+1,sizeOld) - sizeOld + sec;
+        //srcEnd = srcBeg;
+        srcBeg = srcEnd;    // need to start at end of last one
+        sectionBeg = srcBeg;
+        copyToEnd = true;
         //dstBeg = index(sec,sec+1,size) - secLen;
         for (size_t i = Xdone; i < numX; ++i) {
             if (sec <= Xidx[i]) {
-                srcEnd += Xidx[i] - i - sec;
+                srcEnd = sectionBeg + Xidx[i] + Xdone - i - sec;
                 // copy
                 newJ.insert(newJ.end(), Jbegin+srcBeg, Jbegin+srcEnd);
                 newJ.push_back(std::vector<double>());
@@ -349,6 +355,8 @@ int PottsHamiltonian::insertReferenceStates()
                 // insert a whole bunch of nothings
                 newJ.insert(newJ.end(), secLen, std::vector<double>());
                 ++Xdone;    // skip this index in subsequent blocks
+                copyToEnd = false;
+                break;
                 //for (size_t k = 0; k < secLen; ++k)
                 //    newJ.push_back(std::vector<double>());
             }
@@ -356,11 +364,13 @@ int PottsHamiltonian::insertReferenceStates()
             //srcBeg = srcEnd;
             // TODO: consider using a deque and Xidx.pop_front()
         }
-        //end = J.begin() + (size-sec)*(size-sec-1)/2;  // or something
-        //dstEnd = index(sec,sec+1,size);
-        srcEnd = index(sec,sec+1,sizeOld);
-        // copy from begin to end of section
-        newJ.insert(newJ.end(), Jbegin+srcBeg, Jbegin+srcEnd);
+        if (copyToEnd) {
+            //end = J.begin() + (size-sec)*(size-sec-1)/2;  // or something
+            //dstEnd = index(sec,sec+1,size);
+            srcEnd = index(sec-Xdone,sec-Xdone+1,sizeOld);
+            // copy from begin to end of section
+            newJ.insert(newJ.end(), Jbegin+srcBeg, Jbegin+srcEnd);
+        }
     }
 
     int inserted = 0;