[mathicgb] 243/393: Fixed issue where memory consumption would suddenly explode due to converting a submatrix of a very sparse and large matrix into dense format. Very sparse matrices now have their own code-path that avoids this.

[Doug Torrance]

This is an automated email from the git hooks/post-receive script.

dtorrance-guest pushed a commit to branch upstream
in repository mathicgb.

commit c35380f26106377ff17571db4de96b958faf272e
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date:   Sat Apr 13 13:19:27 2013 -0400

    Fixed issue where memory consumption would suddenly explode due to converting a submatrix of a very sparse and large matrix into dense format. Very sparse matrices now have their own code-path that avoids this.
---
 src/mathicgb/F4MatrixReducer.cpp | 89 +++++++++++++++++++++++++++++++++++++++-
 src/mathicgb/F4MatrixReducer.hpp |  2 +-
 2 files changed, 88 insertions(+), 3 deletions(-)

diff --git a/src/mathicgb/F4MatrixReducer.cpp b/src/mathicgb/F4MatrixReducer.cpp
old mode 100644
new mode 100755
index f11a54d..c6e63d5
--- a/src/mathicgb/F4MatrixReducer.cpp
+++ b/src/mathicgb/F4MatrixReducer.cpp
@@ -306,6 +306,84 @@ namespace {
     return std::move(reduced);
   }
 
+  SparseMatrix reduceToEchelonFormSparse(
+    const SparseMatrix& toReduce,
+    const SparseMatrix::Scalar modulus
+  ) {
+    const auto colCount = toReduce.computeColCount();
+
+    const auto noRow = static_cast<SparseMatrix::RowIndex>(-1);
+
+    // pivotRowOfCol[i] is the pivot in column i or noRow
+    // if we have not identified such a pivot so far.
+    std::vector<SparseMatrix::RowIndex> pivotRowOfCol(colCount, noRow);
+
+    DenseRow rowToReduce(colCount);
+
+    // ** Reduce to row echelon form -- every row is a pivot row.
+    SparseMatrix pivots(colCount);
+    for (SparseMatrix::RowIndex row = 0; row < toReduce.rowCount(); ++row) {
+      if (toReduce.emptyRow(row))
+        continue;
+      rowToReduce.clear(colCount);
+      rowToReduce.addRow(toReduce, row);
+
+      SparseMatrix::ColIndex leadingCol = 0;
+      while (true) { // reduce row by previous pivots
+        for (; leadingCol < colCount; ++leadingCol) {
+          auto& entry = rowToReduce[leadingCol];
+          if (entry != 0) {
+            entry %= modulus;
+            if (entry != 0)
+              break;
+          }
+        }
+        if (leadingCol == colCount)
+          break; // The row has been reduced to zero.
+        const auto pivotRow = pivotRowOfCol[leadingCol];
+        if (pivotRow == noRow) { // If the row is a new pivot.
+          rowToReduce.makeUnitary(modulus, leadingCol);
+          pivotRowOfCol[leadingCol] = pivots.rowCount();
+          rowToReduce.appendTo(pivots);
+          break;
+        }
+        rowToReduce.rowReduceByUnitary(pivotRow, pivots, modulus);
+      }
+    }
+
+    // ** Reduce from row echelon form to reduced row echelon form
+
+    SparseMatrix reduced(colCount);
+    auto pivotCol = colCount;
+    // Reduce pivot rows in descending order of leading column. The reduced
+    // pivots go into reduced and we update pivotRowOfCol to refer to the
+    // row indices in reduced as we go along.
+    while (pivotCol != 0) {
+      --pivotCol;
+      const auto row = pivotRowOfCol[pivotCol];
+      if (row == noRow)
+        continue;
+      rowToReduce.clear(colCount);
+      rowToReduce.addRow(pivots, row);
+      MATHICGB_ASSERT(rowToReduce[pivotCol] == 1); // unitary
+      for (auto col = pivotCol + 1; col != colCount; ++col) {
+        auto& entry = rowToReduce[col];
+        if (entry == 0)
+          continue;
+        entry %= modulus;
+        if (entry == 0)
+          continue;
+        const auto pivotRow = pivotRowOfCol[col];
+        if (pivotRow != noRow)
+          rowToReduce.rowReduceByUnitary(pivotRow, reduced, modulus);
+        MATHICGB_ASSERT(entry < modulus);
+      }
+      pivotRowOfCol[pivotCol] = reduced.rowCount();
+      rowToReduce.appendTo(reduced);
+    }
+    return std::move(reduced);
+  }
+
   SparseMatrix reduceToEchelonForm(
     const SparseMatrix& toReduce,
     const SparseMatrix::Scalar modulus
@@ -698,8 +776,15 @@ SparseMatrix F4MatrixReducer::reducedRowEchelonForm(
       return reduceToEchelonFormShrawanDelayedModulus(matrix, mModulus);
     else    
       return reduceToEchelonFormShrawan(matrix, mModulus);
-  } else
-    return reduceToEchelonForm(matrix, mModulus);
+  } else {
+    // todo: actually do some work to determine a good way to determine
+    // when to use the sparse method, or alternatively make some some
+    // sort of hybrid.
+    if ( matrix.computeDensity() < 0.02)
+      return reduceToEchelonFormSparse(matrix, mModulus);
+    else
+      return reduceToEchelonForm(matrix, mModulus);
+  }
 }
 
 SparseMatrix F4MatrixReducer::reducedRowEchelonFormBottomRight(
diff --git a/src/mathicgb/F4MatrixReducer.hpp b/src/mathicgb/F4MatrixReducer.hpp
index 4fc76a5..92b2f33 100755
--- a/src/mathicgb/F4MatrixReducer.hpp
+++ b/src/mathicgb/F4MatrixReducer.hpp
@@ -27,7 +27,7 @@ public:
   /// Returns the reduced row echelon form of matrix.
   SparseMatrix reducedRowEchelonForm(const SparseMatrix& matrix);
 
-  /// Returns the lower right submatrix if the reduced row echelon
+  /// Returns the lower right submatrix of the reduced row echelon
   /// form of matrix. The lower left part is not returned because it is
   /// always zero after row reduction.
   SparseMatrix reducedRowEchelonFormBottomRight(const QuadMatrix& matrix);

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/mathicgb.git