aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-06-22 11:59:25 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:29 +0100
commit93fd0cbb09b71179ff0407a9d04c05882d70c385 (patch)
tree8cff8f5b1e0e4276c68a5bb43276f3f3a9e2669e
parenta1e023343b085f727f368720fc7cda8e0e4a2724 (diff)
Use SVD for matrix solution
It's identical to the previous eigenvalue filtration thing
-rw-r--r--src/post-refinement.c70
1 files changed, 18 insertions, 52 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 64704598..c1bc9314 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -216,57 +216,22 @@ static void apply_shift(struct image *image, int k, double shift)
}
-/* NB This is DIFFERENT to the solve_by_eigenvalue_filtration routing in
- * hrs-scaling.c */
-static gsl_vector *solve_by_eigenvalue_filtration(gsl_vector *v, gsl_matrix *M)
+static void show_eigen(gsl_vector *e_val)
{
- gsl_eigen_symmv_workspace *work;
- gsl_vector *e_val;
- gsl_matrix *e_vec;
- int val, n, frame;
- gsl_vector *shifts;
-
- n = v->size;
- if ( v->size != M->size1 ) return NULL;
- if ( v->size != M->size2 ) return NULL;
-
- /* Diagonalise */
- work = gsl_eigen_symmv_alloc(n);
- e_val = gsl_vector_alloc(n);
- e_vec = gsl_matrix_alloc(n, n);
- val = gsl_eigen_symmv(M, e_val, e_vec, work);
- if ( val ) {
- ERROR("Couldn't diagonalise matrix.\n");
- return NULL;
- }
- gsl_eigen_symmv_free(work);
- val = gsl_eigen_symmv_sort(e_val, e_vec, GSL_EIGEN_SORT_ABS_DESC);
-
- /* Rotate the "solution vector" */
- gsl_vector *rprime;
- rprime = gsl_vector_alloc(n);
- val = gsl_blas_dgemv(CblasTrans, 1.0, e_vec, v, 0.0, rprime);
-
- /* Solve (now easy) */
- gsl_vector *sprime;
- sprime = gsl_vector_alloc(n);
- for ( frame=0; frame<n; frame++ ) {
- double num, den;
- num = gsl_vector_get(rprime, frame);
- den = gsl_vector_get(e_val, frame);
- gsl_vector_set(sprime, frame, num/den);
+ int i;
+ double vmax, vmin;
+ const int n = e_val->size;
+
+ STATUS("Eigenvalues:\n");
+ vmin = +INFINITY;
+ vmax = 0.0;
+ for ( i=0; i<n; i++ ) {
+ double val = gsl_vector_get(e_val, i);
+ STATUS("%i: %e\n", i, val);
+ if ( val > vmax ) vmax = val;
+ if ( val < vmin ) vmin = val;
}
-
- /* Rotate back */
- shifts = gsl_vector_alloc(n);
- val = gsl_blas_dgemv(CblasNoTrans, 1.0, e_vec, sprime, 0.0, shifts);
-
- gsl_vector_free(sprime);
- gsl_vector_free(rprime);
- gsl_matrix_free(e_vec);
- gsl_vector_free(e_val);
-
- return shifts;
+ STATUS("Condition number: %e / %e = %e\n", vmin, vmax, vmin/vmax);
}
@@ -293,6 +258,8 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M)
}
/* "M" is now "U" */
+ show_eigen(s_val);
+
shifts = gsl_vector_calloc(n);
err = gsl_linalg_SV_solve(M, s_vec, s_val, v, shifts);
if ( err ) {
@@ -425,9 +392,8 @@ static double pr_iterate(struct image *image, const RefList *full,
}
max_shift = 0.0;
- //shifts = solve_by_eigenvalue_filtration(v, M);
- shifts = solve_householder(v, M);
- //shifts = solve_svd(v, M);
+ //shifts = solve_householder(v, M);
+ shifts = solve_svd(v, M);
if ( shifts != NULL ) {
for ( param=0; param<NUM_PARAMS; param++ ) {