From bff20f75d4691fe648a9c2e999aeb0d483d28d7c Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Thu, 30 Jul 2015 14:57:08 +0200
Subject: [PATCH] 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 | 9 ++++++++-
1 file changed, 8 insertions(+), 1 deletion(-)
diff --git a/src/tools/LatticeReduction.cpp b/src/tools/LatticeReduction.cpp
index 5cbbb2364..db010fe94 100644
--- a/src/tools/LatticeReduction.cpp
+++ b/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);
--
GitLab