diff options
Diffstat (limited to 'src/post-refinement.c')
-rw-r--r-- | src/post-refinement.c | 51 |
1 files changed, 1 insertions, 50 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c index 2c19811d..0d6aad25 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -358,55 +358,6 @@ static int check_eigen(gsl_vector *e_val, int verbose) } -static gsl_vector *solve_svd_noscale(gsl_vector *v, gsl_matrix *Mp, int *n_filt, - int verbose) -{ - gsl_matrix *s_vec; - gsl_vector *s_val; - int err, n; - gsl_vector *shifts; - gsl_matrix *M; - - n = v->size; - if ( v->size != Mp->size1 ) return NULL; - if ( v->size != Mp->size2 ) return NULL; - - M = gsl_matrix_alloc(n, n); - if ( M == NULL ) return NULL; - gsl_matrix_memcpy(M, Mp); - - 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" */ - - *n_filt = check_eigen(s_val, verbose); - - 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); - gsl_matrix_free(M); - - return shifts; -} - - static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt, int verbose) { @@ -605,7 +556,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, } max_shift = 0.0; - shifts = solve_svd_noscale(v, M, &prdata->n_filtered, verbose); + shifts = solve_svd(v, M, &prdata->n_filtered, verbose); if ( shifts != NULL ) { for ( param=0; param<NUM_PARAMS; param++ ) { |