From ece31b55aa87752f7054eee72ce036641d46356f Mon Sep 17 00:00:00 2001
From: joanvallve <jvallve@iri.upc.edu>
Date: Fri, 17 Dec 2021 17:05:38 +0100
Subject: [PATCH] raising errors if singular matrices
---
sm/lib/gpc/gpc.c | 34 ++++++++++++++++++++++++++++++++++
1 file changed, 34 insertions(+)
diff --git a/sm/lib/gpc/gpc.c b/sm/lib/gpc/gpc.c
index 97cb321..9ce35f8 100644
--- a/sm/lib/gpc/gpc.c
+++ b/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);
--
GitLab