From a18415b78798e53c50c08bdefe83db4a377aa461 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Tue, 5 Mar 2019 13:16:52 -0600 Subject: update quad() to not abort if the matrix is singular --- src/quad.c | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) (limited to 'src/quad.c') diff --git a/src/quad.c b/src/quad.c index c869aed..848a12e 100644 --- a/src/quad.c +++ b/src/quad.c @@ -14,9 +14,16 @@ #include #include #include "misc.h" +#include char quad_err[256]; +void my_handler(const char *reason, const char *file, int line, int gsl_errno) +{ + fprintf(stderr, "gsl: %s:%d: %s: %s\n", file, line, "ERROR", reason); + return; +} + /* Returns an approximate vertex position and time using the QUAD fitter. * * The QUAD fitter was originally developed for SNO and estimates the event @@ -83,6 +90,8 @@ int quad(event *ev, double *pos, double *t0, size_t npoints) static size_t p2[MAX_PMTS]; double tmin; size_t max_tries; + gsl_error_handler_t *old_handler; + int status; wavelength0 = 400.0; n_d2o = get_index_snoman_d2o(wavelength0); @@ -156,7 +165,22 @@ int quad(event *ev, double *pos, double *t0, size_t npoints) gsl_vector_view n_view = gsl_vector_view_array(n,3); gsl_linalg_LU_decomp(&m.matrix, p, &s); - gsl_linalg_LU_invert(&m.matrix, p, &minv.matrix); + + /* Ocassionaly the matrix is singular and we can't invert it. Since we + * don't want our program to quit when this happens, we install a new + * gsl error handler, do the matrix inversion, and then restore the old + * handler. If the matrix inversion failed, then we just continue with + * the next quad point. */ + + /* save original handler, install new handler */ + old_handler = gsl_set_error_handler(&my_handler); + + status = gsl_linalg_LU_invert(&m.matrix, p, &minv.matrix); + + /* restore original handler */ + gsl_set_error_handler(old_handler); + + if (status) continue; gsl_blas_dgemv(CblasNoTrans,1.0,&minv.matrix,&k_view.vector,0.0,&g_view.vector); gsl_blas_dgemv(CblasNoTrans,1.0,&minv.matrix,&n_view.vector,0.0,&h_view.vector); -- cgit