Try lanczos to set up hdcr subspace

lib/Config.h.in

@@ -140,9 +140,6 @@
 /* Define to the one symbol short name of this package. */
 #undef PACKAGE_TARNAME
 
-/* Define to the home page for this package. */
-#undef PACKAGE_URL
-
 /* Define to the version of this package. */
 #undef PACKAGE_VERSION
 

lib/algorithms/CoarsenedMatrix.h

@@ -147,6 +147,36 @@ namespace Grid {
       }
       Orthogonalise();
     }
+
+    virtual void CreateSubspaceLanczos(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) 
+    {
+      // Run a Lanczos with sloppy convergence
+	const int Nstop = nn;
+	const int Nk = nn+10;
+	const int Np = nn+10;
+	const int Nm = Nk+Np;
+	const int MaxIt= 10000;
+	RealD resid = 1.0e-5;
+
+	Chebyshev<FineField> Cheb(0.2,5.,11);
+	ImplicitlyRestartedLanczos<FineField> IRL(hermop,Cheb,Nstop,Nk,Nm,resid,MaxIt);
+
+	FineField noise(FineGrid); gaussian(RNG,noise);
+	std::vector<RealD>     eval(Nm);
+	std::vector<FineField> evec(Nm,FineGrid);
+
+	int Nconv;
+	IRL.calc(eval,evec,
+		 noise,
+		 Nconv);
+
+    	// pull back nn vectors
+	for(int b=0;b<nn;b++){
+	  subspace[b]   = evec[b];
+	}
+	Orthogonalise();
+    }
+
     virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
 
       RealD scale;

tests/Test_dwf_hdcr.cc

@@ -561,6 +561,7 @@ int main (int argc, char ** argv)
   int nb=nbasis/2;
   std::cout<<GridLogMessage << " nbasis/2 = "<<nb<<std::endl;
   Aggregates.CreateSubspace(RNG5,HermDefOp,nb);
+    //  Aggregates.CreateSubspaceLanczos(RNG5,HermDefOp,nb);
   for(int n=0;n<nb;n++){
     G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
     std::cout<<GridLogMessage<<n<<" subspace "<<norm2(Aggregates.subspace[n+nb])<<" "<<norm2(Aggregates.subspace[n]) <<std::endl;