Add public version of tm_cond()

zvezdochiot · Jul 19, 2022 · cc01817 · cc01817
1 parent 77f68f5
commit cc01817
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 3 deletions.
diff --git a/src/lib/tmatrix.h b/src/lib/tmatrix.h
@@ -299,5 +299,12 @@ int tm_rank(tMat *m, int *err);
  * @return Square norm.
  */
 tmVal tm_norm2(tMat* m, int* err);
+/**
+ * @brief Find condition number of the matrix.
+ * @param m matrix object.
+ * @param err error code.
+ * @return condition number.
+ */
+tmVal tm_cond(tMat* m, int* err);
 
 #endif /* T_MATRIX_H */
diff --git a/src/lib/tmatrix_calc.c b/src/lib/tmatrix_calc.c
@@ -14,8 +14,8 @@
 #define PINV_TOL     1E-9
 #define MEM_TMP_VEC 8
 
-/* Condition number */
-tmVal tm_cond(tMat *m, tMat *minv, int *err)
+/* Condition number (internal) */
+tmVal tm_cond_(tMat *m, tMat *minv, int *err)
 {
   int e = 0, rA;
   tmVal cond = 0.0;
@@ -319,7 +319,7 @@ tmVal tm_inv(tMat *dst, tMat *m, int *err)
               lubksb(&tmp,idx,&col);
             }
             // check condition number
-            cond = tm_cond(m, dst, &e);
+            cond = tm_cond_(m, dst, &e);
             if(!e && !(cond > 0.0))
               e = TM_ERR_NO_SOLUTN;
           }
@@ -338,6 +338,25 @@ tmVal tm_inv(tMat *dst, tMat *m, int *err)
   return cond;
 }
 
+/* Find condition number */
+tmVal tm_cond(tMat* m, int* err)
+{
+  int e = 0;
+  tmVal cond = 0.0;
+  /* allocate memory */
+  tMat minv = tm_copy(m, &e); 
+
+  if(!e) {
+    /* condition number is calculated inside the tm_inv function */
+    cond = tm_inv(&minv, m, &e);
+  }
+
+  if(err) *err = e;
+  tm_clear(&minv);
+
+  return cond;
+}
+
 /* Pseudoinverse aux 1 */
 tMat pinva(tMat *src, int* transp, tmVal* tolerance, int* err)
 {

diff --git a/src/tests/tests.c b/src/tests/tests.c
@@ -311,6 +311,9 @@ static char* test_inv()
   cond = tm_inv(&im, &m, &err);
   mu_check("Inv (det):", err);
   mu_assert("Inv (inv): bad condition number", cond > 0 && cond < 1E14);
+
+  mu_assert("Inv (cond): wrong condition number", EQL(cond, tm_cond(&m, &err)));
+  mu_check("Inv (cond):", err);
   pr = tm_simp();
 
   tm_mul(&pr,&m,&im,&err);