diff options
author | Thomas White <taw@physics.org> | 2015-03-13 11:06:28 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2015-03-13 11:06:28 +0100 |
commit | 9225115a425ced041133e849f3816178967e3307 (patch) | |
tree | fd7972cc6ca70394cd530a1ce110dbccfb496c05 /src | |
parent | 9d24d9e4329389efbbdd33da95010005ba9ea585 (diff) |
Move solve_svd() to utils
Diffstat (limited to 'src')
-rw-r--r-- | src/post-refinement.c | 134 |
1 files changed, 0 insertions, 134 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c index 4dae5096..96c0de12 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -362,140 +362,6 @@ static void apply_shift(Crystal *cr, int k, double shift) } -static int check_eigen(gsl_vector *e_val, int verbose) -{ - int i; - double vmax, vmin; - const int n = e_val->size; - const double max_condition = 1e6; - int n_filt = 0; - - if ( verbose ) STATUS("Eigenvalues:\n"); - vmin = +INFINITY; - vmax = 0.0; - for ( i=0; i<n; i++ ) { - double val = gsl_vector_get(e_val, i); - if ( verbose ) STATUS("%i: %e\n", i, val); - if ( val > vmax ) vmax = val; - if ( val < vmin ) vmin = val; - } - - for ( i=0; i<n; i++ ) { - double val = gsl_vector_get(e_val, i); - if ( val < vmax/max_condition ) { - gsl_vector_set(e_val, i, 0.0); - n_filt++; - } - } - - vmin = +INFINITY; - vmax = 0.0; - for ( i=0; i<n; i++ ) { - double val = gsl_vector_get(e_val, i); - if ( val == 0.0 ) continue; - if ( val > vmax ) vmax = val; - if ( val < vmin ) vmin = val; - } - if ( verbose ) { - STATUS("Condition number: %e / %e = %5.2f\n", - vmax, vmin, vmax/vmin); - STATUS("%i out of %i eigenvalues filtered.\n", n_filt, n); - } - - return n_filt; -} - - -static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt, - int verbose) -{ - gsl_matrix *s_vec; - gsl_vector *s_val; - int err, n; - gsl_vector *shifts; - gsl_vector *SB; - gsl_vector *SinvX; - gsl_matrix *S; /* rescaling matrix due to Bricogne */ - gsl_matrix *AS; - gsl_matrix *SAS; - int i; - gsl_matrix *SAS_copy; - - n = v->size; - if ( v->size != M->size1 ) return NULL; - if ( v->size != M->size2 ) return NULL; - - /* Calculate the rescaling matrix S */ - S = gsl_matrix_calloc(n, n); - for ( i=0; i<n; i++ ) { - double sii = pow(gsl_matrix_get(M, i, i), -0.5); - gsl_matrix_set(S, i, i, sii); - } - - /* Calculate the matrix SAS, which we will be (not) inverting */ - AS = gsl_matrix_calloc(n, n); - SAS = gsl_matrix_calloc(n, n); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, M, S, 0.0, AS); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, S, AS, 0.0, SAS); - gsl_matrix_free(AS); - - /* Calculate the vector SB, which is the RHS of the equation */ - SB = gsl_vector_calloc(n); - gsl_blas_dgemv(CblasNoTrans, 1.0, S, v, 0.0, SB); - - if ( verbose ) { - STATUS("The equation after rescaling:\n"); - show_matrix_eqn(SAS, SB); - } - - SAS_copy = gsl_matrix_alloc(SAS->size1, SAS->size2); - gsl_matrix_memcpy(SAS_copy, SAS); - - /* Do the SVD */ - s_val = gsl_vector_calloc(n); - s_vec = gsl_matrix_calloc(n, n); - err = gsl_linalg_SV_decomp_jacobi(SAS, s_vec, s_val); - if ( err ) { - if ( verbose ) ERROR("SVD failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - gsl_matrix_free(SAS); - gsl_matrix_free(S); - return NULL; - } - /* "SAS" is now "U" */ - - /* Filter the eigenvalues */ - *n_filt = check_eigen(s_val, verbose); - - gsl_matrix_free(SAS_copy); - - /* Solve the equation SAS.SinvX = SB */ - SinvX = gsl_vector_calloc(n); - err = gsl_linalg_SV_solve(SAS, s_vec, s_val, SB, SinvX); - gsl_vector_free(SB); - gsl_matrix_free(SAS); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - - if ( err ) { - ERROR("Matrix solution failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(S); - gsl_vector_free(SinvX); - return NULL; - } - - /* Calculate S.SinvX to get X, the shifts */ - shifts = gsl_vector_calloc(n); - gsl_blas_dgemv(CblasNoTrans, 1.0, S, SinvX, 0.0, shifts); - - gsl_matrix_free(S); - gsl_vector_free(SinvX); - - return shifts; -} - - /* Perform one cycle of post refinement on 'image' against 'full' */ static double pr_iterate(Crystal *cr, const RefList *full, PartialityModel pmodel, int *n_filtered) |