raising errors if singular matrices

sm/lib/gpc/gpc.c

@@ -20,11 +20,25 @@
 #include <gsl/gsl_math.h>
 #include "gpc.h"
 #include "gpc_utils.h"
+#include "../../csm/logging.h"
 
 /* Note that we use static values here so we don't need to evaluate that all the time */
 #define M(matrix, rows, col) static gsl_matrix*matrix=0; if(!matrix) matrix = gsl_matrix_alloc(rows,col);
 #define MF(matrix) /*gsl_matrix_free(matrix)*/
 
+static int
+singular (const gsl_matrix * LU)
+{
+  size_t i, n = LU->size1;
+
+  for (i = 0; i < n; i++)
+    {
+      double u = gsl_matrix_get (LU, i, i);
+      if (u == 0) return 1;
+    }
+
+ return 0;
+}
 
 int gpc_solve(int K, const struct gpc_corr*c,
 	const double*x0, const double *cov_x0,
@@ -131,6 +145,11 @@ int gpc_solve(int K, const struct gpc_corr*c,
 	 
 	/*	mS = mD - mB.trans * mA.inv * mB;
 	  temp22b = inv(A) */
+	if ( singular(mA) )
+	{
+	    sm_error("gpc_solve: mA is singular!\n");
+	    return 0;
+	}
 	m_inv(mA, temp22b); 
 	/* temp22c = inv(A) * mB           */
 	m_mult(temp22b, mB, temp22c);
@@ -141,6 +160,11 @@ int gpc_solve(int K, const struct gpc_corr*c,
 	m_add(mD,temp22b,mS);
 	
 	/* mSa = mS.inv * mS.det; */
+    if ( singular(mS) )
+    {
+        sm_error("gpc_solve: mS is singular!\n");
+        return 0;
+    }
 	m_inv(mS, mSa);
 	m_scale(m_det(mS), mSa);
 	
@@ -163,6 +187,11 @@ int gpc_solve(int K, const struct gpc_corr*c,
 	m_trans(g1, g1t);
 	m_trans(g2, g2t);
 	m_trans(mB, mBt);
+    if ( singular(mA) )
+    {
+        sm_error("gpc_solve: mA is singular!\n");
+        return 0;
+    }
 	m_inv(mA, mAi);
 
 	M(m1t,1,2);	M(m1,2,1);
@@ -255,6 +284,11 @@ int gpc_solve(int K, const struct gpc_corr*c,
 
 	m_scale(2*lambda, W);
 	gsl_matrix_add(bigM,W);
+    if ( singular(bigM) )
+    {
+        sm_error("gpc_solve: bigM is singular!\n");
+        return 0;
+    }
 	m_inv(bigM, temp44);
 	m_mult(temp44, g, x);
 	m_scale(-1.0, x);