+#ifndef matrixreduction_h_aPeinoXeiFeethuz
+#define matrixreduction_h_aPeinoXeiFeethuz
+
+#include <iostream>
+
+#include "booleanmatrix.h"
+
+
+namespace libstick {
+
+/** Takes the boundary matrix b and transforms it inplace into a reduced matrix
+ * b. The matrix v is required to be a unit matrix having the same dimensions
+ * as b. It is maintained in order to leave the product b*inverse(v) as an
+ * invariant. Hence, the resulting b is equal to the product of the boundary
+ * matrix times v. */
+template<class IT, class D>
+void reduceBoundaryMatrix(BooleanColRowMatrix<IT> &b, BooleanColMatrix_base<IT, D> &v) {
+ assert(b.size() == v.width());
+
+ // Make every column reduced, i.e., it is a zero-column or contains only one 1.
+ for (unsigned c=0; c < b.size(); ++c) {
+ const typename BooleanColRowMatrix<IT>::colbase::column_type &col = b.getColumn(c);
+ if (col.size() == 0)
+ continue;
+
+ // The column is not a zero-column, take the lowest 1, and reduce the other columns at the row
+ IT r = col.back();
+ assert(b.get(r, c));
+
+ // Get a copy of the r-th row
+ typename BooleanColRowMatrix<IT>::rowbase::row_type row(b.getRow(r).size());
+ copy(b.getRow(r).begin(), b.getRow(r).end(), row.begin());
+ assert(row.size() >= 1);
+
+ // Get rid of 1s at that row right of column c.
+ typename BooleanColRowMatrix<IT>::row_type::const_iterator it = lower_bound(row.begin(), row.end(), c);
+ for(; it != row.end(); ++it) {
+ if (*it <= c)
+ continue;
+
+ assert(b.get(r, *it));
+ b.addColumn(*it, col);
+ v.addColumn(*it, v.getColumn(c));
+ assert(!b.get(r, *it));
+ }
+ }
+}
+
+}
+
+#endif