Changed tolerance in lattice reduction

Lattice reduction is used in plumed to convert box to its
"least skewed" form. I changed a bit the threashold since it
was giving problems with high optimization flags.
Additionally, I added some warning when too many iterations
are performed (which indicates a numerical issue)

src/tools/LatticeReduction.cpp

@@ -21,10 +21,11 @@
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
 #include "Exception.h"
 #include "LatticeReduction.h"
+#include <stdio.h>
 
 namespace PLMD{
 
-const double epsilon=1e-28;
+const double epsilon=1e-14;
 
 void LatticeReduction::sort(Vector v[3]){
   for(int i=0;i<3;i++) for(int j=i+1;j<3;j++) if(modulo2(v[i])>modulo2(v[j])){
@@ -36,6 +37,7 @@ void LatticeReduction::sort(Vector v[3]){
 void LatticeReduction::reduce(Vector&a,Vector&b){
   double ma=modulo2(a);
   double mb=modulo2(b);
+  unsigned counter=0;
   while(true){
     if(mb>ma){
       Vector t(a); a=b; b=t;
@@ -44,6 +46,8 @@ void LatticeReduction::reduce(Vector&a,Vector&b){
     a-=b*floor(dotProduct(a,b)/modulo2(b)+0.5);
     ma=modulo2(a);
     if(mb<=ma+epsilon) break;
+    counter++;
+    if(counter%10000==0) fprintf(stderr,"WARNING: LatticeReduction::reduce stuck after %u iterations\n",counter);
   }
 
   Vector t(a); a=b; b=t;
@@ -100,6 +104,7 @@ void LatticeReduction::reduceFast(Tensor&t){
   v[0]=t.getRow(0);
   v[1]=t.getRow(1);
   v[2]=t.getRow(2);
+  unsigned counter=0;
   while(true){
     sort(v);
     reduce(v[0],v[1]);
@@ -129,6 +134,8 @@ void LatticeReduction::reduceFast(Tensor&t){
       }
     }
     if(modulo2(best)+epsilon>=modulo2(v[2])) break;
+    counter++;
+    if(counter%10000==0) fprintf(stderr,"WARNING: LatticeReduction::reduceFast stuck after %u iterations\n",counter);
     v[2]=best;
   }
   sort(v);