Commits

Martin Albrecht committed 2ffea05

fixing solving (for systems which are consistent)

Comments (0)

Files changed (5)

src/packedmatrix.h

  *
  * \param DST May be NULL for automatic creation.
  * \param A Source matrix.
- *
- * \wordoffset
  */
 
 mzd_t *mzd_copy(mzd_t *DST, const mzd_t *A);
 #include "trsm.h"
 #include "permutation.h"
 
-void mzd_solve_left(mzd_t *A, mzd_t *B, const int cutoff, const int inconsistency_check) {    
+int mzd_solve_left(mzd_t *A, mzd_t *B, const int cutoff, const int inconsistency_check) {    
   if(A->ncols > B->nrows)
     m4ri_die("mzd_solve_left: A ncols (%d) need to be lower than B nrows (%d).\n", A->ncols, B->nrows);
 
-  _mzd_solve_left (A, B, cutoff, inconsistency_check);
+  return _mzd_solve_left (A, B, cutoff, inconsistency_check);
 }
  
-void mzd_pluq_solve_left (mzd_t *A, size_t rank, 
-                          mzp_t *P, mzp_t *Q, 
-                          mzd_t *B, const int cutoff, const int inconsistency_check) 
-{
+int mzd_pluq_solve_left (mzd_t *A, size_t rank, 
+                         mzp_t *P, mzp_t *Q, 
+                         mzd_t *B, const int cutoff, const int inconsistency_check) {
   if(A->ncols > B->nrows)
     m4ri_die("mzd_pluq_solve_left: A ncols (%d) need to be lower than B nrows (%d).\n", A->ncols, B->nrows);
   if(P->length != A->nrows)
   if(Q->length != A->ncols)
       m4ri_die("mzd_pluq_solve_left: A ncols (%d) need to match Q size (%d).\n", A->ncols, P->length);
 
-  _mzd_pluq_solve_left (A, rank, P, Q, B, cutoff, inconsistency_check);
+  return _mzd_pluq_solve_left (A, rank, P, Q, B, cutoff, inconsistency_check);
 }
 
-void _mzd_pluq_solve_left (mzd_t *A, size_t rank, 
-                           mzp_t *P, mzp_t *Q, 
-                           mzd_t *B, const int cutoff, const int inconsistency_check) {
+int _mzd_pluq_solve_left (mzd_t *A, size_t rank, 
+                          mzp_t *P, mzp_t *Q, 
+                          mzd_t *B, const int cutoff, const int inconsistency_check) {
   /** A is supposed to store L lower triangular and U upper triangular
    *  B is modified in place 
    *  (Bi's in the comments are just modified versions of B)
    *  4) Q B5 = B4
    */
 
+  int retval = 0;
+
   /* P B2 = B1 or B2 = P^T B1*/
   mzd_apply_p_left(B, P);
   
   mzd_t *Y1 = mzd_init_window(B,0,0,rank,B->ncols);
   mzd_trsm_lower_left(LU, Y1, cutoff);
   
-  if (inconsistency_check) {
-    /* Check for inconsistency */
-    /** FASTER without this check
-     * 
-     * update with the lower part of L 
+  if (inconsistency_check) { /* Check for inconsistency */    
+    /** FASTER without this check; update with the lower part of L
      */
-    mzd_t *H = mzd_init_window(A, rank, 0, A->nrows, rank);
-    mzd_t *Y2 = mzd_init_window(B,rank,0,B->nrows,B->ncols);
+    mzd_t *H  = mzd_init_window(A, rank, 0, A->nrows, rank);
+    mzd_t *Y2 = mzd_init_window(B, rank, 0, A->nrows, B->ncols);
+    if(A->nrows < B->nrows) {
+      mzd_t *Y3 = mzd_init_window(B, A->nrows, 0, B->nrows, B->ncols);
+      mzd_set_ui(Y3, 0);
+      mzd_free_window(Y3);
+    }
+
     mzd_addmul(Y2, H, Y1, cutoff);
     /*
      * test whether Y2 is the zero matrix
      */
     if( !mzd_is_zero(Y2) ) {
-      //printf("inconsistent system of size %llu x %llu\n", Y2->nrows, Y2->ncols);
-      //printf("Y2=");
-      //mzd_print(Y2);
+      retval = -1;
     }
     mzd_free_window(H);
     mzd_free_window(Y2);
   mzd_free_window(Y1);
   
   if (!inconsistency_check) {
-    /** Default is to set the indefined bits to zero 
-     * if inconsistency has been checked then 
-     *    Y2 bits are already all zeroes
-     * thus this clearing is not needed
+    /** Default is to set the undefined bits to zero if inconsistency
+     * has been checked then Y2 bits are already all zeroes thus this
+     * clearing is not needed
      */
     for(size_t i = rank; i<B->nrows; i++) {
       for(size_t j=0; j<B->ncols; j+=RADIX) {
     }
   }
   /* Q B5 = B4 or B5 = Q^T B4*/
-  mzd_apply_p_right(B, Q);
+  mzd_apply_p_left_trans(B, Q);
 
   /* P L U Q B5 = B1 */
+  return retval;
 }
 
-void _mzd_solve_left (mzd_t *A, mzd_t *B, const int cutoff, const int inconsistency_check) {
+int _mzd_solve_left (mzd_t *A, mzd_t *B, const int cutoff, const int inconsistency_check) {
   /**
    *  B is modified in place 
    *  (Bi's in the comments are just modified versions of B)
   /* PLUQ = A */
   size_t rank = _mzd_pluq(A, P, Q, cutoff);  
   /* 2, 3, 4, 5 */
-  mzd_pluq_solve_left(A, rank, P, Q, B, cutoff, inconsistency_check);
+  int retval = mzd_pluq_solve_left(A, rank, P, Q, B, cutoff, inconsistency_check);
   
   mzp_free(P);
   mzp_free(Q);
+  return retval;
 }
 
 mzd_t *mzd_kernel_left_pluq(mzd_t *A, const int cutoff) {
  * \param A Input matrix (overwritten).
  * \param B Input matrix, being overwritten by the solution matrix X
  * \param cutoff Minimal dimension for Strassen recursion (default: 0).
- * \param inconsistency_check decide wether or not to check for
- *        incosistency (faster without but output not defined if
- *        system is not consistent).
+ * \param inconsistency_check decide wether or not to perform a quick
+ *        check for incosistency (faster without but output not
+ *        defined if system is not consistent).
+ * \return 0 if a solution was found, -1 otherwise
  *
+ * \warn This function may return wrong solutions if no solution
+ * exists even with inconsistency_check set to 1. To check the result
+ * perform a matrix multiplication.
  */
-void mzd_solve_left(mzd_t *A, mzd_t *B, const int cutoff, 
-                    const int inconsistency_check);
+int mzd_solve_left(mzd_t *A, mzd_t *B, const int cutoff, const int inconsistency_check);
 
 /**
  * \brief Solves (P L U Q) X = B
  * \param Q Input column permutation matrix.
  * \param B Input matrix, being overwritten by the solution matrix X.
  * \param cutoff Minimal dimension for Strassen recursion (default: 0).
- * \param inconsistency_check decide whether or not to check for
- *        incosistency (faster without but output not defined if
- *        system is not consistent).
+ * \param inconsistency_check decide whether or not to perform a quick
+ *        check for incosistency (faster without but output not
+ *        defined if system is not consistent).  \return 0 if a
+ *        solution was found, -1 otherwise
+ * \return 0 if a solution was found, -1 otherwise
+ *
+ * \warn This function may return wrong solutions if no solution
+ * exists even with inconsistency_check set to 1. To check the result
+ * perform a matrix multiplication.
  */
-void mzd_pluq_solve_left (mzd_t *A, size_t rank, 
-                          mzp_t *P, mzp_t *Q, 
-                          mzd_t *B, const int cutoff, const int inconsistency_check);
+int mzd_pluq_solve_left (mzd_t *A, size_t rank, 
+                         mzp_t *P, mzp_t *Q, 
+                         mzd_t *B, const int cutoff, const int inconsistency_check);
 
 /**
  * \brief  Solves (P L U Q) X = B
  * \param Q Input column permutation matrix.
  * \param B Input matrix, being overwritten by the solution matrix X.
  * \param cutoff Minimal dimension for Strassen recursion (default: 0).
- * \param inconsistency_check decide whether or not to check for
- *        incosistency (faster without but output not defined if
- *        system is not consistent).
- *
+ * \param inconsistency_check decide whether or not to perform a quick
+ *        check for incosistency (faster without but output not
+ *        defined if system is not consistent).  \return 0 if a
+ *        solution was found, -1 otherwise
+ * \return 0 if a solution was found, -1 otherwise
  */
-void _mzd_pluq_solve_left(mzd_t *A, size_t rank, 
-                          mzp_t *P, mzp_t *Q, 
-                          mzd_t *B, const int cutoff, const int inconsistency_check);
+int _mzd_pluq_solve_left(mzd_t *A, size_t rank, 
+                         mzp_t *P, mzp_t *Q, 
+                         mzd_t *B, const int cutoff, const int inconsistency_check);
 
 /**
  * \brief Solves A X = B with A and B matrices.
  * \param A Input matrix.
  * \param B Input matrix, being overwritten by the solution matrix X.
  * \param cutoff Minimal dimension for Strassen recursion (default: 0).
- * \param inconsistency_check decide whether or not to check for
- *        incosistency (faster without but output not defined if
- *        system is not consistent).
+ * \param inconsistency_check decide whether or not to perform a quick
+ *        check for incosistency (faster without but output not
+ *        defined if system is not consistent).  \return 0 if a
+ *        solution was found, -1 otherwise
+ * \return 0 if a solution was found, -1 otherwise
+ *
+ * \warn This function may return wrong solutions if no solution
+ * exists even with inconsistency_check set to 1. To check the result
+ * perform a matrix multiplication.
  */
-void _mzd_solve_left(mzd_t *A, mzd_t *B, const int cutoff, const int inconsistency_check);
+int _mzd_solve_left(mzd_t *A, mzd_t *B, const int cutoff, const int inconsistency_check);
 
 
 /**
  *
  * \sa mzd_pluq()
  *
- * \return X
+ * \return X, NULL if kernel is empty
  */
 
 mzd_t *mzd_kernel_left_pluq(mzd_t *A, const int cutoff);

testsuite/test_pluq.c

 
   mzd_t* Acopy = mzd_copy (NULL,A);
 
-
-
   mzp_t* Pt = mzp_init(m);
   mzp_t* Q = mzp_init(n);
   int r = mzd_pluq(A, Pt, Q, 0);

testsuite/test_solve.c

 
   mzd_t *Acopy = mzd_copy(NULL, A);
   size_t r = mzd_echelonize_m4ri(Acopy,0,0);
-  printf("solve_left m: %4zu, n: %4zu, r: %4zu da: %2zu db: %2zu ",m, n, r, offsetA, offsetB);
+  printf("solve_left m: %4zu, n: %4zu, r: %4zu da: %4zu db: %4zu ",m, n, r, offsetA, offsetB);
   mzd_free(Acopy);
   Acopy = mzd_copy(NULL, A);
     
-  mzd_solve_left(A, B, 0, 1);
+  int consistency = mzd_solve_left(A, B, 0, 1);
 
   //copy B
   mzd_t *X = mzd_init(B->nrows,B->ncols);
   
   int status = 0;
   
-  if(r==m) {
-    for ( i=0; i<m; ++i)
-      for ( j=0; j<n; ++j){
-        if (mzd_read_bit (Z,i,j)){
-          status = 1;
-        }
-      }
-    if (!status)
+  if(consistency == 0) {
+    status = 1 - mzd_is_zero(Z);
+    if (status == 0) {
       printf("passed\n");
-    else
+    } else if (r != n) {
+      printf("failed (but expected)\n");
+      status = 0;
+    } else {
       printf("FAILED\n");
-
+    }
   } else {
-    printf("check skipped r!=m\n");
+    printf("skipped (no solution)\n");
   }
   mzd_free(Bcopy);
   mzd_free(B1);
   status += test_pluq_solve_left( 80,  20,  0,  0);
   status += test_pluq_solve_left( 80,  80,  0,  0);
 
+  status += test_pluq_solve_left (10, 20, 0, 0);
+  status += test_pluq_solve_left (10, 80, 0, 0);
+  status += test_pluq_solve_left (10, 20, 0, 0);
+  status += test_pluq_solve_left (10, 80, 0, 0);
+  status += test_pluq_solve_left (70, 20, 0, 0);
+  status += test_pluq_solve_left (70, 80, 0, 0);
+  status += test_pluq_solve_left (70, 20, 0, 0);
+  status += test_pluq_solve_left (70, 80, 0, 0);
+  status += test_pluq_solve_left (770, 1600, 0, 0);
+  status += test_pluq_solve_left (1764, 1345, 0, 0);
+
   status += test_pluq_solve_left(  2,   4,  0,  1);
   status += test_pluq_solve_left(  4,   1,  0,  1);
   status += test_pluq_solve_left( 10,  20,  0,  1);
   status += test_pluq_solve_left( 80,  20,  0, 15);
   status += test_pluq_solve_left( 80,  80,  0, 15);
 
-/*   status += test_pluq_solve_left (10, 20, 15, 0); */
-/*   status += test_pluq_solve_left (10, 80, 15, 0); */
-/*   status += test_pluq_solve_left (10, 20, 15, 20); */
-/*   status += test_pluq_solve_left (10, 80, 15, 20); */
-/*   status += test_pluq_solve_left (70, 20, 15, 0); */
-/*   status += test_pluq_solve_left (70, 80, 15, 0); */
-/*   status += test_pluq_solve_left (70, 20, 15, 20); */
-/*   status += test_pluq_solve_left (70, 80, 15, 20); */
-/*   status += test_pluq_solve_left (770, 1600, 75, 89); */
-/*   status += test_pluq_solve_left (1764, 1345, 198, 123); */
+  status += test_pluq_solve_left (10, 20, 15, 0);
+  status += test_pluq_solve_left (10, 80, 15, 0);
+  status += test_pluq_solve_left (10, 20, 15, 20);
+  status += test_pluq_solve_left (10, 80, 15, 20);
+  status += test_pluq_solve_left (70, 20, 15, 0);
+  status += test_pluq_solve_left (70, 80, 15, 0);
+  status += test_pluq_solve_left (70, 20, 15, 20);
+  status += test_pluq_solve_left (70, 80, 15, 20);
+  status += test_pluq_solve_left (770, 1600, 75, 89);
+  status += test_pluq_solve_left (1764, 1345, 198, 123);
 
   if (!status) {
     printf("All tests passed.\n");