Merge pull request #22 from j-danner/master

fix sanity check for underdetermined systems + add test case

Author: malb

m4ri/solve.c

@@ -32,9 +32,9 @@ int mzd_solve_left(mzd_t *A, mzd_t *B, int const cutoff, int const inconsistency
     m4ri_die("mzd_solve_left: A ncols (%d) must be smaller than B nrows (%d).\n", A->ncols,
              B->nrows);
 
-  if (A->nrows != B->nrows)
-    m4ri_die("mzd_solve_left: A nrows (%d) must be equal to B nrows (%d).\n", A->nrows,
-             B->nrows);
+  if (B->nrows != MAX(A->ncols,A->nrows))
+    m4ri_die("mzd_solve_left: B nrows (%d) must be equal to max of A nrows (%d) and A ncols (%d).\n", B->nrows,
+             A->nrows, A->ncols);
 
   return _mzd_solve_left(A, B, cutoff, inconsistency_check);
 }
@@ -120,6 +120,12 @@ int _mzd_pluq_solve_left(mzd_t const *A, rci_t rank, mzp_t const *P, mzp_t const
 }
 
 int _mzd_solve_left(mzd_t *A, mzd_t *B, int const cutoff, int const inconsistency_check) {
+  if (inconsistency_check && B->nrows > A->nrows) {
+    mzd_t const *Bpad = mzd_init_window_const(B, A->nrows+1, 0, B->nrows, B->ncols);
+    if(!mzd_is_zero(Bpad)) return -1;
+    mzd_free_window((mzd_t *) Bpad);
+  }
+
   /**
    *  B is modified in place
    *  (Bi's in the comments are just modified versions of B)

m4ri/solve.h

@@ -36,6 +36,8 @@
  * \brief Solves A X = B with A and B matrices.
  *
  * The solution X is stored inplace on B.
+ * If A->nrows < A->ncols, the matrix A is implicitly padded with zeros to
+ * match B->nrows.
  *
  * \param A Input matrix (overwritten).
  * \param B Input matrix, being overwritten by the solution matrix X

tests/test_solve.c

@@ -59,19 +59,39 @@ int test_solve_left_random(rci_t mA, rci_t nA, rci_t nB, int consistent) {
   mzd_t *secret = mzd_init(nA, nB);
   mzd_randomize(A);
   mzd_randomize(secret);
-  mzd_t *B; 
-  if (consistent) {
-    B = mzd_mul(NULL, A, secret, 0); 
+  mzd_t *B;
+  if(consistent) {
+    B = mzd_mul(NULL, A, secret, 0);
     assert(B->nrows == A->nrows);
     assert(B->ncols == nB);
   }
   else {
-    B = mzd_init(mB, nB);
+    B = mzd_init(mA, nB);
     mzd_randomize(B);
   }
+  mzd_free(secret);
+
   // copy A & B
   mzd_t *Acopy = mzd_copy(NULL, A);
   mzd_t *Bcopy = mzd_copy(NULL, B);
+
+  // add rows to B, s.t. X fits into B
+  if(B->nrows < mB) {
+    mzd_t *padding = mzd_init(nA-mA, nB);
+    if(!consistent)
+      mzd_randomize(padding);
+    mzd_t *B_padded = mzd_stack(NULL, B, padding);
+    mzd_free(B);
+    B = B_padded;
+    mzd_free(padding);
+
+    mzd_t const *Bpad = mzd_init_window_const(B, A->nrows+1, 0, B->nrows, B->ncols);
+    mzd_print(Bpad);
+    mzd_free_window((mzd_t *) Bpad);
+  }
+  assert(B->nrows == mB);
+  assert(B->ncols == nB);
+
   int consistency = !mzd_solve_left(A, B, 0, 1);
 
   if (consistent && !consistency) {
@@ -83,8 +103,8 @@ int test_solve_left_random(rci_t mA, rci_t nA, rci_t nB, int consistent) {
     printf("skipped (OK, no solution found)\n");
     return 0;
   }
-  // copy B
-  mzd_t *X = mzd_submatrix(NULL, B, 0, 0, A->ncols, B->ncols);
+
+  mzd_t const *X = mzd_init_window_const(B, 0, 0, A->ncols, B->ncols);
   mzd_t *B1 = mzd_mul(NULL, Acopy, X, 0);
   mzd_t *Z  = mzd_add(NULL, Bcopy, B1);
 
@@ -99,8 +119,8 @@ int test_solve_left_random(rci_t mA, rci_t nA, rci_t nB, int consistent) {
   mzd_free(B1);
   mzd_free(Z);
   mzd_free(A);
+  mzd_free((mzd_t*) X);
   mzd_free(B);
-  mzd_free(X);
   return status;
 }
 
@@ -114,12 +134,14 @@ int main() {
   status += test_solve_left_random(1100, 1100, 1000, 1);
   status += test_solve_left_random(1000, 1000, 1100, 1);
   status += test_solve_left_random(1100, 1000, 1100, 1);
+  status += test_solve_left_random(1000, 1100, 1000, 1);
 
   status += test_solve_left_random(1100, 1000, 1000, 0);
   status += test_solve_left_random(1000, 1000, 1000, 0);
   status += test_solve_left_random(1100, 1100, 1000, 0);
   status += test_solve_left_random(1000, 1000, 1100, 0);
   status += test_solve_left_random(1100, 1000, 1100, 0);
+  status += test_solve_left_random(1000, 1100, 1000, 0);
 
 
   for (size_t i = 0; i < 100; i++) {