linalg: fix some bugs in row echelon form calculation

src/linalg.rs

@@ -7,6 +7,7 @@ pub enum Error {
 
 pub fn gauss_reduction<A: LinalgScalar + PartialEq>(array: &mut Array2<A>) -> Result<(), Error> {
     let (n, m) = array.dim();
+    assert!(n <= m);
 
     // Reduce to upper triangular with ones on diagonal
     for j in 0..n {
@@ -67,9 +68,9 @@ pub fn gauss_reduction<A: LinalgScalar + PartialEq>(array: &mut Array2<A>) -> Re
 pub fn row_echelon_form<A: LinalgScalar + PartialEq>(array: &mut Array2<A>) {
     let (n, m) = array.dim();
 
-    // Reduce to upper triangular with ones on diagonal
+    let mut j = 0;
     let mut k = 0;
-    for j in 0..n {
+    while j < m && k < n {
         // Find non-zero element in current column, at or below row k
         let Some(s) = array
             .slice(s![k.., j])
@@ -79,6 +80,7 @@ pub fn row_echelon_form<A: LinalgScalar + PartialEq>(array: &mut Array2<A>) {
         else {
             // All the elements at or below row k are zero. Done with this
             // column.
+            j += 1;
             continue;
         };
 
@@ -102,6 +104,7 @@ pub fn row_echelon_form<A: LinalgScalar + PartialEq>(array: &mut Array2<A>) {
             }
         }
 
+        j += 1;
         k += 1;
     }
 }

src/systematic.rs

@@ -39,6 +39,12 @@ pub fn parity_to_systematic(h: &SparseMatrix) -> Result<SparseMatrix, Error> {
         a[[j, k]] = GF2::one();
     }
     linalg::row_echelon_form(&mut a);
+    // Check that the matrix has full rank by checking that there is a non-zero
+    // element in the last row (we start looking by the end, since chances are
+    // higher to find a non-zero element there).
+    if !(0..m).rev().any(|j| a[[n - 1, j]] != GF2::zero()) {
+        return Err(Error::NotFullRank);
+    }
     // write point for columns that do not "go down" in the row echelon form
     let mut k = 0;
     let mut h_new = SparseMatrix::new(n, m);
@@ -65,10 +71,7 @@ pub fn parity_to_systematic(h: &SparseMatrix) -> Result<SparseMatrix, Error> {
                 break;
             }
         }
-        if !found {
-            // No column "going down" found. The matrix doesn't have full rank.
-            return Err(Error::NotFullRank);
-        }
+        assert!(found);
     }
     // Insert remaining columns at the write point
     for j in j0..m {