Simplified lanczos, added Eigen diagonalisation.

Curious if we can deprecate dependencly on BLAS.
Will see when we get 48^3 running on our BG/Q port

lib/algorithms/iterative/BlockConjugateGradient.h

@@ -56,11 +56,8 @@ class BlockConjugateGradient : public OperatorFunction<Field> {
   Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
   
   BlockConjugateGradient(BlockCGtype cgtype,int _Orthog,RealD tol, Integer maxit, bool err_on_no_conv = true)
-    : Tolerance(tol),
+    : Tolerance(tol), CGtype(cgtype),   blockDim(_Orthog),  MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv)
-    CGtype(cgtype),
+  {};
-    blockDim(_Orthog),
-    MaxIterations(maxit),
-    ErrorOnNoConverge(err_on_no_conv){};
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Thin QR factorisation (google it)

lib/algorithms/iterative/EigenSort.h

@@ -1,81 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/algorithms/iterative/EigenSort.h
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
-#ifndef GRID_EIGENSORT_H
-#define GRID_EIGENSORT_H
-
-
-namespace Grid {
-    /////////////////////////////////////////////////////////////
-    // Eigen sorter to begin with
-    /////////////////////////////////////////////////////////////
-
-template<class Field>
-class SortEigen {
- private:
-  
-//hacking for testing for now
- private:
-  static bool less_lmd(RealD left,RealD right){
-    return left > right;
-  }  
-  static bool less_pair(std::pair<RealD,Field const*>& left,
-                        std::pair<RealD,Field const*>& right){
-    return left.first > (right.first);
-  }  
-  
-  
- public:
-
-  void push(DenseVector<RealD>& lmd,
-            DenseVector<Field>& evec,int N) {
-    DenseVector<Field> cpy(lmd.size(),evec[0]._grid);
-    for(int i=0;i<lmd.size();i++) cpy[i] = evec[i];
-    
-    DenseVector<std::pair<RealD, Field const*> > emod(lmd.size());    
-    for(int i=0;i<lmd.size();++i)
-      emod[i] = std::pair<RealD,Field const*>(lmd[i],&cpy[i]);
-
-    partial_sort(emod.begin(),emod.begin()+N,emod.end(),less_pair);
-
-    typename DenseVector<std::pair<RealD, Field const*> >::iterator it = emod.begin();
-    for(int i=0;i<N;++i){
-      lmd[i]=it->first;
-      evec[i]=*(it->second);
-      ++it;
-    }
-  }
-  void push(DenseVector<RealD>& lmd,int N) {
-    std::partial_sort(lmd.begin(),lmd.begin()+N,lmd.end(),less_lmd);
-  }
-  bool saturated(RealD lmd, RealD thrs) {
-    return fabs(lmd) > fabs(thrs);
-  }
-};
-
-}
-#endif

tests/solver/Test_dwf_lanczos.cc

@@ -92,16 +92,15 @@ int main (int argc, char ** argv)
 
   
   std::vector<RealD>          eval(Nm);
-  FermionField    src(FrbGrid); gaussian(RNG5rb,src);
+  FermionField    src(FrbGrid); 
+  gaussian(RNG5rb,src);
   std::vector<FermionField> evec(Nm,FrbGrid);
   for(int i=0;i<1;i++){
-    std::cout << i<<" / "<< Nm<< " grid pointer "<<evec[i]._grid<<std::endl;
+    std::cout << GridLogMessage <<i<<" / "<< Nm<< " grid pointer "<<evec[i]._grid<<std::endl;
   };
 
   int Nconv;
-  IRL.calc(eval,evec,
+  IRL.calc(eval,evec,src,Nconv);
-	   src,
-	   Nconv);
 
 
   Grid_finalize();