Commit 8ae23835 authored by Thomas White's avatar Thomas White
Browse files

partialator: Fail a bit more gracefully if equations can't be solved

parent d541875f
......@@ -263,6 +263,10 @@ static double iterate_scale(struct image *images, int n,
e_val = gsl_vector_alloc(n);
e_vec = gsl_matrix_alloc(n, n);
val = gsl_eigen_symmv(M, e_val, e_vec, work);
if ( val ) {
ERROR("Couldn't diagonalise matrix.\n");
return 0.0;
}
gsl_eigen_symmv_free(work);
val = gsl_eigen_symmv_sort(e_val, e_vec, GSL_EIGEN_SORT_ABS_DESC);
......
......@@ -22,6 +22,7 @@
#include <getopt.h>
#include <assert.h>
#include <pthread.h>
#include <gsl/gsl_errno.h>
#include "utils.h"
#include "hdf5-file.h"
......@@ -274,6 +275,8 @@ int main(int argc, char *argv[])
}
STATUS("There are %i patterns to process\n", n_total_patterns);
gsl_set_error_handler_off();
images = malloc(n_total_patterns * sizeof(struct image));
if ( images == NULL ) {
ERROR("Couldn't allocate memory for images.\n");
......
......@@ -320,13 +320,17 @@ static double pr_iterate(struct image *image, const RefList *full,
//show_matrix_eqn(M, v, NUM_PARAMS);
shifts = gsl_vector_alloc(NUM_PARAMS);
gsl_linalg_HH_solve(M, v, shifts);
max_shift = 0.0;
for ( param=0; param<NUM_PARAMS; param++ ) {
double shift = gsl_vector_get(shifts, param);
apply_shift(image, param, shift);
if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
if ( gsl_linalg_HH_solve(M, v, shifts) ) {
ERROR("Couldn't solve normal equations!\n");
} else {
for ( param=0; param<NUM_PARAMS; param++ ) {
double shift = gsl_vector_get(shifts, param);
apply_shift(image, param, shift);
if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
}
}
gsl_matrix_free(M);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment