# Commits

committed 9277f00

changed API and updated docs: mzd_reduce_ mzd_echelonize_

• Participants
• Parent commits 4b3b7a6
• Branches default

# File src/brilliantrussian.c

}
}

-int mzd_reduce_m4ri(packedmatrix *A, int full, int k, packedmatrix *T, size_t *L) {
+int mzd_echelonize_m4ri(packedmatrix *A, int full, int k, packedmatrix *T, size_t *L) {
/**
-   * The algorithm works as follows:
-   *
-   * Step 1.Denote the first column to be processed in a given
+   * \par General algorithm
+   * \li Step 1.Denote the first column to be processed in a given
* iteration as \f\$a_i\f\$. Then, perform Gaussian elimination on the
* first \f\$3k\f\$ rows after and including the \f\$i\f\$-th row to
* produce an identity matrix in \f\$a_{i,i} ... a_{i+k-1,i+k-1},\f\$
* and zeroes in \f\$a_{i+k,i} ... a_{i+3k-1,i+k-1}\f\$.
*
-   * Step 2. Construct a table consisting of the \f\$2^k\f\$ binary strings of
+   * \li Step 2. Construct a table consisting of the \f\$2^k\f\$ binary strings of
* length k in a Gray code.  Thus with only \f\$2^k\f\$ vector
* additions, all possible linear combinations of these k rows
* have been precomputed.
*
-   *
-   * Step 3. One can rapidly process the remaining rows from \f\$i +
+   * \li Step 3. One can rapidly process the remaining rows from \f\$i +
* 3k\f\$ until row \f\$m\f\$ (the last row) by using the table. For
* example, suppose the \f\$j\f\$-th row has entries \f\$a_{j,i}
* ... a_{j,i+k-1}\f\$ in the columns being processed. Selecting the
* remaining columns from \f\$ i + k\f\$ to n in the appropriate way,
* as if Gaussian elimination had been performed.
*
-   * Step 4. While the above form of the algorithm will reduce a
+   * \li Step 4. While the above form of the algorithm will reduce a
* system of boolean linear equations to unit upper triangular form,
* and thus permit a system to be solved with back substitution, the
* M4RI algorithm can also be used to invert a matrix, or put the
* system into reduced row echelon form (RREF). Simply run Step 3
* on rows \f\$0 ... i-1\f\$ as well as on rows \f\$i + 3k
* ... m\f\$. This only affects the complexity slightly, changing the
-   * 2.5 coeffcient to 3
+   * 2.5 coeffcient to 3.
+   *
+   * \attention This function implements a variant of the algorithm
+   * described above.
*/

const size_t ncols = A->ncols;
return r;
}

-void mzd_top_reduce_m4ri(packedmatrix *A, int k, packedmatrix *T, size_t *L) {
+void mzd_top_echelonize_m4ri(packedmatrix *A, int k, packedmatrix *T, size_t *L) {
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));

-  mzd_reduce_m4ri(big, TRUE, k, T, L);
+  mzd_echelonize_m4ri(big, TRUE, k, T, L);

for(i=0; i < size; i++) {
if (!mzd_read_bit(big, i,i )) {

# File src/brilliantrussian.h

/**
* \file brilliantrussian.h
- * \brief Matrix operations using Gray codes.
+ * \brief M4RI and M4RM.
*
* \author Gregory Bard <bard@fordham.edu>
* \author Martin Albrecht <M.R.Albrecht@rhul.ac.uk>
packedmatrix *T2, size_t *L2, packedmatrix *T3, size_t *L3);

/**
- * \brief Perform matrix reduction using the 'Method of the Four
- * Russians' (M4RI) or Kronrod-Method.
+ * \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 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
*/

-int mzd_reduce_m4ri(packedmatrix *M, int full, int k, packedmatrix *T, size_t *L);
+int mzd_echelonize_m4ri(packedmatrix *M, int full, int k, packedmatrix *T, size_t *L);

/**
* \brief Given a matrix in upper triangular form compute the reduced row
* \wordoffset
*/

-void mzd_top_reduce_m4ri(packedmatrix *M, int k, packedmatrix *T, size_t *L);
+void mzd_top_echelonize_m4ri(packedmatrix *M, int k, packedmatrix *T, size_t *L);

/**
* \brief Invert the matrix M using Konrod's method. To avoid

# File src/lqup.h

/**
* \file lqup.h
*
- * \brief PLUQ matrix decomposition routines
+ * \brief PLUQ matrix decomposition routines.
*
* \author Clement Pernet <clement.pernet@gmail.com>
- *
- * \internal
*/

* Computes the transposed PLUQ matrix decomposition using a block
* recursive algorithm.
*
- * If (P,L,U,Q) satisfy P^T LU Q^T = A, it returns (P, L, U, Q).
+ * If (P,L,U,Q) satisfy PLUQ = A, it returns (P, L, U, Q).
*
- * The Row echelon form (not reduced) can be read from the upper
- * triangular matrix U.
+ * 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.
+ * This is the wrapper function including bounds checks. See
+ * _mzd_pluq() for implementation details.
*
* The matrix L and U are stored in place over A.  U is represented
* by the matrix Q U Q^T
/**
* \brief PLUQ matrix decomposition.
*
- * Computes the transposed PLUQ matrix decomposition using a block
- * recursive algorithm.
+ * Computes the PLUQ matrix decomposition using a block recursive
+ * algorithm.
*
- * If (L,Q,U,P) satisfy PLUQ = A^T, it returns (L^T, Q^T, U^T, P^T).
+ * If (P,L,U,Q) satisfy PLUQ = A, this routine returns (P,L,U,Q).
*
- * The Row echelon form (not reduced) can be read from the upper
- * triangular matrix L^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.  L^T is represented
- * by the matrix Q^T L^T Q
+ * The matrix L and U are stored in place over A.
*
* \param A Input matrix
* \param P Output row permutation matrix
*
* Computes the PLUQ matrix decomposition using the naive algorithm.
*
- * If (L,Q,U,P) satisfy PLUQ = A, it returns (L, Q, U, P).
+ * If (P,L,U,Q) satisfy PLUQ = A, it returns (P, L, U, Q).
*
- * The Row echelon form (not reduced) can be read from the upper
- * triangular matrix L.
- *
- * The matrix L and U are stored in place over A.  L is represented by
- * the matrix Q L Q.
+ * 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
+ *
+ * \sa mzd_pluq()
*/

size_t _mzd_pluq_naive(packedmatrix *A, permutation * P, permutation * Q);

/**
- * \brief Compute the (reduced) row echelon form.
+ * \brief (Reduced) row echelon form using PLUQ factorisation.
*
* \param A Matrix.
* \param full Return the reduced row echelon form, not only upper triangular form.
*
* \wordoffset
+ * \sa mzd_pluq()
*/

# File src/m4ri.h

* http://m4ri.sagemath.org for details.
*
* \example testsuite/test_multiplication.c
- * \example testsuite/test_elimination.c
*/

#include <stdio.h>

# File src/misc.h

*
* \param n Integer
*
-* \warn Does not handle multiples of RADIX correctly
+* \warning Does not handle multiples of RADIX correctly
*/

/**
* \brief Return handle for locale memory management cache.
- *
- * \todo Make thread safe.
+ *
+ * \attention Not thread safe.
*/

static inline mm_block *m4ri_mmc_handle(void) {

# File src/packedmatrix.c

return pivots;
}

-int mzd_reduce_naive(packedmatrix *m, int full) {
+int mzd_echelonize_naive(packedmatrix *m, int full) {
return mzd_gauss_delayed(m, 0, full);
}

H = mzd_concat(NULL, A, I);

-  x = mzd_reduce_naive(H, TRUE);
+  x = mzd_echelonize_naive(H, TRUE);

if (x == FALSE) {
mzd_free(H);

# File src/packedmatrix.h

static inline void mzd_row_swap_offset(packedmatrix *M, const size_t rowa, const size_t rowb, const size_t offset) {
size_t i;
-  /** \todo: this is pathetic/test code **/
+  /**
+   * \todo this is pathetic/test code
+   **/
for(i=offset; i<M->ncols; i++) {
const BIT temp = mzd_read_bit(M, rowa, i);
mzd_write_bit(M, rowa, i, mzd_read_bit(M, rowb, i));
* destrow, but only begins at the column coloffset.
*
* \param M Matrix
- * \param sourcerow Index of source row
- * \param destrow Index of target row
+ * \param dstrow Index of target row
+ * \param srcrow Index of source row
* \param coloffset Column offset
*/

/**
* \brief Gaussian elimination.
*
- * This will do Gaussian elimination on the matrix m.  If  full =
- *  FALSE, then it will do triangular style elimination, and if
- *  full = TRUE, it will do Gauss-Jordan style, or full elimination.
+ * This will do Gaussian elimination on the matrix m.  If full=FALSE,
+ *  then it will do triangular style elimination, and if full=TRUE,
+ *  it will do Gauss-Jordan style, or full elimination.
*
* \param M Matrix
* \param full Gauss-Jordan style or upper triangular form only.
*
* \wordoffset
+ *
+ * \sa mzd_echelonize_m4ri(), mzd_echelonize_pluq()
*/

-int mzd_reduce_naive(packedmatrix *M, const int full);
+int mzd_echelonize_naive(packedmatrix *M, const int full);

/**
* \brief Return TRUE if A == B.

# File src/permutation.h

/**
* \file permutation.h
*
- * \brief Permutation matrices for the M4RI library.
+ * \brief Permutation matrices.
*
* \author Martin Albrecht <M.R.Albrecht@rhul.ac.uk>
*

# File src/pluq_mmpf.c

if(found==0) {
/* undo clearing */
-          /* alternative idea: copy out the stripe and do ordinary
-           * Gaussian elimination then use this stripe as a lookup but
-           * add stuff to the actual matrix too?*/
+          /* this is brain dead! */
+          /* don't undo it, P encodes length already applied anyway */
for(l = curr_pos; l != 0; l--)
if(mzd_read_bit(A, i, start_col + l - 1))
mzd_row_add_offset(A, i, start_row + l - 1, start_col + l);

# File src/pluq_mmpf.h

/**
* \file pluq_mmpf.h
- * \brief PLUQ factorisation using Gray codes
+ * \brief PLUQ factorization using Gray codes.
+ *
*
* \author Martin Albrecht <M.R.Albrecht@rhul.ac.uk>
+ *
+ * \example testsuite/test_lqup.c
*/

#include "permutation.h"

/**
- * Perform PLUQ factorization on A using Gray codes.
+ * \brief PLUQ matrix decomposition of A using Gray codes.
*
* \param A Matrix.
- * \param k Size of Gray code tables.
* \param P Preallocated row permutation.
* \param Q Preallocated column permutation.
+ * \param k Size of Gray code tables.
*
* \wordoffset
*/
size_t _mzd_pluq_mmpf(packedmatrix *A, permutation * P, permutation * Q, int k);

/**
- * Perform PLUQ factorization on a submatrix of up to dimension k
- * starting at (r,c).
+ * \brief PLUQ matrix decomposition of a submatrix of up to dimension
+ * k starting at (r,c).
*
* \param A Matrix.
* \param r Row Offset.

# File src/solve.h

/**
* \file solve.h
*
- * \brief system solving with matrix routines.
+ * \brief System solving with matrix routines.
*
* \author Jean-Guillaume Dumas <Jean-Guillaume.Dumas@imag.fr>
*
+ * \attention This file is currently broken.
*/
#ifndef SOLVE_H
#define SOLVE_H
#include "packedmatrix.h"

/**
- * \brief System solving with matrix.
+ * \brief Solves A X = B with A and B matrices.
*
- * Solves A X = B where A, and B are input matrices.  The solution X
- * is stored inplace on B.
+ * The solution X is stored inplace on B.
*
* \param A Input matrix.
* \param B Input matrix, being overwritten by the solution matrix X
const int inconsistency_check);

/**
- * \brief System solving with matrix.
+ * \brief Solves (P L U Q) X = B
+ *
+ * A is an input matrix supposed to store both:
+ * \li  an upper right triangular matrix U
+ * \li  a lower left unitary triangular matrix L.
*
- * Solves (P L U Q) X = B where
- * A is an input matrix supposed to store both:
- *      an upper right triangular matrix U
- *      a lower left unitary triangular matrix L
- * P and Q are permutation matrices
- * B is the input matrix to be solved.
* The solution X is stored inplace on B
+ *
* This version assumes that the matrices are at an even position on
* the RADIX grid and that their dimension is a multiple of RADIX.
*
packedmatrix *B, const int cutoff, const int inconsistency_check);

/**
- * \brief System solving with matrix.
+ * \brief  Solves (P L U Q) X = B
*
- * Solves (P L U Q) X = B where
* A is an input matrix supposed to store both:
- *      an upper right triangular matrix U
- *      a lower left unitary triangular matrix L
- * P and Q are permutation matrices
- * B is the input matrix to be solved.
- * The solution X is stored inplace on B
+ * \li an upper right triangular matrix U
+ * \li a lower left unitary triangular matrix L.
+
+ * The solution X is stored inplace on B.
+ *
* This version assumes that the matrices are at an even position on
* the RADIX grid and that their dimension is a multiple of RADIX.
*
packedmatrix *B, const int cutoff, const int inconsistency_check);

/**
- * \brief System solving with matrix.
+ * \brief Solves A X = B with A and B matrices.
*
- * Solves A X = B where A, and B are input matrices.
- * The solution X is stored inplace on B
+ * The solution X is stored inplace on B.
+ *
* This version assumes that the matrices are at an even position on
* the RADIX grid and that their dimension is a multiple of RADIX.
*
* \param inconsistency_check decide whether or not to check for
*        incosistency (faster without but output not defined if
*        system is not consistent).
- *
*/
void _mzd_solve_left(packedmatrix *A, packedmatrix *B, const int cutoff, const int inconsistency_check);

-
#endif // SOLVE_H

# File src/strassen.c

#ifdef HAVE_OPENMP
packedmatrix *_mzd_mul_mp_even(packedmatrix *C, packedmatrix *A, packedmatrix *B, int cutoff) {
/**
-   * \todo: make sure not to overwrite crap after ncols and before width*RADIX
+   * \todo make sure not to overwrite crap after ncols and before width*RADIX
*/
size_t a,b,c;
size_t anr, anc, bnr, bnc;

packedmatrix *_mzd_addmul_even(packedmatrix *C, packedmatrix *A, packedmatrix *B, int cutoff) {
/**
-   * \todo: make sure not to overwrite crap after ncols and before width*RADIX
+   * \todo make sure not to overwrite crap after ncols and before width*RADIX
*/

size_t a,b,c;

# File src/trsm.c

|B0 |B1 |
|___|___|
\endverbatim
-
-   * U00 and B0 are possibly located at uneven locations.
-   * Their column dimension is lower than 64
-   * The first column of U01, U11, B1 are aligned to words.
+     * \li U00 and B0 are possibly located at uneven locations.
+   * \li Their column dimension is lower than 64.
+   * \li The first column of U01, U11, B1 are aligned at words.
*/
packedmatrix *B0  = mzd_init_window (B,  0,  0, mb, n1);
packedmatrix *B1  = mzd_init_window (B,  0, n1, mb, nb);
else{
/**
\verbatim
-
|\
| \
|  \
|B0  |B1  |
|____|____|
\endverbatim
-
-   * L00 and B0 are possibly located at uneven locations.
-   * Their column dimension is lower than 64
-   * The first column of L10, L11, B1 are aligned to words.
+   * \li L00 and B0 are possibly located at uneven locations.
+   * \li Their column dimension is lower than 64.
+   * \li The first column of L10, L11, B1 are aligned to words.
*/
packedmatrix *B0  = mzd_init_window (B,  0,  0, mb, n1);
packedmatrix *B1  = mzd_init_window (B,  0, n1, mb, nb);
|L10 |L11\  |      |
|____|____\ |______|
\endverbatim
-    * L00 L10 B0 and B1 are possibly located at uneven locations.
-    * Their column dimension is lower than 64
-    * The first column of L01, L11, B1 are aligned to words.
+    * \li L00 L10 B0 and B1 are possibly located at uneven locations.
+    * \li Their column dimension is lower than 64.
+    * \li The first column of L01, L11, B1 are aligned to words.
*/

packedmatrix *B0  = mzd_init_window (B,  0,  0, m1, nb);
\ | |      |
\| |______|
\endverbatim
-     * U00, B0 and B1 are possibly located at uneven locations.
-     * Their column dimension is greater than 64
-     * The first column of U01, U11, B0 and B1 are aligned to words.
+     * \li U00, B0 and B1 are possibly located at uneven locations.
+     * \li Their column dimension is greater than 64
+     * \li The first column of U01, U11, B0 and B1 are aligned to words.
*/

packedmatrix *B0  = mzd_init_window (B,  0,  0, m1, nb);

# File src/trsm.h

#include "packedmatrix.h"

/**
- * \brief TRiangular System solving with Matrix.
+ * \brief Solves X U = B with X and B matrices and U upper triangular.
*
- * Solves X U = B where X and B are matrices, and U is upper
- * triangular.  X is stored inplace on B.
+ * X is stored inplace on B.
*
+ * \attention Note, that the 'right' variants of TRSM are slower than
+ * the 'left' variants.
+ *
* This is the wrapper function including bounds checks. See
- * _mzd_trsm_upper_right for implementation details.
- *
+ * _mzd_trsm_upper_right() for implementation details.
+ *
* \param U Input upper triangular matrix.
* \param B Input matrix, being overwritten by the solution matrix X
* \param cutoff Minimal dimension for Strassen recursion.
void mzd_trsm_upper_right(packedmatrix *U, packedmatrix *B, const int cutoff);

/**
- * \brief TRiangular System solving with Matrix.
+ * \brief Solves X U = B with X and B matrices and U upper triangular.
*
- * Solves X U = B where X and B are matrices, and U is upper triangular.
+ * X is stored inplace on B.
+ *
+ * \attention Note, that the 'right' variants of TRSM are slower than
+ * the 'left' variants.
*
* \param U Input upper triangular matrix.
* \param B Input matrix, being overwritten by the solution matrix X
* \param cutoff Minimal dimension for Strassen recursion.
- *
- * \internal
*/
void _mzd_trsm_upper_right(packedmatrix *U, packedmatrix *B, const int cutoff);

/**
- * \brief TRiangular System solving with Matrix.
+ * \brief Solves X L = B with X and B matrices and L lower triangular.
*
- * Solves X L = B where X and B are matrices, and L is lower triangular.
- * X is stored inplace on B
- * *
+ * X is stored inplace on B.
+ *
* This is the wrapper function including bounds checks. See
- * _mzd_trsm_upper_right for implementation details.
+ * _mzd_trsm_upper_right() for implementation details.
+ *
+ * \attention Note, that the 'right' variants of TRSM are slower than the 'left'
+ * variants.
*
* \param L Input upper triangular matrix.
* \param B Input matrix, being overwritten by the solution matrix X
void mzd_trsm_lower_right(packedmatrix *L, packedmatrix *B, const int cutoff);

/**
- * \brief TRiangular System solving with matrix.
+ * \brief Solves X L = B with X and B with matrices and L lower
+ * triangular.
+ *
+ * This version assumes that the matrices are at an even position on
+ * the RADIX grid and that their dimension is a multiple of RADIX.  X
+ * is stored inplace on B.
*
- * Solves X L = B where X and B are matrices, and L is lower triangular.
- * This version assumes that the matrices are at an even position on
- * the RADIX grid and that their dimension is a multiple of RADIX.
- * X is stored inplace on B
+ * \attention Note, that the 'right' variants of TRSM are slower than
+ * the 'left' variants.
*
* \param L Input lower triangular matrix.
* \param B Input matrix, being overwritten by the solution matrix X
void _mzd_trsm_lower_right(packedmatrix *L, packedmatrix *B, const int cutoff);

/**
- * \brief TRiangular System solving with matrix.
- *
- * Solves L X = B where X and B are matrices, and L is lower triangular.
+ * \brief Solves L X = B with X and B matrices and L lower triangular.
+ *
* X is stored inplace on B.
*
* This is the wrapper function including bounds checks. See
- * _mzd_trsm_lower_left for implementation details.
+ * _mzd_trsm_lower_left() for implementation details.
*
* \param L Input lower triangular matrix.
* \param B Input matrix, being overwritten by the solution matrix X
void mzd_trsm_lower_left(packedmatrix *L, packedmatrix *B, const int cutoff);

/**
- * \brief TRiangular System solving with Matrix.
- *
- * Solves L X = B where X and B are matrices, and L is lower
- * triangular.  X is stored inplace on B.
+ * \brief Solves L X = B with X and B matrices and L lower triangular.
+ *
+ * X is stored inplace on B.
*
* \param L Input lower triangular matrix.
* \param B Input matrix, being overwritten by the solution matrix X
void _mzd_trsm_lower_left(packedmatrix *L, packedmatrix *B, const int cutoff);

/**
- * \brief TRiangular System solving with matrix.
+ * \brief Solves U X = B with X and B matrices and U upper triangular.
*
- * Solves U X = B where X and B are matrices, and U is upper triangular.
- *  X is stored inplace on B
+ * X is stored inplace on B.
*
* This is the wrapper function including bounds checks. See
- * _mzd_trsm_upper_left for implementation details.
+ * _mzd_trsm_upper_left() for implementation details.
*
* \param U Input upper triangular matrix.
* \param B Input matrix, being overwritten by the solution matrix X
void mzd_trsm_upper_left(packedmatrix *U, packedmatrix *B, const int cutoff);

/**
- * \brief TRiangular System solving with matrix.
+ * \brief Solves U X = B with X and B matrices and U upper triangular.
*
- * Solves U X = B where X and B are matrices, and U is upper triangular.
- * X is stored inplace on B
+ * X is stored inplace on B.
*
* \param U Input upper triangular matrix.
* \param B Input matrix, being overwritten by the solution matrix X

# File testsuite/bench_elimination.c

}
wt = walltime(&clockZero);
t = cpucycles();
-  size_t r = mzd_reduce_m4ri(A, 1, 0, NULL, NULL);
+  size_t r = mzd_echelonize_m4ri(A, 1, 0, NULL, NULL);
printf("m: %5d, n: %5d, r: %5d, cpu cycles: %llu wall time: %lf\n",m, n, r, cpucycles() - t, walltime(&wt));

mzd_free(A);

# File testsuite/test_elimination.c

E = mzd_copy(NULL, A);
F = mzd_copy(NULL, A);

-  mzd_reduce_m4ri(A, 1, 0, NULL, NULL);
+  mzd_echelonize_m4ri(A, 1, 0, NULL, NULL);

-  mzd_reduce_m4ri(B, 1, 8, NULL, NULL);
+  mzd_echelonize_m4ri(B, 1, 8, NULL, NULL);

-  mzd_reduce_m4ri(C, 0, 0, NULL, NULL);
-  mzd_top_reduce_m4ri(C, 0, NULL, NULL);
+  mzd_echelonize_m4ri(C, 0, 0, NULL, NULL);
+  mzd_top_echelonize_m4ri(C, 0, NULL, NULL);

-  mzd_reduce_m4ri(D, 0, 4, NULL, NULL);
-  mzd_top_reduce_m4ri(D, 4, NULL, NULL);
+  mzd_echelonize_m4ri(D, 0, 4, NULL, NULL);
+  mzd_top_echelonize_m4ri(D, 4, NULL, NULL);

-  mzd_reduce_naive(E, 1);
+  mzd_echelonize_naive(E, 1);

-  mzd_reduce_naive(F, 0);
-  mzd_top_reduce_m4ri(F, 0, NULL, NULL);
+  mzd_echelonize_naive(F, 0);
+  mzd_top_echelonize_m4ri(F, 0, NULL, NULL);

if(mzd_equal(A, B) != TRUE) {
printf("A != B ");