From 93fd0cbb09b71179ff0407a9d04c05882d70c385 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Wed, 22 Jun 2011 11:59:25 +0200
Subject: Use SVD for matrix solution

It's identical to the previous eigenvalue filtration thing
---
 src/post-refinement.c | 70 +++++++++++++--------------------------------------
 1 file changed, 18 insertions(+), 52 deletions(-)

(limited to 'src/post-refinement.c')

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++ ) {
-- 
cgit v1.2.3