checkpoint documentation and bug fixes

TestSuite/inputs/input23.inp

@@ -8,7 +8,7 @@ Connectors 1 1.0
 SolverOptions=hasQuantumNumbers,nowft,nosu2,,
 Version=55460ffd4e0d2587072e9595d7d4ed211c66e83e
 OutputFile=data23.txt
-InfiniteLoopKeptStates=64
-FiniteLoops 4  7 100 0 -7 100 0 -7 100 0 7 100 0 
+InfiniteLoopKeptStates=60
+FiniteLoops 3  7 100 0 -5 100 0 -2 100 0 
 TargetQuantumNumbers 2 0.5 0.5
 

TestSuite/inputs/input24.inp

@@ -8,7 +8,9 @@ Connectors 1 1.0
 SolverOptions=hasQuantumNumbers,nowft,nosu2,,checkpoint
 Version=55460ffd4e0d2587072e9595d7d4ed211c66e83e
 OutputFile=data24.txt
-InfiniteLoopKeptStates=64
-FiniteLoops 4  7 200 0 -7 200 0 -7 200 0 7 200 0 
+InfiniteLoopKeptStates=60
+FiniteLoops 2   -7 200 0 7 200 0 
 TargetQuantumNumbers 2 0.5 0.5
-checkpoint 1 3 oracles/dataValidated23.txt   
+CheckpointEnabled=1
+CheckpointFilename=data23.txt
+

TestSuite/inputs/processing12.txt

@@ -0,0 +1,7 @@
+energy
+gprof
+#observables
+#C
+#N
+#Sz
+dmrg

doc/manual.tex

@@ -849,7 +849,47 @@ \subsection{Wave Function Transformation}
 the old system-environment permutation $P^{SE}$, and the new system permutation $P^S$.
 A sum should be assumed for all indices except $\eta_p$.
 \subsection{SU(2) Reduced Operators}
+
 \subsection{Checkpointing}
+Let's say you want to first run 3 moves on a 16-site lattice:
+\verb=7 100 0 -5 100 0 -2 100 0=, and then continue it with 2 more
+\verb=-7 200 0 7 200 0=. Then one needs to set up input files as in 
+\verb=TestSuite/inputs/input23.inp=
+and
+\verb=TestSuite/inputs/input24.inp=.
+When you run \verb=configure.pl= make sure that you select \verb=DiskStacks= as
+the stacks will need to be saved to disk for resuming.
+Then, compile and do the first run:
+\begin{verbatim}
+cd src/
+perl configure.pl < ../TestSuite/inputs/model23.spec
+(or say ModelHeisenberg and DiskStack)
+make
+./dmrg ../TestSuite/inputs/input23.inp
+\end{verbatim}
+and now, one should see files
+\verb=SystemStackdata.txt= and \verb=EnvironStackdata.txt= as well
+as a fatter than usual data23.txt. You are ready for the continuation of this run now with:
+\begin{verbatim}
+./dmrg ../TestSuite/inputs/input24.inp
+\end{verbatim}
+Note that the continued run's input (\verb=input24.inp=) has the option checkpoint, the CheckpointEnabled tag
+and the CheckpointFilename tag.
+Continued (raw) results will be in file \verb=data24= as usual.
+
+Note the following caveats and ``todos'':
+\begin{itemize}
+\item Run will be (slightly?) slower due to the use of DiskStacks.
+\item There's no check (yet) of finite loops for consistency while checkpoint is in use.
+Therefore, make sure the second run is starting at a point on the lattice where the
+previous to-be-continued run left off.
+\item Bases are stored ("checkpointed") into disk and accumulate without being purged.
+\item Stacks are written to disk and never purged.
+\end{itemize} 
+The last two caveats have the positive effect of allowing to resume at any point of a previous run.
+On the other hand, a option to purge and leave only the most recent resume checkpoint would be useful, as it
+would reduce disk usage, and loading times later.
+
 \subsection{Distributed Parallelization}
 \subsection{Shared-memory Parallelization}
 

src/Engine/Checkpoint.h

@@ -77,6 +77,8 @@ DISCLOSED WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
 /*! \file Checkpoint.h
  *
  *  checkpointing functions
+ *  this class also owns the stacks since they 
+ *  are so related to checkpointing
  */
 #ifndef CHECKPOINT_H
 #define CHECKPOINT_H
@@ -112,7 +114,19 @@ namespace Dmrg {
 			void load(BasisType &pS,BasisType &pE)
 			{
 				typename IoType::In ioTmp(parameters_.checkpoint.filename);
-				BasisType pS1(ioTmp,"#CHKPOINTSYSTEM",parameters_.checkpoint.index);
+				size_t loop = ioTmp.count("#NAME=#CHKPOINTSYSTEM");
+				if (loop<1) {
+					std::cerr<<"There are no resumable loops in file "<<parameters_.checkpoint.filename<<"\n";
+					throw std::runtime_error("Checkpoint::load(...)\n");
+				}
+				loop--;
+				if (loop<parameters_.checkpoint.index) {
+					std::cerr<<"There are "<<loop<<" resumable checkpoints\n";
+					std::cerr<<"But you requested to go back "<<parameters_.checkpoint.index<<"\n";
+					throw std::runtime_error("Checkpoint::load(...)\n");
+				}
+				loop -= parameters_.checkpoint.index;
+				BasisType pS1(ioTmp,"#CHKPOINTSYSTEM",loop);
 				pS=pS1;
 				BasisType pE1(ioTmp,"#CHKPOINTENVIRON");
 				pE=pE1;

src/Engine/Diagonalization.h

@@ -137,7 +137,7 @@ namespace Dmrg {
 
 
 			std::ostringstream msg;
-			msg<<"Setting up ham. pse.size="<<pSE.size();
+			msg<<"Setting up Hamiltonian basis of size="<<pSE.size();
 			progress_.printline(msg,std::cout);
 
 			TargetVectorType tmpVec;

src/Engine/DmrgSolver.h

@@ -219,7 +219,6 @@ namespace Dmrg {
 				infiniteDmrgLoop(S,X,Y,E,pS,pE,psi);
 			}
 
-			stepCurrent_ = pS.block().size()-1;
 			finiteDmrgLoops(S,E,pS,pE,psi);
 
 			std::ostringstream msg2;
@@ -300,8 +299,29 @@ namespace Dmrg {
        					TargettingType& psi)
 		{
 			useReflection_=false; // disable reflection symmetry for finite loop if it was enabled:
-
 			if (parameters_.options.find("nofiniteloops")!=std::string::npos) return;
+			if (parameters_.finiteLoop.size()==0)
+				throw std::runtime_error("finiteDmrgLoops(...): there are no finite loops! (and nofiniteloops is not set)\n");
+
+			// set initial site to add to either system or environment:
+			// this is a bit tricky and has been a source of endless bugs
+			// basically we have pS on the left and pE on the right, 
+			// and we need to determine which site is to be added
+			// let us set the initial direction first:
+			size_t direction = EXPAND_SYSTEM;
+			if (parameters_.finiteLoop[0].stepLength<0) direction=EXPAND_ENVIRON;
+			// all right, now we can get the actual site to add:
+
+			std::vector<size_t> siteToAdd(1,pE.block()[0]); // left-most site of pE
+			if (direction==EXPAND_ENVIRON) {
+				siteToAdd[0] = pS.block()[pS.block().size()-1]; // right-most site of pS
+			}
+			// now stepCurrent_ is such that sitesIndices_[stepCurrent_] = siteToAdd
+			// so:
+			int sc = utils::isInVector(sitesIndices_,siteToAdd);
+			if (sc<0) throw std::runtime_error("finiteDmrgLoops(...): internal error: siteIndices_\n");
+			stepCurrent_ = sc; // phew!!, that's all folks, now bugs, go away!!
+
 			for (size_t i=0;i<parameters_.finiteLoop.size();i++)  {
 				std::ostringstream msg;
 				msg<<"Finite loop number "<<i;

src/Engine/IoSimple.h

@@ -287,7 +287,7 @@ namespace Dmrg {
 					return sc;
 				}
 
-				std::pair<std::string,size_t> advance(std::string const &s,int level=0)
+				std::pair<std::string,size_t> advance(std::string const &s,int level=0,bool beQuiet=false)
 				{
 					std::string temp="NOTFOUND";
 					std::string tempSaved="NOTFOUND";
@@ -311,14 +311,34 @@ namespace Dmrg {
 					}
 					//std::cerr<<"count="<<c<<"\n";
 					if (!found && tempSaved=="NOTFOUND") {
-						std::cerr<<"Not found "<<s<<" in file "<<filename_;
-						std::cerr<<" level="<<level<<" counter="<<counter<<"\n";
+						if (!beQuiet) {
+							std::cerr<<"Not found "<<s<<" in file "<<filename_;
+							std::cerr<<" level="<<level<<" counter="<<counter<<"\n";
+						}
 						throw std::runtime_error("IoSimple::In::read()\n");
 					}
 					//std::cerr<<"------------\n";
 					return std::pair<std::string,size_t>(tempSaved,counter);
 				}
 
+				size_t count(const std::string& s)
+				{
+					size_t i = 0;
+					while(i<1000) {
+						try {
+							advance(s,0,true);
+							i++;
+						} catch (std::exception& e) {
+							rewind();
+							return i;
+						}
+					}
+					std::string ss = "IoSimple::count(...): too many "
+						+s+" in file "+filename_+"\n";
+					throw std::runtime_error(s.c_str());
+
+				}
+
 				template<typename T>
 				void read(std::vector<std::pair<T,T> > &x,std::string const &s,int level=0)
 				{

src/Engine/ParametersDmrgSolver.h

@@ -156,7 +156,12 @@ namespace Dmrg {
 		void load(IoInputType& io)
 		{
 			io.readline(enabled,"CheckpointEnabled=");
-			io.readline(index,"CheckpointIndex=");
+			try {
+				io.readline(index,"CheckpointIndex=");
+			} catch (std::exception& e) {
+				io.rewind();
+				index=0;
+			}
 			io.readline(filename,"CheckpointFilename=");
 		}
 	};

src/Engine/WaveFunctionTransformation.h

@@ -359,7 +359,7 @@ namespace Dmrg {
 					dmrgWave_.ws=wsStack_.top();
 					wsStack_.pop();
 				} else {
-					std::cerr<<"PUSHING STACK ERROR S\n";
+					//std::cerr<<"PUSHING STACK ERROR S\n";
 					throw std::runtime_error("System Stack is empty\n");
 				}
 			}
@@ -370,7 +370,7 @@ namespace Dmrg {
 					weStack_.pop();
 					//std::cerr<<"CHANGED We taken from stack\n";
 				} else {
-					std::cerr<<"PUSHING STACK ERROR E\n";
+					//std::cerr<<"PUSHING STACK ERROR E\n";
 					throw std::runtime_error("Environ Stack is empty\n");
 				}
 			}

src/Models/HubbardOneBand/ModelHubbard.h

@@ -141,7 +141,7 @@ namespace Dmrg {
 
 		size_t orbitals() const { return NUMBER_OF_ORBITALS; }
 
-		size_t hilbertSize() const { return pow(2,2*NUMBER_OF_ORBITALS); }
+		size_t hilbertSize() const { return (size_t)pow(2,2*NUMBER_OF_ORBITALS); }
 
 		//! find creation operator matrices for (i,sigma) in the natural basis, find quantum numbers and number of electrons
 		//! for each state in the basis

src/configure.pl

@@ -685,12 +685,14 @@ sub createObserverDriver
 	
 	MyConcurrency concurrency(argc,argv);
 	
-	if (argc!=3 && concurrency.name()=="serial") {
-		std::string s = "The observer driver takes two  arguments now: \\n";
+	if (argc<2) {
+		std::string s = "The observer driver takes at least one  argumentw: \\n";
 		s = s + "(i) the name of the input file";
-		s = s+ " and (ii) a comma-separated list of options of what to compute\\n";
+		s = s+ " and, optionally, (ii) a comma-separated list of options of what to compute\\n";
 		throw std::runtime_error(s);
 	}
+	std::string options = "cc,nn,szsz";
+	if (argc>2) options = argv[2];
 	
 	//Setup the Geometry
 	typedef Geometry<RealType> GeometryType;
@@ -712,11 +714,11 @@ sub createObserverDriver
 	if (hasTimeEvolution)
 		mainLoop<ParametersModelType,GeometryType,MyConcurrency,IoInputType,$modelName,
 			ModelHelperLocal,InternalProductOnTheFly,VectorWithOffsets,
-			$targetting,MySparseMatrixReal>(mp,geometry,hasTimeEvolution,concurrency,io,dmrgSolverParams.filename,argv[2]);
+			$targetting,MySparseMatrixReal>(mp,geometry,hasTimeEvolution,concurrency,io,dmrgSolverParams.filename,options);
 	else
 		mainLoop<ParametersModelType,GeometryType,MyConcurrency,IoInputType,$modelName,
 			ModelHelperLocal,InternalProductOnTheFly,VectorWithOffset,
-			$targetting,MySparseMatrixReal>(mp,geometry,hasTimeEvolution,concurrency,io,dmrgSolverParams.filename,argv[2]);
+			$targetting,MySparseMatrixReal>(mp,geometry,hasTimeEvolution,concurrency,io,dmrgSolverParams.filename,options);
 } // main
 
 EOF