From 0c1f39250f05185f82e943c1622bbda411f8962f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 11 Oct 2013 16:34:23 +0200 Subject: Add solve_svd_noscale() --- src/post-refinement.c | 55 ++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 52 insertions(+), 3 deletions(-) diff --git a/src/post-refinement.c b/src/post-refinement.c index 43cc002b..2c19811d 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -314,7 +314,7 @@ static void apply_shift(Crystal *cr, int k, double shift) } -static int check_eigen(gsl_matrix *M, gsl_vector *e_val, int verbose) +static int check_eigen(gsl_vector *e_val, int verbose) { int i; double vmax, vmin; @@ -358,6 +358,55 @@ static int check_eigen(gsl_matrix *M, 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) { @@ -418,7 +467,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt, /* "SAS" is now "U" */ /* Filter the eigenvalues */ - *n_filt = check_eigen(SAS_copy, s_val, verbose); + *n_filt = check_eigen(s_val, verbose); gsl_matrix_free(SAS_copy); @@ -556,7 +605,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, } max_shift = 0.0; - shifts = solve_svd(v, M, &prdata->n_filtered, verbose); + shifts = solve_svd_noscale(v, M, &prdata->n_filtered, verbose); if ( shifts != NULL ) { for ( param=0; param