debuging singular matrix

sm/lib/gpc/gpc.c

@@ -26,20 +26,6 @@
 #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,
 	double *x_out) 

sm/lib/gpc/gpc_utils.c

@@ -30,6 +30,11 @@ void m_inv(const gsl_matrix*A, gsl_matrix*invA) {
 	/* Make LU decomposition of matrix m */
 	int s;
 	gsl_linalg_LU_decomp (m, perm, &s);
+	if ( singular (m))
+	{
+	    fprintf(stderr, "m_inv: m is singular!\n");
+	    return;
+	}
 	/* Invert the matrix m */
 	gsl_linalg_LU_invert (m, perm, invA);
 	gsl_permutation_free(perm);
@@ -132,3 +137,15 @@ void m_display(const char*str, gsl_matrix*m) {
 	}
 }
 
+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;
+}

sm/lib/gpc/gpc_utils.h

@@ -28,4 +28,6 @@ int poly_greatest_real_root(unsigned int n, const double*a, double *root);
 
 void m_display(const char*str, gsl_matrix*m);
 
+int singular (const gsl_matrix * LU);
+
 #endif