Commits

Martin Albrecht  committed dca9c55

improved documentation (added docs on return values) and removed redundant parameters from mzd_echelonize_m4ri()

  • Participants
  • Parent commits 592f751

Comments (0)

Files changed (8)

File src/brilliantrussian.c

   }
 }
 
-int mzd_echelonize_m4ri(packedmatrix *A, int full, int k, packedmatrix *T, size_t *L) {
+size_t mzd_echelonize_m4ri(packedmatrix *A, int full, int k) {
   /**
    * \par General algorithm
    * \li Step 1.Denote the first column to be processed in a given
   return r;
 }
 
-void mzd_top_echelonize_m4ri(packedmatrix *A, int k, packedmatrix *T, size_t *L) {
+void mzd_top_echelonize_m4ri(packedmatrix *A, int k) {
   const size_t ncols = A->ncols; 
   size_t r = 0;
   size_t c = 0;
   size_t *L=(size_t *)m4ri_mm_malloc(twokay * sizeof(size_t));
   packedmatrix *answer;
   
-  mzd_echelonize_m4ri(big, TRUE, k, T, L);
+  mzd_echelonize_m4ri(big, TRUE, k);
   
   for(i=0; i < size; i++) {
     if (!mzd_read_bit(big, i,i )) {

File src/brilliantrussian.h

 void mzd_make_table( packedmatrix *M, size_t r, size_t c, int k, packedmatrix *T, size_t *L);
 
 /**
- * \brief The function looks up k bits from position i,startcol in each row
- * and adds the appropriate row from T to the row i. 
+ * \brief The function looks up k bits from position i,startcol in
+ * each row and adds the appropriate row from T to the row i.
  *
  * This process is iterated for i from startrow to stoprow
  * (exclusive).
 /**
  * \brief Matrix elimination using the 'Method of the Four Russians'
  * (M4RI).
- * 
+ *
  * \param M Matrix to be reduced.
  * \param full Return the reduced row echelon form, not only upper triangular form.
  * \param k M4RI parameter, may be 0 for auto-choose.
- * \param T Preallocated table, may be NULL for automatic creation.
- * \param L Preallocated lookup table, may be NULL for automatic creation.
  *
  * \example testsuite/test_elimination.c
  * \example testsuite/bench_elimination.c
  * 
  * \wordoffset
+ *
+ * \return Rank of A.
  */
 
-int mzd_echelonize_m4ri(packedmatrix *M, int full, int k, packedmatrix *T, size_t *L);
+size_t mzd_echelonize_m4ri(packedmatrix *M, int full, int k);
 
 /**
  * \brief Given a matrix in upper triangular form compute the reduced row
  * \param L Preallocated lookup table, may be NULL for automatic creation.
  *
  * \wordoffset
+ *
+ * \note This function isn't as optimized as it should be.
  */
 
-void mzd_top_echelonize_m4ri(packedmatrix *M, int k, packedmatrix *T, size_t *L);
+void mzd_top_echelonize_m4ri(packedmatrix *M, int k);
 
 /**
- * \brief Invert the matrix M using Konrod's method. To avoid
- * recomputing the identity matrix over and over again, I may be
- * passed in as identity parameter.
+ * \brief Invert the matrix M using Konrod's method. 
+ * 
+ * To avoid recomputing the identity matrix over and over again, I may
+ * be passed in as identity parameter.
  *
  * \param M Matrix to be reduced.
  * \param I Identity matrix.
  * \param k M4RI parameter, may be 0 for auto-choose.
  *
  * \wordoffset
+ *
+ * \return Inverse of M
+ *
+ * \note This function allocates a new matrix for the inverse which
+ * must be free'd using mzd_free() once not needed anymore.
  */
 
 packedmatrix *mzd_invert_m4ri(packedmatrix *M, packedmatrix *I, int k);
  * \param k M4RI parameter, may be 0 for auto-choose.
  *
  * \wordoffset
+ *
+ * \return Pointer to C.
  */
 
 packedmatrix *mzd_mul_m4rm(packedmatrix *C, packedmatrix *A, packedmatrix *B, int k);
  * \param k M4RI parameter, may be 0 for auto-choose.
  *
  * \wordoffset
+ *
+ * \return Pointer to C.
  */
 
 packedmatrix *mzd_addmul_m4rm(packedmatrix *C, packedmatrix *A, packedmatrix *B, int k);
  * \author William Hart -- block matrix implementation, use of several Gray code tables, general speed-ups
  *
  * \wordoffset
+ *
+ * \return Pointer to C.
  */
 
 packedmatrix *_mzd_mul_m4rm(packedmatrix *C, packedmatrix *A, packedmatrix *B, int k, int clear);
  * \brief PLUQ matrix decomposition routines.
  *
  * \author Clement Pernet <clement.pernet@gmail.com>
+ * 
+ * \note This file should be called pluq.h and will be renamed in the
+ * future.
  */
 
 
  *
  * If (P,L,U,Q) satisfy PLUQ = A, it returns (P^T, L, U, Q^T).
  *
+ * P and Q must be preallocated but they don't have to be
+ * identity permutations. If cutoff is zero a value is chosen
+ * automatically. It is recommended to set cutoff to zero for most
+ * applications.
+ *
  * The row echelon form (not reduced) can be read from the upper
  * triangular matrix U. See mzd_echelonize_pluq() for details.
  * 
  * This is the wrapper function including bounds checks. See
  * _mzd_pluq() for implementation details.
  *
- * \param A Input matrix
- * \param P Output row permutation matrix
- * \param Q Output column permutation matrix
+ * \param A Input m x n matrix
+ * \param P Output row permutation of length m
+ * \param Q Output column permutation matrix of length n
  * \param cutoff Minimal dimension for Strassen recursion.
  *
+ * \sa _mzd_pluq() _mzd_pluq_mmpf() mzd_echelonize_pluq()
+ *
+ * \example testsuite/test_lqup.c
+ *
+ * \wordoffset
+ * \return Rank of A.
  */
 
 size_t mzd_pluq(packedmatrix *A, permutation *P, permutation * Q, const int cutoff);
 /**
  * \brief PLUQ matrix decomposition.
  *
- * Computes the PLUQ matrix decomposition using a block recursive
- * algorithm.
+ * See mzd_pluq() for details.
  *
- * If (P,L,U,Q) satisfy PLUQ = A, this routine returns (P^T,L,U,Q^T).
- *
- * The row echelon form (not reduced) can be read from the upper
- * triangular matrix U*Q. See mzd_echelonize_pluq() for (reduced) row
- * echelon forms using PLUQ factorisation.
- * 
- * The matrix L and U are stored in place over A.
- * 
  * \param A Input matrix
  * \param P Output row permutation matrix
  * \param Q Output column permutation matrix
  * \param cutoff Minimal dimension for Strassen recursion.
+ *
+ * \sa mzd_pluq()
+ *
+ * \wordoffset
+ * \return Rank of A.
  */
 
 size_t _mzd_pluq(packedmatrix *A, permutation * P, permutation * Q, const int cutoff);
 /**
  * \brief PLUQ matrix decomposition (naive base case).
  *
- * Computes the PLUQ matrix decomposition using the naive algorithm.
- *
- * If (P,L,U,Q) satisfy PLUQ = A, it returns (P^T, L, U, Q^T). 
- * 
- * The matrix L and U are stored in place over A.
+ * See mzd_pluq() for details.
  * 
  * \param A Input matrix
  * \param P Output row permutation matrix
  * \param Q Output column permutation matrix
  *
  * \sa mzd_pluq()
+ *
+ * \wordoffset
+ * \return Rank of A.
  */
 
 size_t _mzd_pluq_naive(packedmatrix *A, permutation * P, permutation * Q);
  * \param full Return the reduced row echelon form, not only upper triangular form.
  *
  * \wordoffset
+ *
  * \sa mzd_pluq()
+ *
+ * \return Rank of A.
  */
 
 

File src/permutation.c

   m4ri_mm_free(condemned);
 }
 
+void mzp_set_ui(permutation *P, unsigned int value) {
+  size_t i;
+  for (i=0; i<P->length; i++) {
+    P->values[i] = i;
+  }
+}
+
 void mzd_apply_p_left(packedmatrix *A, permutation *P) {
   size_t i;
   if(A->ncols == 0)

File src/permutation.h

 void mzp_free_window(permutation* condemned);
 
 /**
+ * \brief Set the permutation P to the identity permutation. The only
+ * allowed value is 1.
+ *
+ *
+ * \param P Permutation
+ * \param value 1
+ *
+ * \note This interface was chosen to be similar to mzd_set_ui().
+ */
+
+void mzp_set_ui(permutation *P, unsigned int value);
+
+
+/**
  * Apply the permutation P to A from the left.
  *
  * This is equivalent to row swaps walking from length-1 to 0.

File src/pluq_mmpf.h

  * \param k Size of Gray code tables.
  *
  * \wordoffset
+ *
+ * \return Rank of A.
  */
 
 size_t _mzd_pluq_mmpf(packedmatrix *A, permutation * P, permutation * Q, int k);
  * \param P Preallocated row permutation.
  * \param Q Preallocated column permutation.
  * \param todo Preallocated temporary buffer.
+ *
+ * \retval kbar Maximum k for which a pivot could be found.
  */
 
 size_t _mzd_pluq_submatrix(packedmatrix *A, size_t r, size_t c, int k, permutation *P, permutation *Q, size_t *todo);

File testsuite/test_elimination.c

   packedmatrix *G = mzd_copy(NULL, A);
 
   /* M4RI k=auto */
-  mzd_echelonize_m4ri(A, 1, 0, NULL, NULL);
+  mzd_echelonize_m4ri(A, 1, 0);
 
   /* M4RI k=8 */
-  mzd_echelonize_m4ri(B, 1, 8, NULL, NULL);
+  mzd_echelonize_m4ri(B, 1, 8);
 
   /* M4RI Upper Triangular k=auto*/
-  mzd_echelonize_m4ri(C, 0, 0, NULL, NULL);
-  mzd_top_echelonize_m4ri(C, 0, NULL, NULL);
+  mzd_echelonize_m4ri(C, 0, 0);
+  mzd_top_echelonize_m4ri(C, 0);
 
   /* M4RI Upper Triangular k=4*/
-  mzd_echelonize_m4ri(D, 0, 4, NULL, NULL);
-  mzd_top_echelonize_m4ri(D, 4, NULL, NULL);
+  mzd_echelonize_m4ri(D, 0, 4);
+  mzd_top_echelonize_m4ri(D, 4);
 
   /* Gauss */
   mzd_echelonize_naive(E, 1);
 
   /* Gauss Upper Triangular */
   mzd_echelonize_naive(F, 0);
-  mzd_top_echelonize_m4ri(F, 0, NULL, NULL);
+  mzd_top_echelonize_m4ri(F, 0);
 
   /* PLUQ */
   mzd_echelonize_pluq(G, 1);

File testsuite/test_solve.c

   }
 
   packedmatrix *Acopy = mzd_copy(NULL, A);
-  size_t r = mzd_echelonize_m4ri(Acopy,0,0,NULL,NULL);
+  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);
   mzd_free(Acopy);
   Acopy = mzd_copy(NULL, A);