aboutsummaryrefslogtreecommitdiff
path: root/src/post-refinement.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-06-21 15:26:31 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:29 +0100
commitd01cb93fa904b6f8f5b32953e547776f1797a14f (patch)
treea9c6b63285956fe83f33e348732aa5aa3f3519f5 /src/post-refinement.c
parent019234358834610a0d7bb97a3587b8b0e1c779aa (diff)
"Pluggable" matrix solvers
Diffstat (limited to 'src/post-refinement.c')
-rw-r--r--src/post-refinement.c99
1 files changed, 96 insertions, 3 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index f8b66593..5be7cf41 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -19,6 +19,8 @@
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_blas.h>
#include "image.h"
#include "post-refinement.h"
@@ -214,6 +216,83 @@ 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)
+{
+ 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);
+ }
+
+ /* 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;
+}
+
+
+static gsl_vector *solve_householder(gsl_vector *v, gsl_matrix *M)
+{
+ int n, err;
+ gsl_vector *shifts;
+
+ n = v->size;
+ if ( v->size != M->size1 ) return NULL;
+ if ( v->size != M->size2 ) return NULL;
+
+ shifts = gsl_vector_alloc(n);
+
+ err = gsl_linalg_HH_solve(M, v, shifts);
+
+ if ( err != 0 ) {
+ ERROR("Failed to solve equations: %s\n", gsl_strerror(err));
+ gsl_vector_free(shifts);
+ return NULL;
+ }
+
+ return shifts;
+}
+
+
/* Perform one cycle of post refinement on 'image' against 'full' */
static double pr_iterate(struct image *image, const RefList *full,
const char *sym)
@@ -226,6 +305,7 @@ static double pr_iterate(struct image *image, const RefList *full,
RefListIterator *iter;
RefList *reflections;
double max_shift;
+ int nref = 0;
reflections = image->reflections;
@@ -293,20 +373,33 @@ static double pr_iterate(struct image *image, const RefList *full,
}
+ nref++;
+
}
//show_matrix_eqn(M, v, NUM_PARAMS);
- shifts = gsl_vector_alloc(NUM_PARAMS);
+ STATUS("%i reflections were scalable\n", nref);
+ if ( nref == 0 ) {
+ ERROR("No reflections left to scale!\n");
+ return 0.0;
+ }
+
max_shift = 0.0;
- if ( gsl_linalg_HH_solve(M, v, shifts) == 0 ) {
+ //shifts = solve_by_eigenvalue_filtration(v, M);
+ shifts = solve_householder(v, M);
+ if ( shifts != NULL ) {
for ( param=0; param<NUM_PARAMS; param++ ) {
double shift = gsl_vector_get(shifts, param);
apply_shift(image, param, shift);
+ //STATUS("Shift %i: %e\n", param, shift);
if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
}
- } /* else problem with the equations, so leave things as they were */
+ } else {
+ ERROR("Problem solving equations.\n");
+ /* Leave things as they were */
+ }
gsl_matrix_free(M);
gsl_vector_free(v);