diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/partialator.c | 25 | ||||
| -rw-r--r-- | src/post-refinement.c | 22 | ||||
| -rw-r--r-- | src/post-refinement.h | 2 | ||||
| -rw-r--r-- | src/scaling-report.c | 3 |
4 files changed, 45 insertions, 7 deletions
diff --git a/src/partialator.c b/src/partialator.c index 88912226..ad6d9a49 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -309,6 +309,16 @@ static void display_progress(int n_images, int n_crystals) } +static const char *str_flags(Crystal *cr) +{ + if ( crystal_get_user_flag(cr) ) { + return "N"; + } + + return "-"; +} + + int main(int argc, char *argv[]) { int c; @@ -671,6 +681,21 @@ int main(int argc, char *argv[]) /* Output results */ write_reflist(outfile, full); + /* Dump parameters */ + FILE *fh; + fh = fopen("partialator.params", "w"); + if ( fh == NULL ) { + ERROR("Couldn't open partialator.params!\n"); + } else { + for ( i=0; i<n_crystals; i++ ) { + fprintf(fh, "%4i %5.2f %8.5e %s\n", i, + crystal_get_osf(crystals[i]), + crystal_get_image(crystals[i])->div, + str_flags(crystals[i])); + } + fclose(fh); + } + /* Clean up */ for ( i=0; i<n_crystals; i++ ) { reflist_free(crystal_get_reflections(crystals[i])); diff --git a/src/post-refinement.c b/src/post-refinement.c index 94c8a816..88c609c1 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -314,13 +314,12 @@ static void apply_shift(Crystal *cr, int k, double shift) } -static int check_eigen(gsl_vector *e_val) +static int check_eigen(gsl_matrix *M, gsl_vector *e_val, int verbose) { int i; double vmax, vmin; const int n = e_val->size; const double max_condition = 1e6; - const int verbose = 0; int n_filt = 0; if ( verbose ) STATUS("Eigenvalues:\n"); @@ -374,6 +373,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt, gsl_matrix *AS; gsl_matrix *SAS; int i; + gsl_matrix *SAS_copy; n = v->size; if ( v->size != M->size1 ) return NULL; @@ -393,6 +393,18 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt, 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); @@ -408,11 +420,9 @@ 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(s_val); + *n_filt = check_eigen(SAS_copy, s_val, verbose); - /* 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); + gsl_matrix_free(SAS_copy); /* Solve the equation SAS.SinvX = SB */ SinvX = gsl_vector_calloc(n); diff --git a/src/post-refinement.h b/src/post-refinement.h index b964975e..63b6ff66 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -56,8 +56,8 @@ enum { REF_CSX, REF_CSY, REF_CSZ, - NUM_PARAMS, REF_DIV, + NUM_PARAMS, REF_R, }; diff --git a/src/scaling-report.c b/src/scaling-report.c index acb007ae..17268489 100644 --- a/src/scaling-report.c +++ b/src/scaling-report.c @@ -889,6 +889,9 @@ void sr_iteration(SRContext *sr, int iteration, struct srdata *d) cairo_restore(sr->cr); } + + STATUS("%i filtered total, %5.2f filtered per crystal\n", + d->n_filtered, (double)d->n_filtered / d->n); } |