Minor changes fixes

lib/algorithms/iterative/ConjugateGradient.h

@@ -52,6 +52,7 @@ class ConjugateGradient : public OperatorFunction<Field> {
         MaxIterations(maxit),
         ErrorOnNoConverge(err_on_no_conv){};
 
+
   void operator()(LinearOperatorBase<Field> &Linop, const Field &src,
                   Field &psi) {
     psi.checkerboard = src.checkerboard;

lib/algorithms/iterative/ConjugateGradientMixedPrec.h

@@ -60,6 +60,7 @@ namespace Grid {
     }
   
     void operator() (const FieldD &src_d_in, FieldD &sol_d){
+
       TotalInnerIterations = 0;
 	
       GridStopWatch TotalTimer;

lib/algorithms/iterative/SimpleLanczos.h

@@ -51,12 +51,16 @@ void LAPACK_dstegr(char *jobz, char *range, int *n, double *d, double *e,
 // eliminate temorary vector in calc()
 #define MEM_SAVE
 
-namespace Grid {
+namespace Grid
+{
 
-struct Bisection {
+  struct Bisection
+  {
 
 #if 0
-static void get_eig2(int row_num,std::vector<RealD> &ALPHA,std::vector<RealD> &BETA, std::vector<RealD> & eig)
+    static void get_eig2 (int row_num, std::vector < RealD > &ALPHA,
+			  std::vector < RealD > &BETA,
+			  std::vector < RealD > &eig)
     {
       int i, j;
         std::vector < RealD > evec1 (row_num + 3);
@@ -64,7 +68,8 @@ static void get_eig2(int row_num,std::vector<RealD> &ALPHA,std::vector<RealD> &B
       RealD eps2;
         ALPHA[1] = 0.;
         BETHA[1] = 0.;
-  for(i=0;i<row_num-1;i++) {
+      for (i = 0; i < row_num - 1; i++)
+	{
 	  ALPHA[i + 1] = A[i * (row_num + 1)].real ();
 	  BETHA[i + 2] = A[i * (row_num + 1) + 1].real ();
 	}
@@ -76,17 +81,22 @@ static void get_eig2(int row_num,std::vector<RealD> &ALPHA,std::vector<RealD> &B
       int begin = 1;
       int end = row_num;
       int swapped = 1;
-  while(swapped) {
+      while (swapped)
+	{
 	  swapped = 0;
-    for(i=begin;i<end;i++){
+	  for (i = begin; i < end; i++)
-      if(mag(evec2[i])>mag(evec2[i+1]))	{
+	    {
+	      if (mag (evec2[i]) > mag (evec2[i + 1]))
+		{
 		  swap (evec2 + i, evec2 + i + 1);
 		  swapped = 1;
 		}
 	    }
 	  end--;
-    for(i=end-1;i>=begin;i--){
+	  for (i = end - 1; i >= begin; i--)
-      if(mag(evec2[i])>mag(evec2[i+1]))	{
+	    {
+	      if (mag (evec2[i]) > mag (evec2[i + 1]))
+		{
 		  swap (evec2 + i, evec2 + i + 1);
 		  swapped = 1;
 		}
@@ -94,10 +104,14 @@ static void get_eig2(int row_num,std::vector<RealD> &ALPHA,std::vector<RealD> &B
 	  begin++;
 	}
 
-  for(i=0;i<row_num;i++){
+      for (i = 0; i < row_num; i++)
-    for(j=0;j<row_num;j++) {
+	{
-      if(i==j) H[i*row_num+j]=evec2[i+1];
+	  for (j = 0; j < row_num; j++)
-      else H[i*row_num+j]=0.;
+	    {
+	      if (i == j)
+		H[i * row_num + j] = evec2[i + 1];
+	      else
+		H[i * row_num + j] = 0.;
 	    }
 	}
     }
@@ -109,9 +123,7 @@ static void bisec(std::vector<RealD> &c,
 		       int m1,
 		       int m2,
 		       RealD eps1,
-		  RealD relfeh,
+		       RealD relfeh, std::vector < RealD > &x, RealD & eps2)
-		  std::vector<RealD> &x,
-		  RealD &eps2)
     {
       std::vector < RealD > wu (n + 2);
 
@@ -121,53 +133,75 @@ static void bisec(std::vector<RealD> &c,
       b[1] = 0.0;
       xmin = c[n] - fabs (b[n]);
       xmax = c[n] + fabs (b[n]);
-  for(i=1;i<n;i++){
+      for (i = 1; i < n; i++)
+	{
 	  h = fabs (b[i]) + fabs (b[i + 1]);
-    if(c[i]+h>xmax) xmax= c[i]+h;
+	  if (c[i] + h > xmax)
-    if(c[i]-h<xmin) xmin= c[i]-h;
+	    xmax = c[i] + h;
+	  if (c[i] - h < xmin)
+	    xmin = c[i] - h;
 	}
       xmax *= 2.;
 
       eps2 = relfeh * ((xmin + xmax) > 0.0 ? xmax : -xmin);
-  if(eps1<=0.0) eps1=eps2;
+      if (eps1 <= 0.0)
+	eps1 = eps2;
       eps2 = 0.5 * eps1 + 7.0 * (eps2);
       x0 = xmax;
-  for(i=m1;i<=m2;i++){
+      for (i = m1; i <= m2; i++)
+	{
 	  x[i] = xmax;
 	  wu[i] = xmin;
 	}
 
-  for(k=m2;k>=m1;k--){
+      for (k = m2; k >= m1; k--)
+	{
 	  xu = xmin;
 	  i = k;
-    do{
+	  do
-      if(xu<wu[i]){
+	    {
+	      if (xu < wu[i])
+		{
 		  xu = wu[i];
 		  i = m1 - 1;
 		}
 	      i--;
-    }while(i>=m1);
+	    }
-    if(x0>x[k]) x0=x[k];
+	  while (i >= m1);
-    while((x0-xu)>2*relfeh*(fabs(xu)+fabs(x0))+eps1){
+	  if (x0 > x[k])
+	    x0 = x[k];
+	  while ((x0 - xu) > 2 * relfeh * (fabs (xu) + fabs (x0)) + eps1)
+	    {
 	      x1 = (xu + x0) / 2;
 
 	      a = 0;
 	      q = 1.0;
-      for(i=1;i<=n;i++){
+	      for (i = 1; i <= n; i++)
-	q=c[i]-x1-((q!=0.0)? b[i]*b[i]/q:fabs(b[i])/relfeh);
+		{
-	if(q<0) a++;
+		  q =
+		    c[i] - x1 -
+		    ((q != 0.0) ? b[i] * b[i] / q : fabs (b[i]) / relfeh);
+		  if (q < 0)
+		    a++;
 		}
 //      printf("x1=%0.14e a=%d\n",x1,a);
-      if(a<k){
+	      if (a < k)
-	if(a<m1){
+		{
+		  if (a < m1)
+		    {
 		      xu = x1;
 		      wu[m1] = x1;
-	}else {
+		    }
+		  else
+		    {
 		      xu = x1;
 		      wu[a + 1] = x1;
-	  if(x[a]>x1) x[a]=x1;
+		      if (x[a] > x1)
+			x[a] = x1;
 		    }
-      }else x0=x1;
+		}
+	      else
+		x0 = x1;
 	    }
 	  printf ("x0=%0.14e xu=%0.14e k=%d\n", x0, xu, k);
 	  x[k] = (x0 + xu) / 2;
@@ -180,8 +214,8 @@ static void bisec(std::vector<RealD> &c,
 /////////////////////////////////////////////////////////////
 
 
-template<class Field> 
+  template < class Field > class SimpleLanczos
-    class SimpleLanczos {
+  {
 
     const RealD small = 1.0e-16;
   public:
@@ -209,50 +243,60 @@ public:
     /////////////////////////
     // Constructor
     /////////////////////////
-    void init(void){};
+    void init (void)
-    void Abort(int ff, DenseVector<RealD> &evals,  DenseVector<DenseVector<RealD> > &evecs);
+    {
+    };
+    void Abort (int ff, DenseVector < RealD > &evals,
+		DenseVector < DenseVector < RealD > >&evecs);
 
-    SimpleLanczos(
+    SimpleLanczos (LinearOperatorBase < Field > &Linop,	// op
-				LinearOperatorBase<Field> &Linop, // op
 		   OperatorFunction < Field > &poly,	// polynmial
 		   int _Nstop,	// sought vecs
 		   int _Nk,	// sought vecs
 		   int _Nm,	// spare vecs
 		   RealD _eresid,	// resid in lmdue deficit 
 		   int _Niter):	// Max iterations
+     
       _Linop (Linop),
       _poly (poly),
-      Nstop(_Nstop),
+      Nstop (_Nstop), Nk (_Nk), Nm (_Nm), eresid (_eresid), Niter (_Niter)
-      Nk(_Nk),
-      Nm(_Nm),
-      eresid(_eresid),
-      Niter(_Niter)
     {
-      Np = Nm-Nk; assert(Np>0);
+      Np = Nm - Nk;
+      assert (Np > 0);
     };
 
     /////////////////////////
     // Sanity checked this routine (step) against Saad.
     /////////////////////////
-    void RitzMatrix(DenseVector<Field>& evec,int k){
+    void RitzMatrix (DenseVector < Field > &evec, int k)
+    {
 
-      if(1) return;
+      if (1)
+	return;
 
       GridBase *grid = evec[0]._grid;
       Field w (grid);
       std::cout << GridLogMessage << "RitzMatrix " << std::endl;
-      for(int i=0;i<k;i++){
+      for (int i = 0; i < k; i++)
-	_poly(_Linop,evec[i],w);
+	{
+	  _Linop.HermOp (evec[i], w);
+//      _poly(_Linop,evec[i],w);
 	  std::cout << GridLogMessage << "[" << i << "] ";
-	for(int j=0;j<k;j++){
+	  for (int j = 0; j < k; j++)
+	    {
 	      ComplexD in = innerProduct (evec[j], w);
-	  if ( fabs((double)i-j)>1 ) { 
+	      if (fabs ((double) i - j) > 1)
-	    if (abs(in) >1.0e-9 )  { 
+		{
+		  if (abs (in) > 1.0e-9)
+		    {
 		      std::cout << GridLogMessage << "oops" << std::endl;
 		      abort ();
-	    } else 
+		    }
+		  else
 		    std::cout << GridLogMessage << " 0 ";
-	  } else { 
+		}
+	      else
+		{
 		  std::cout << GridLogMessage << " " << in << " ";
 		}
 	    }
@@ -262,17 +306,18 @@ public:
 
     void step (DenseVector < RealD > &lmd,
 	       DenseVector < RealD > &lme,
-	    Field& last,
+	       Field & last, Field & current, Field & next, uint64_t k)
-	    Field& current, 
-	  	Field & next,
-		uint64_t k)
     {
-	  if (lmd.size()<=k) lmd.resize(k+Nm);
+      if (lmd.size () <= k)
-	  if (lme.size()<=k) lme.resize(k+Nm);
+	lmd.resize (k + Nm);
+      if (lme.size () <= k)
+	lme.resize (k + Nm);
 
 
-      _poly(_Linop,current,next );   // 3. wk:=Avk−βkv_{k−1}
+//      _poly(_Linop,current,next );   // 3. wk:=Avk−βkv_{k−1}
-      if(k>0){
+      _Linop.HermOp (current, next);	// 3. wk:=Avk−βkv_{k−1}
+      if (k > 0)
+	{
 	  next -= lme[k - 1] * last;
 	}
 //      std::cout<<GridLogMessage << "<last|next>" << innerProduct(last,next) <<std::endl;
@@ -289,10 +334,14 @@ public:
 
       int interval = Nm / 100 + 1;
       if ((k % interval) == 0)
-	std::cout<<GridLogMessage << k << " : alpha = " << zalph << " beta "<<beta<<std::endl;
+	std::
+	  cout << GridLogMessage << k << " : alpha = " << zalph << " beta " <<
+	  beta << std::endl;
       const RealD tiny = 1.0e-20;
-      if ( beta < tiny ) { 
+      if (beta < tiny)
-	std::cout<<GridLogMessage << " beta is tiny "<<beta<<std::endl;
+	{
+	  std::cout << GridLogMessage << " beta is tiny " << beta << std::
+	    endl;
 	}
       lmd[k] = alph;
       lme[k] = beta;
@@ -303,10 +352,7 @@ public:
 		    DenseVector < RealD > &lme,
 		    int Nk,
 		    int Nm,
-		   DenseVector<RealD>& Qt,
+		    DenseVector < RealD > &Qt, RealD Dsh, int kmin, int kmax)
-		   RealD Dsh, 
-		   int kmin,
-		   int kmax)
     {
       int k = kmin - 1;
       RealD x;
@@ -325,7 +371,8 @@ public:
       x = -s * lme[k + 1];
       lme[k + 1] = c * lme[k + 1];
 
-      for(int i=0; i<Nk; ++i){
+      for (int i = 0; i < Nk; ++i)
+	{
 	  RealD Qtmp1 = Qt[i + Nm * k];
 	  RealD Qtmp2 = Qt[i + Nm * (k + 1)];
 	  Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
@@ -333,7 +380,8 @@ public:
 	}
 
       // Givens transformations
-      for(int k = kmin; k < kmax-1; ++k){
+      for (int k = kmin; k < kmax - 1; ++k)
+	{
 
 	  RealD Fden = 1.0 / hypot (x, lme[k - 1]);
 	  RealD c = lme[k - 1] * Fden;
@@ -348,12 +396,14 @@ public:
 	  lme[k] = c * s * (tmpa1 - tmpa2) + (c * c - s * s) * tmpb;
 	  lme[k - 1] = c * lme[k - 1] - s * x;
 
-	if(k != kmax-2){
+	  if (k != kmax - 2)
+	    {
 	      x = -s * lme[k + 1];
 	      lme[k + 1] = c * lme[k + 1];
 	    }
 
-	for(int i=0; i<Nk; ++i){
+	  for (int i = 0; i < Nk; ++i)
+	    {
 	      RealD Qtmp1 = Qt[i + Nm * k];
 	      RealD Qtmp2 = Qt[i + Nm * (k + 1)];
 	      Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
@@ -368,11 +418,10 @@ public:
 #else
 #define LAPACK_INT long long
 #endif
-    void diagonalize_lapack(DenseVector<RealD>& lmd,
+    void diagonalize_lapack (DenseVector < RealD > &lmd, DenseVector < RealD > &lme, int N1,	// all
-		     DenseVector<RealD>& lme, 
-		     int N1, // all
 			     int N2,	// get
-		     GridBase *grid) {
+			     GridBase * grid)
+    {
       const int size = Nm;
       LAPACK_INT NN = N1;
       double evals_tmp[NN];
@@ -380,13 +429,19 @@ public:
       double EE[NN];
       for (int i = 0; i < NN; i++)
 	for (int j = i - 1; j <= i + 1; j++)
-      if ( j < NN && j >= 0 ) {
+	  if (j < NN && j >= 0)
-        if (i==j) DD[i] = lmd[i];
+	    {
-        if (i==j) evals_tmp[i] = lmd[i];
+	      if (i == j)
-        if (j==(i-1)) EE[j] = lme[j];
+		DD[i] = lmd[i];
+	      if (i == j)
+		evals_tmp[i] = lmd[i];
+	      if (j == (i - 1))
+		EE[j] = lme[j];
 	    }
       LAPACK_INT evals_found;
-  LAPACK_INT lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ;
+      LAPACK_INT lwork =
+	((18 * NN) >
+	 (1 + 4 * NN + NN * NN) ? (18 * NN) : (1 + 4 * NN + NN * NN));
       LAPACK_INT liwork = 3 + NN * 10;
       LAPACK_INT iwork[liwork];
       double work[lwork];
@@ -406,28 +461,31 @@ public:
       int interval = (NN / total) + 1;
       double vl = 0.0, vu = 0.0;
       LAPACK_INT il = interval * node + 1, iu = interval * (node + 1);
-  if (iu > NN)  iu=NN;
+      if (iu > NN)
+	iu = NN;
       double tol = 0.0;
-    if (1) {
+      if (1)
+	{
 	  memset (evals_tmp, 0, sizeof (double) * NN);
-      if ( il <= NN){
+	  if (il <= NN)
+	    {
 	      printf ("total=%d node=%d il=%d iu=%d\n", total, node, il, iu);
 #ifdef USE_MKL
 	      dstegr (&jobz, &range, &NN,
 #else
 	      LAPACK_dstegr (&jobz, &range, &NN,
 #endif
-            (double*)DD, (double*)EE,
+			     (double *) DD, (double *) EE, &vl, &vu, &il, &iu,	// these four are ignored if second parameteris 'A'
-            &vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A'
 			     &tol,	// tolerance
 			     &evals_found, evals_tmp, (double *) NULL, &NN,
-            isuppz,
+			     isuppz, work, &lwork, iwork, &liwork, &info);
-            work, &lwork, iwork, &liwork,
+	      for (int i = iu - 1; i >= il - 1; i--)
-            &info);
+		{
-        for (int i = iu-1; i>= il-1; i--){
+		  printf ("node=%d evals_found=%d evals_tmp[%d] = %g\n", node,
-          printf("node=%d evals_found=%d evals_tmp[%d] = %g\n",node,evals_found, i - (il-1),evals_tmp[i - (il-1)]);
+			  evals_found, i - (il - 1), evals_tmp[i - (il - 1)]);
 		  evals_tmp[i] = evals_tmp[i - (il - 1)];
-          if (il>1) evals_tmp[i-(il-1)]=0.;
+		  if (il > 1)
+		    evals_tmp[i - (il - 1)] = 0.;
 		}
 	    }
 	  {
@@ -442,9 +500,7 @@ public:
 
     void diagonalize (DenseVector < RealD > &lmd,
 		      DenseVector < RealD > &lme,
-		     int N2,
+		      int N2, int N1, GridBase * grid)
-		     int N1,
-		     GridBase *grid)
     {
 
 #ifdef USE_LAPACK
@@ -466,25 +522,27 @@ public:
       return nn;
     }
 
-    void orthogonalize(Field& w,
+    void orthogonalize (Field & w, DenseVector < Field > &evec, int k)
-		       DenseVector<Field>& evec,
-		       int k)
     {
       double t0 = -usecond () / 1e6;
       typedef typename Field::scalar_type MyComplex;
       MyComplex ip;
 
-      if ( 0 ) {
+      if (0)
-	for(int j=0; j<k; ++j){
+	{
+	  for (int j = 0; j < k; ++j)
+	    {
 	      normalise (evec[j]);
-	  for(int i=0;i<j;i++){
+	      for (int i = 0; i < j; i++)
+		{
 		  ip = innerProduct (evec[i], evec[j]);	// are the evecs normalised? ; this assumes so.
 		  evec[j] = evec[j] - ip * evec[i];
 		}
 	    }
 	}
 
-      for(int j=0; j<k; ++j){
+      for (int j = 0; j < k; ++j)
+	{
 	  ip = innerProduct (evec[j], w);	// are the evecs normalised? ; this assumes so.
 	  w = w - ip * evec[j];
 	}
@@ -493,24 +551,27 @@ public:
       OrthoTime += t0;
     }
 
-    void setUnit_Qt(int Nm, DenseVector<RealD> &Qt) {
+    void setUnit_Qt (int Nm, DenseVector < RealD > &Qt)
-      for(int i=0; i<Qt.size(); ++i) Qt[i] = 0.0;
+    {
-      for(int k=0; k<Nm; ++k) Qt[k + k*Nm] = 1.0;
+      for (int i = 0; i < Qt.size (); ++i)
+	Qt[i] = 0.0;
+      for (int k = 0; k < Nm; ++k)
+	Qt[k + k * Nm] = 1.0;
     }
 
 
-    void calc(
+    void calc (DenseVector < RealD > &eval, const Field & src, int &Nconv)
-		DenseVector<RealD>& eval,
-	      const Field& src,
-	      int& Nconv)
     {
 
       GridBase *grid = src._grid;
 //      assert(grid == src._grid);
 
-	std::cout<<GridLogMessage << " -- Nk = " << Nk << " Np = "<< Np << std::endl;
+      std::
+	cout << GridLogMessage << " -- Nk = " << Nk << " Np = " << Np << std::
+	endl;
       std::cout << GridLogMessage << " -- Nm = " << Nm << std::endl;
-	std::cout<<GridLogMessage << " -- size of eval   = " << eval.size() << std::endl;
+      std::cout << GridLogMessage << " -- size of eval   = " << eval.
+	size () << std::endl;
 
 //      assert(c.size() && Nm == eval.size());
 
@@ -530,9 +591,12 @@ public:
       // (uniform vector) Why not src??
       //      evec[0] = 1.0;
       current = src;
-	std:: cout<<GridLogMessage <<"norm2(src)= " << norm2(src)<<std::endl;
+      std::cout << GridLogMessage << "norm2(src)= " << norm2 (src) << std::
+	endl;
       normalise (current);
-	std:: cout<<GridLogMessage <<"norm2(evec[0])= " << norm2(current) <<std::endl;
+      std::
+	cout << GridLogMessage << "norm2(evec[0])= " << norm2 (current) <<
+	std::endl;
 
       // Initial Nk steps
       OrthoTime = 0.;
@@ -541,7 +605,12 @@ public:
 
       uint64_t iter = 0;
 
+      bool initted = false;
+      std::vector < RealD > low (Nstop * 10);
+      std::vector < RealD > high (Nstop * 10);
+      RealD cont = 0.;
       while (1) {
+	  cont = 0.;
 	  std::vector < RealD > lme2 (Nm);
 	  std::vector < RealD > lmd2 (Nm);
 	  for (uint64_t k = 0; k < Nm; ++k, iter++) {
@@ -550,8 +619,12 @@ while(1){
 	      current = next;
 	    }
 	  double t1 = usecond () / 1e6;
-	std::cout<<GridLogMessage <<"IRL::Initial steps: "<<t1-t0<< "seconds"<<std::endl; t0=t1;
+	  std::cout << GridLogMessage << "IRL::Initial steps: " << t1 -
-	std::cout<<GridLogMessage <<"IRL::Initial steps:OrthoTime "<<OrthoTime<< "seconds"<<std::endl;
+	    t0 << "seconds" << std::endl;
+	  t0 = t1;
+	  std::
+	    cout << GridLogMessage << "IRL::Initial steps:OrthoTime " <<
+	    OrthoTime << "seconds" << std::endl;
 
 	  // getting eigenvalues
 	  lmd2.resize (iter + 2);
@@ -561,25 +634,64 @@ while(1){
 	      lme2[k + 2] = lme[k];
 	    }
 	  t1 = usecond () / 1e6;
-	std::cout<<GridLogMessage <<"IRL:: copy: "<<t1-t0<< "seconds"<<std::endl; t0=t1;
+	  std::cout << GridLogMessage << "IRL:: copy: " << t1 -
+	    t0 << "seconds" << std::endl;
+	  t0 = t1;
 	  {
 	    int total = grid->_Nprocessors;
 	    int node = grid->_processor;
 	    int interval = (Nstop / total) + 1;
 	    int iu = (iter + 1) - (interval * node + 1);
 	    int il = (iter + 1) - (interval * (node + 1));
+	    std::vector < RealD > eval2 (iter + 3);
 	    RealD eps2;
-	  Bisection::bisec(lmd2,lme2,iter,il,iu,1e-16,1e-10, eval,eps2);
+	    Bisection::bisec (lmd2, lme2, iter, il, iu, 1e-16, 1e-10, eval2,
+			      eps2);
 //        diagonalize(eval2,lme2,iter,Nk,grid);
-	for(int i=il;i<=iu;i++)
+	    RealD diff = 0.;
-	printf("eval[%d]=%0.14e\n",i,eval[i]);
+	    for (int i = il; i <= iu; i++) {
+		if (initted)
+		  diff =
+		    fabs (eval2[i] - high[iu-i]) / (fabs (eval2[i]) +
+						      fabs (high[iu-i]));
+		if (initted && (diff > eresid))
+		  cont = 1.;
+		if (initted)
+		  printf ("eval[%d]=%0.14e %0.14e, %0.14e\n", i, eval2[i],
+			  high[iu-i], diff);
+		high[iu-i] = eval2[i];
+	      }
+	    il = (interval * node + 1);
+	    iu = (interval * (node + 1));
+	    Bisection::bisec (lmd2, lme2, iter, il, iu, 1e-16, 1e-10, eval2,
+			      eps2);
+	    for (int i = il; i <= iu; i++) {
+		if (initted)
+		  diff =
+		    fabs (eval2[i] - low[i]) / (fabs (eval2[i]) +
+						fabs (low[i]));
+		if (initted && (diff > eresid))
+		  cont = 1.;
+		if (initted)
+		  printf ("eval[%d]=%0.14e %0.14e, %0.14e\n", i, eval2[i],
+			  low[i], diff);
+		low[i] = eval2[i];
+	      }
 	    t1 = usecond () / 1e6;
-	std::cout<<GridLogMessage <<"IRL:: diagonalize: "<<t1-t0<< "seconds"<<std::endl; t0=t1;
+	    std::cout << GridLogMessage << "IRL:: diagonalize: " << t1 -
+	      t0 << "seconds" << std::endl;
+	    t0 = t1;
 	  }
 
 	  for (uint64_t k = 0; k < Nk; ++k) {
 //          eval[k] = eval2[k];
 	    }
+	  if (initted)
+	    {
+	      grid->GlobalSumVector (&cont, 1);
+	      if (cont < 1.) return;
+	    }
+	  initted = true;
 	}
 
     }
@@ -594,28 +706,35 @@ while(1){
    Q.e_1 = y and Q is unitary.
 **/
     template < class T >
-  static T orthQ(DenseMatrix<T> &Q, DenseVector<T> y){
+      static T orthQ (DenseMatrix < T > &Q, DenseVector < T > y)
+    {
       int N = y.size ();	//Matrix Size
       Fill (Q, 0.0);
       T tau;
-    for(int i=0;i<N;i++){
+      for (int i = 0; i < N; i++)
+	{
 	  Q[i][0] = y[i];
 	}
       T sig = conj (y[0]) * y[0];
       T tau0 = fabs (sqrt (sig));
 
-    for(int j=1;j<N;j++){
+      for (int j = 1; j < N; j++)
+	{
 	  sig += conj (y[j]) * y[j];
 	  tau = abs (sqrt (sig));
 
-      if(abs(tau0) > 0.0){
+	  if (abs (tau0) > 0.0)
+	    {
 
 	      T gam = conj ((y[j] / tau) / tau0);
-	for(int k=0;k<=j-1;k++){  
+	      for (int k = 0; k <= j - 1; k++)
+		{
 		  Q[k][j] = -gam * y[k];
 		}
 	      Q[j][j] = tau0 / tau;
-      } else {
+	    }
+	  else
+	    {
 	      Q[j - 1][j] = 1.0;
 	    }
 	  tau0 = tau;
@@ -628,7 +747,8 @@ while(1){
 	Q.e_k = y and Q is unitary.
 **/
     template < class T >
-  static T orthU(DenseMatrix<T> &Q, DenseVector<T> y){
+      static T orthU (DenseMatrix < T > &Q, DenseVector < T > y)
+    {
       T tau = orthQ (Q, y);
       SL (Q);
       return tau;
@@ -645,37 +765,45 @@ say con = 2
 
 **/
 
-template<class T>
+    template < class T > static void Lock (DenseMatrix < T > &H,	///Hess mtx     
-static void Lock(DenseMatrix<T> &H, 	///Hess mtx	
 					   DenseMatrix < T > &Q,	///Lock Transform
 					   T val,	///value to be locked
 					   int con,	///number already locked
-		 RealD small,
+					   RealD small, int dfg, bool herm)
-		 int dfg,
-		 bool herm)
     {
       //ForceTridiagonal(H);
 
       int M = H.dim;
-  DenseVector<T> vec; Resize(vec,M-con);
+      DenseVector < T > vec;
+      Resize (vec, M - con);
 
-  DenseMatrix<T> AH; Resize(AH,M-con,M-con);
+      DenseMatrix < T > AH;
+      Resize (AH, M - con, M - con);
       AH = GetSubMtx (H, con, M, con, M);
 
-  DenseMatrix<T> QQ; Resize(QQ,M-con,M-con);
+      DenseMatrix < T > QQ;
+      Resize (QQ, M - con, M - con);
 
-  Unity(Q);   Unity(QQ);
+      Unity (Q);
+      Unity (QQ);
 
-  DenseVector<T> evals; Resize(evals,M-con);
+      DenseVector < T > evals;
-  DenseMatrix<T> evecs; Resize(evecs,M-con,M-con);
+      Resize (evals, M - con);
+      DenseMatrix < T > evecs;
+      Resize (evecs, M - con, M - con);
 
       Wilkinson < T > (AH, evals, evecs, small);
 
       int k = 0;
       RealD cold = abs (val - evals[k]);
-  for(int i=1;i<M-con;i++){
+      for (int i = 1; i < M - con; i++)
+	{
 	  RealD cnew = abs (val - evals[i]);
-    if( cnew < cold ){k = i; cold = cnew;}
+	  if (cnew < cold)
+	    {
+	      k = i;
+	      cold = cnew;
+	    }
 	}
       vec = evecs[k];
 
@@ -686,80 +814,110 @@ static void Lock(DenseMatrix<T> &H, 	///Hess mtx
       AH = Hermitian (QQ) * AH;
       AH = AH * QQ;
 
-  for(int i=con;i<M;i++){
+      for (int i = con; i < M; i++)
-    for(int j=con;j<M;j++){
+	{
+	  for (int j = con; j < M; j++)
+	    {
 	      Q[i][j] = QQ[i - con][j - con];
 	      H[i][j] = AH[i - con][j - con];
 	    }
 	}
 
-  for(int j = M-1; j>con+2; j--){
+      for (int j = M - 1; j > con + 2; j--)
+	{
 
-    DenseMatrix<T> U; Resize(U,j-1-con,j-1-con);
+	  DenseMatrix < T > U;
-    DenseVector<T> z; Resize(z,j-1-con); 
+	  Resize (U, j - 1 - con, j - 1 - con);
+	  DenseVector < T > z;
+	  Resize (z, j - 1 - con);
 	  T nm = norm (z);
-    for(int k = con+0;k<j-1;k++){
+	  for (int k = con + 0; k < j - 1; k++)
+	    {
 	      z[k - con] = conj (H (j, k + 1));
 	    }
 	  normalise (z);
 
 	  RealD tmp = 0;
-    for(int i=0;i<z.size()-1;i++){tmp = tmp + abs(z[i]);}
+	  for (int i = 0; i < z.size () - 1; i++)
+	    {
+	      tmp = tmp + abs (z[i]);
+	    }
 
-    if(tmp < small/( (RealD)z.size()-1.0) ){ continue;}	
+	  if (tmp < small / ((RealD) z.size () - 1.0))
+	    {
+	      continue;
+	    }
 
 	  tau = orthU (U, z);
 
-    DenseMatrix<T> Hb; Resize(Hb,j-1-con,M);	
+	  DenseMatrix < T > Hb;
+	  Resize (Hb, j - 1 - con, M);
 
-    for(int a = 0;a<M;a++){
+	  for (int a = 0; a < M; a++)
-      for(int b = 0;b<j-1-con;b++){
+	    {
+	      for (int b = 0; b < j - 1 - con; b++)
+		{
 		  T sum = 0;
-	for(int c = 0;c<j-1-con;c++){
+		  for (int c = 0; c < j - 1 - con; c++)
+		    {
 		      sum += H[a][con + 1 + c] * U[c][b];
 		    }		//sum += H(a,con+1+c)*U(c,b);}
 		  Hb[b][a] = sum;
 		}
 	    }
 
-    for(int k=con+1;k<j;k++){
+	  for (int k = con + 1; k < j; k++)
-      for(int l=0;l<M;l++){
+	    {
+	      for (int l = 0; l < M; l++)
+		{
 		  H[l][k] = Hb[k - 1 - con][l];
 		}
 	    }			//H(Hb[k-1-con][l] , l,k);}}
 
-    DenseMatrix<T> Qb; Resize(Qb,M,M);	
+	  DenseMatrix < T > Qb;
+	  Resize (Qb, M, M);
 
-    for(int a = 0;a<M;a++){
+	  for (int a = 0; a < M; a++)
-      for(int b = 0;b<j-1-con;b++){
+	    {
+	      for (int b = 0; b < j - 1 - con; b++)
+		{
 		  T sum = 0;
-	for(int c = 0;c<j-1-con;c++){
+		  for (int c = 0; c < j - 1 - con; c++)
+		    {
 		      sum += Q[a][con + 1 + c] * U[c][b];
 		    }		//sum += Q(a,con+1+c)*U(c,b);}
 		  Qb[b][a] = sum;
 		}
 	    }
 
-    for(int k=con+1;k<j;k++){
+	  for (int k = con + 1; k < j; k++)
-      for(int l=0;l<M;l++){
+	    {
+	      for (int l = 0; l < M; l++)
+		{
 		  Q[l][k] = Qb[k - 1 - con][l];
 		}
 	    }			//Q(Qb[k-1-con][l] , l,k);}}
 
-    DenseMatrix<T> Hc; Resize(Hc,M,M);	
+	  DenseMatrix < T > Hc;
+	  Resize (Hc, M, M);
 
-    for(int a = 0;a<j-1-con;a++){
+	  for (int a = 0; a < j - 1 - con; a++)
-      for(int b = 0;b<M;b++){
+	    {
+	      for (int b = 0; b < M; b++)
+		{
 		  T sum = 0;
-	for(int c = 0;c<j-1-con;c++){
+		  for (int c = 0; c < j - 1 - con; c++)
+		    {
 		      sum += conj (U[c][a]) * H[con + 1 + c][b];
 		    }		//sum += conj( U(c,a) )*H(con+1+c,b);}
 		  Hc[b][a] = sum;
 		}
 	    }
 
-    for(int k=0;k<M;k++){
+	  for (int k = 0; k < M; k++)
-      for(int l=con+1;l<j;l++){
+	    {
+	      for (int l = con + 1; l < j; l++)
+		{
 		  H[l][k] = Hc[k][l - 1 - con];
 		}
 	    }			//H(Hc[k][l-1-con] , l,k);}}
@@ -773,4 +931,3 @@ static void Lock(DenseMatrix<T> &H, 	///Hess mtx
 
 }
 #endif
-