update quad() to not abort if the matrix is singular

1 files changed, 25 insertions, 1 deletions

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 <gsl/gsl_sort.h>
 #include <gsl/gsl_cblas.h>
 #include "misc.h"
+#include <gsl/gsl_errno.h>
 
 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);