aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-06-21 17:46:09 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:29 +0100
commitdd214964cef3b552d1d4a7abf47eb5893f3d7ca0 (patch)
treef29ccdc4e50c7c6b4c21ed6e6a745231537dfd0c
parentd01cb93fa904b6f8f5b32953e547776f1797a14f (diff)
Add SVD matrix solution
-rw-r--r--src/post-refinement.c41
1 files changed, 41 insertions, 0 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 5be7cf41..7301b57b 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -270,6 +270,46 @@ static gsl_vector *solve_by_eigenvalue_filtration(gsl_vector *v, gsl_matrix *M)
}
+static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M)
+{
+ gsl_matrix *s_vec;
+ gsl_vector *s_val;
+ int err, n;
+ gsl_vector *shifts;
+
+ n = v->size;
+ if ( v->size != M->size1 ) return NULL;
+ if ( v->size != M->size2 ) return NULL;
+
+ s_val = gsl_vector_calloc(n);
+ s_vec = gsl_matrix_calloc(n, n);
+
+ err = gsl_linalg_SV_decomp_jacobi(M, s_vec, s_val);
+ if ( err ) {
+ ERROR("SVD failed: %s\n", gsl_strerror(err));
+ gsl_matrix_free(s_vec);
+ gsl_vector_free(s_val);
+ return NULL;
+ }
+ /* "M" is now "U" */
+
+ shifts = gsl_vector_calloc(n);
+ err = gsl_linalg_SV_solve(M, s_vec, s_val, v, shifts);
+ if ( err ) {
+ ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
+ gsl_matrix_free(s_vec);
+ gsl_vector_free(s_val);
+ gsl_vector_free(shifts);
+ return NULL;
+ }
+
+ gsl_matrix_free(s_vec);
+ gsl_vector_free(s_val);
+
+ return shifts;
+}
+
+
static gsl_vector *solve_householder(gsl_vector *v, gsl_matrix *M)
{
int n, err;
@@ -387,6 +427,7 @@ 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);
if ( shifts != NULL ) {
for ( param=0; param<NUM_PARAMS; param++ ) {