diff options
| -rw-r--r-- | src/post-refinement.c | 41 |
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++ ) { |