aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-06-22 16:01:17 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:30 +0100
commitd72134811702f08888035a743b9d0609610c4de7 (patch)
tree1a0b2be3e809258f09b0c1632c794f751bb608ed
parentdc9d2400e3c1ab3a3362cd2aa9243260f3dfdbf4 (diff)
Avoid shifts along eigenvectors with small eigenvalues
-rw-r--r--src/post-refinement.c19
1 files changed, 18 insertions, 1 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index a6dd5f75..a48a2052 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -221,6 +221,7 @@ static void show_eigen(gsl_vector *e_val)
int i;
double vmax, vmin;
const int n = e_val->size;
+ const double max_condition = 1e6;
STATUS("Eigenvalues:\n");
vmin = +INFINITY;
@@ -231,7 +232,23 @@ static void show_eigen(gsl_vector *e_val)
if ( val > vmax ) vmax = val;
if ( val < vmin ) vmin = val;
}
- STATUS("Condition number: %e / %e = %e\n", vmin, vmax, vmin/vmax);
+
+ for ( i=0; i<n; i++ ) {
+ double val = gsl_vector_get(e_val, i);
+ if ( val < vmax/max_condition ) {
+ gsl_vector_set(e_val, i, 0.0);
+ }
+ }
+
+ vmin = +INFINITY;
+ vmax = 0.0;
+ for ( i=0; i<n; i++ ) {
+ double val = gsl_vector_get(e_val, i);
+ if ( val == 0.0 ) continue;
+ if ( val > vmax ) vmax = val;
+ if ( val < vmin ) vmin = val;
+ }
+ STATUS("Condition number: %e / %e = %5.2f\n", vmax, vmin, vmax/vmin);
}