diff options
| author | Thomas White <taw@physics.org> | 2011-06-22 16:01:17 +0200 |
|---|---|---|
| committer | Thomas White <taw@physics.org> | 2012-02-22 15:27:30 +0100 |
| commit | d72134811702f08888035a743b9d0609610c4de7 (patch) | |
| tree | 1a0b2be3e809258f09b0c1632c794f751bb608ed | |
| parent | dc9d2400e3c1ab3a3362cd2aa9243260f3dfdbf4 (diff) | |
Avoid shifts along eigenvectors with small eigenvalues
| -rw-r--r-- | src/post-refinement.c | 19 |
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); } |