[mathicgb] 107/393: Now using a single concurrent hash map for building the matrix. It makes no difference single-core, but now using more threads than cores is actually a win - I'm seeing 12% speed-up of 3 threads on dual core compared to 2 threads where before there was no advantage. I suspect that this is caused by the third thread no longer having to create its own hash table and also now the third thread can use the idle core if the second thread is waiting while the first thread locked the hash table for writes. Hyperthreading may have an impact too, though I am skeptical of that.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:42 UTC 2015
This is an automated email from the git hooks/post-receive script.
dtorrance-guest pushed a commit to branch upstream
in repository mathicgb.
commit d25ac48ddcbdc06efd2fad70baa4f644df275655
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Wed Nov 7 19:44:01 2012 +0100
Now using a single concurrent hash map for building the matrix. It makes no difference single-core, but now using more threads than cores is actually a win - I'm seeing 12% speed-up of 3 threads on dual core compared to 2 threads where before there was no advantage. I suspect that this is caused by the third thread no longer having to create its own hash table and also now the third thread can use the idle core if the second thread is waiting while the first thread locked the hash tab [...]
---
src/mathicgb/F4MatrixBuilder.cpp | 371 +++++++++++++++---------------------
src/mathicgb/F4MatrixBuilder.hpp | 48 ++---
src/mathicgb/FixedSizeMonomialMap.h | 133 +++++++++++--
src/mathicgb/MonomialMap.hpp | 29 ++-
src/mathicgb/QuadMatrix.cpp | 99 +++++++++-
src/mathicgb/QuadMatrix.hpp | 5 +
src/mathicgb/QuadMatrixBuilder.cpp | 264 ++++---------------------
src/mathicgb/QuadMatrixBuilder.hpp | 127 +++---------
src/test/QuadMatrixBuilder.cpp | 118 ++++++------
src/test/SparseMatrix.cpp | 2 +-
10 files changed, 531 insertions(+), 665 deletions(-)
diff --git a/src/mathicgb/F4MatrixBuilder.cpp b/src/mathicgb/F4MatrixBuilder.cpp
index 905281a..735eba5 100755
--- a/src/mathicgb/F4MatrixBuilder.cpp
+++ b/src/mathicgb/F4MatrixBuilder.cpp
@@ -5,85 +5,43 @@
#include <omp.h>
#endif
-MATHICGB_NO_INLINE QuadMatrixBuilder::LeftRightColIndex
- F4MatrixBuilder::findOrCreateColumn
-(
+MATHICGB_NO_INLINE
+std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+F4MatrixBuilder::findOrCreateColumn(
const const_monomial monoA,
- const const_monomial monoB,
- QuadMatrixBuilder& builder
+ const const_monomial monoB
) {
MATHICGB_ASSERT(!monoA.isNull());
MATHICGB_ASSERT(!monoB.isNull());
- const auto col = builder.findColumnProduct(monoA, monoB);
- if (col.valid())
- return col;
- return createColumn(builder, monoA, monoB);
+ const auto col = ColReader(mMap).findProduct(monoA, monoB);
+ if (col.first != 0)
+ return std::make_pair(*col.first, col.second);
+ return createColumn(monoA, monoB);
}
-MATHICGB_INLINE QuadMatrixBuilder::LeftRightColIndex
- F4MatrixBuilder::findOrCreateColumn
-(
+MATHICGB_INLINE
+std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+F4MatrixBuilder::findOrCreateColumn(
const const_monomial monoA,
const const_monomial monoB,
- const ColReader& colMap,
- QuadMatrixBuilder& builder
+ const ColReader& colMap
) {
MATHICGB_ASSERT(!monoA.isNull());
MATHICGB_ASSERT(!monoB.isNull());
const auto col = colMap.findProduct(monoA, monoB);
- if (col == 0)
- return findOrCreateColumn(monoA, monoB, builder);
- return *col;
-}
-
-MATHICGB_NO_INLINE const std::pair<
- QuadMatrixBuilder::LeftRightColIndex,
- QuadMatrixBuilder::LeftRightColIndex
-> F4MatrixBuilder::findOrCreateTwoColumns
-(
- const const_monomial monoA1,
- const const_monomial monoA2,
- const const_monomial monoB,
- QuadMatrixBuilder& builder
-) {
- MATHICGB_ASSERT(!monoA1.isNull());
- MATHICGB_ASSERT(!monoA2.isNull());
- MATHICGB_ASSERT(!monoB.isNull());
- MATHICGB_ASSERT(!ring().monomialEQ(monoA1, monoA2));
- auto colPair = builder.findTwoColumnsProduct(monoA1, monoA2, monoB);
- if (!colPair.first.valid()) {
- colPair.first = createColumn(builder, monoA1, monoB);
- // Syncing builder to mBuilder could have created a col for monoA2*monoB.
- if (&mBuilder != &builder && !colPair.second.valid())
- colPair.second = builder.findColumnProduct(monoA2, monoB);
- }
- if (!colPair.second.valid())
- colPair.second = createColumn(builder, monoA2, monoB);
- MATHICGB_ASSERT(colPair == std::make_pair(
- findOrCreateColumn(monoA1, monoB, builder),
- findOrCreateColumn(monoA2, monoB, builder)));
- return colPair;
+ if (col.first == 0)
+ return findOrCreateColumn(monoA, monoB);
+ return std::make_pair(*col.first, col.second);
}
-MATHICGB_INLINE const std::pair<
- QuadMatrixBuilder::LeftRightColIndex,
- QuadMatrixBuilder::LeftRightColIndex
-> F4MatrixBuilder::findOrCreateTwoColumns
-(
+MATHICGB_NO_INLINE
+void F4MatrixBuilder::createTwoColumns(
const const_monomial monoA1,
const const_monomial monoA2,
- const const_monomial monoB,
- const ColReader& colMap,
- QuadMatrixBuilder& builder
+ const const_monomial monoB
) {
- MATHICGB_ASSERT(!monoA1.isNull());
- MATHICGB_ASSERT(!monoA2.isNull());
- MATHICGB_ASSERT(!monoB.isNull());
- MATHICGB_ASSERT(!ring().monomialEQ(monoA1, monoA2));
- auto colPair = colMap.findTwoProducts(monoA1, monoA2, monoB);
- if (colPair.first == 0 || colPair.second == 0)
- return findOrCreateTwoColumns(monoA1, monoA2, monoB, builder);
- return std::make_pair(*colPair.first, *colPair.second);
+ createColumn(monoA1, monoB);
+ createColumn(monoA2, monoB);
}
F4MatrixBuilder::F4MatrixBuilder(
@@ -92,9 +50,14 @@ F4MatrixBuilder::F4MatrixBuilder(
const size_t memoryQuantum
):
mThreadCount(threadCount),
+ mTmp(basis.ring().allocMonomial()),
mBasis(basis),
- mBuilder(basis.ring(), memoryQuantum),
- mTmp(basis.ring().allocMonomial())
+ mMap(basis.ring()),
+ mMonomialsLeft(),
+ mMonomialsRight(),
+ mBuilder(basis.ring(), mMap, mMonomialsLeft, mMonomialsRight, memoryQuantum),
+ mLeftColCount(0),
+ mRightColCount(0)
{
MATHICGB_ASSERT(threadCount >= 1);
@@ -172,21 +135,17 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
// Process pending rows until we are done. Note that the methods
// we are calling here can add more items to mTodo.
+ QuadMatrixBuilder mainBuilder(
+ ring(), mMap, mMonomialsLeft, mMonomialsRight, mBuilder.memoryQuantum()
+ );
+
#ifdef _OPENMP
- struct ThreadData {
- QuadMatrixBuilder* builder;
- monomial tmp;
- };
MATHICGB_ASSERT(mThreadCount >= 1);
- std::vector<ThreadData> threadData(mThreadCount);
- if (mThreadCount == 1) {
- threadData[0].builder = &mBuilder;
- threadData[0].tmp = mTmp;
- } else {
- for (size_t i = 0; i < threadData.size(); ++i) {
- threadData[i].builder = new QuadMatrixBuilder(ring());
- threadData[i].tmp = new exponent[ring().maxMonomialByteSize()];
- }
+ std::vector<QuadMatrixBuilder*> threadData(mThreadCount);
+ for (size_t i = 0; i < threadData.size(); ++i) {
+ threadData[i] = i == 0 ? &mainBuilder : new QuadMatrixBuilder(
+ ring(), mMap, mMonomialsLeft, mMonomialsRight, mBuilder.memoryQuantum()
+ );
}
#endif
@@ -209,10 +168,9 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
const size_t taskIndex = taskOMP;
#ifdef _OPENMP
MATHICGB_ASSERT(omp_get_thread_num() < threadData.size());
- const ThreadData& td = threadData[omp_get_thread_num()];
- QuadMatrixBuilder& builder = *td.builder;
+ QuadMatrixBuilder& builder = *threadData[omp_get_thread_num()];
#else
- QuadMatrixBuilder& builder = mBuilder;
+ QuadMatrixBuilder& builder = mainBuilder;
#endif
const RowTask task = currentTasks[taskIndex];
@@ -228,137 +186,91 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
#ifdef _OPENMP
#pragma omp flush
- MATHICGB_ASSERT(mThreadCount > 1 || threadData.back().builder == &mBuilder);
- if (mThreadCount > 1) {
- for (auto it = threadData.begin(); it != threadData.end(); ++it) {
- MATHICGB_ASSERT(&mBuilder != it->builder);
- QuadMatrix qm(it->builder->buildMatrixAndClear());
- mBuilder.takeRowsFrom(std::move(qm));
- delete it->builder;
+ for (auto it = threadData.begin(); it != threadData.end(); ++it) {
+ if (&mainBuilder != *it) {
+ mainBuilder.takeRowsFrom((*it)->buildMatrixAndClear());
+ delete *it;
}
}
threadData.clear();
#endif
- mBuilder.sortColumnsLeftRightParallel(mBasis.order(), mThreadCount);
- matrix = mBuilder.buildMatrixAndClear();
+ matrix = mainBuilder.buildMatrixAndClear();
+ {
+ ColReader reader(mMap);
+ matrix.leftColumnMonomials.clear();
+ matrix.rightColumnMonomials.clear();
+ const auto end = reader.end();
+ for (auto it = reader.begin(); it != end; ++it) {
+ const auto p = *it;
+ monomial copy = ring().allocMonomial();
+ ring().monomialCopy(p.second, copy);
+ auto& monos = p.first.left() ?
+ matrix.leftColumnMonomials : matrix.rightColumnMonomials;
+ const auto index = p.first.index();
+ if (monos.size() <= index)
+ monos.resize(index + 1);
+ MATHICGB_ASSERT(monos[index].isNull());
+ monos[index] = copy;
+ }
+ }
+#ifdef MATHICGB_DEBUG
+ for (size_t side = 0; side < 2; ++side) {
+ auto& monos = side == 0 ?
+ matrix.leftColumnMonomials : matrix.rightColumnMonomials;
+ for (auto it = monos.begin(); it != monos.end(); ++it) {
+ MATHICGB_ASSERT(!it->isNull());
+ }
+ }
+#endif
+ matrix.sortColumnsLeftRightParallel(mThreadCount);
+ mMap.clearNonConcurrent();
}
-F4MatrixBuilder::LeftRightColIndex
+std::pair<F4MatrixBuilder::LeftRightColIndex, ConstMonomial>
F4MatrixBuilder::createColumn(
- QuadMatrixBuilder& builder,
const const_monomial monoA,
const const_monomial monoB
) {
MATHICGB_ASSERT(!monoA.isNull());
MATHICGB_ASSERT(!monoB.isNull());
- // Must be declared up here since we cannot return out of an omp critical
- // section.
- QuadMatrixBuilder::LeftRightColIndex newCol;
-#pragma omp critical
+ std::lock_guard<std::mutex> lock(mCreateColumnLock);
+ // see if the column exists now after we have synchronized
{
- ring().monomialMult(monoA, monoB, mTmp);
-
- MATHICGB_ASSERT(ring().hashValid(mTmp));
- MATHICGB_ASSERT(!builder.findColumn(mTmp).valid());
-#ifndef _OPENMP
- MATHICGB_ASSERT(&builder == &mBuilder);
-#endif
-
-#ifdef _OPENMP
- // builder does not have a column for mTmp but mBuilder might
- if (&mBuilder != &builder)
- newCol = mBuilder.findColumn(mTmp);
- if (!newCol.valid())
-#endif
- {
- // we need to add a new column to mBuilder
-
- // look for a reducer of mTmp
- size_t reducerIndex = mBasis.divisor(mTmp);
- if (reducerIndex == static_cast<size_t>(-1)) {
- newCol = mBuilder.createColumnRight(mTmp);
- } else {
- newCol = mBuilder.createColumnLeft(mTmp);
-
- // schedule the reducer to be added as a row
- RowTask task = {};
- task.addToTop = true;
- task.poly = &mBasis.poly(reducerIndex);
- task.multiply = ring().allocMonomial();
- ring().monomialDivideToNegative
- (mTmp, task.poly->getLeadMonomial(), task.multiply);
- MATHICGB_ASSERT(ring().hashValid(task.multiply));
- mTodo.push_back(task);
- }
- }
- MATHICGB_ASSERT(newCol.valid());
- MATHICGB_ASSERT(ring().monomialEQ(mBuilder.monomialOfCol(newCol), mTmp));
- MATHICGB_ASSERT(newCol == mBuilder.findColumn(mTmp));
-
-#ifdef _OPENMP
- // If we are running in parallel we need to catch builder up to date with
- // with respect to the columns in mBuilder. We have just ensured that
- // mBuilder has a column for mTmp so now builder will also get that column.
-
- if (&builder != &mBuilder) {
- // add missing left columns
- if (builder.leftColCount() != mBuilder.leftColCount()) {
- const auto colCount = mBuilder.leftColCount();
- MATHICGB_ASSERT(builder.leftColCount() < colCount);
- for (auto col = builder.leftColCount(); col < colCount; ++col) {
- const const_monomial monoToCopy = mBuilder.monomialOfLeftCol(col);
- MATHICGB_ASSERT(!builder.findColumn(monoToCopy).valid());
- builder.createColumnLeft(monoToCopy);
- MATHICGB_ASSERT(builder.findColumn(monoToCopy) ==
- mBuilder.findColumn(monoToCopy));
- MATHICGB_ASSERT(ring().monomialEQ
- (mBuilder.monomialOfLeftCol(col), builder.monomialOfLeftCol(col)));
- }
- }
-#ifdef MATHICGB_DEBUG
- MATHICGB_ASSERT(builder.leftColCount() == mBuilder.leftColCount());
- for (SparseMatrix::ColIndex col = 0;
- col < builder.leftColCount(); ++col) {
- MATHICGB_ASSERT(ring().monomialEQ
- (mBuilder.monomialOfLeftCol(col), builder.monomialOfLeftCol(col)));
- }
-#endif
-
- // add missing right columns
- if (builder.rightColCount() != mBuilder.rightColCount()) {
- const auto colCount = mBuilder.rightColCount();
- MATHICGB_ASSERT(builder.rightColCount() < colCount);
- for (auto col = builder.rightColCount(); col < colCount; ++col) {
- const const_monomial monoToCopy = mBuilder.monomialOfRightCol(col);
- MATHICGB_ASSERT(!builder.findColumn(monoToCopy).valid());
- builder.createColumnRight(monoToCopy);
- MATHICGB_ASSERT(builder.findColumn(monoToCopy) ==
- mBuilder.findColumn(monoToCopy));
- MATHICGB_ASSERT(ring().monomialEQ
- (mBuilder.monomialOfRightCol(col), builder.monomialOfRightCol(col)));
- }
- }
-#ifdef MATHICGB_DEBUG
- MATHICGB_ASSERT(builder.rightColCount() == mBuilder.rightColCount());
- for (SparseMatrix::ColIndex col = 0;
- col < builder.rightColCount(); ++col) {
- MATHICGB_ASSERT(ring().monomialEQ
- (mBuilder.monomialOfRightCol(col), builder.monomialOfRightCol(col)));
- }
-#endif
-
- MATHICGB_ASSERT(ring().monomialEQ(mBuilder.monomialOfCol(mBuilder.findColumn(mTmp)), mTmp));
- MATHICGB_ASSERT(newCol == mBuilder.findColumn(mTmp));
+ const auto found(ColReader(mMap).findProduct(monoA, monoB));
+ if (found.first != 0)
+ return std::make_pair(*found.first, found.second);
+ }
- MATHICGB_ASSERT(builder.findColumn(mTmp).valid());
- MATHICGB_ASSERT(ring().monomialEQ(builder.monomialOfCol(builder.findColumn(mTmp)), mTmp));
- MATHICGB_ASSERT(newCol == builder.findColumn(mTmp));
- }
-#endif
+ // The column really does not exist, so we need to create it
+ ring().monomialMult(monoA, monoB, mTmp);
+ MATHICGB_ASSERT(ring().hashValid(mTmp));
+
+ // look for a reducer of mTmp
+ const size_t reducerIndex = mBasis.divisor(mTmp);
+ const bool insertLeft = (reducerIndex != static_cast<size_t>(-1));
+ if (insertLeft) {
+ RowTask task = {}; // schedule the new reducer to be added as a row
+ task.addToTop = true;
+ task.poly = &mBasis.poly(reducerIndex);
+ task.multiply = ring().allocMonomial();
+ ring().monomialDivideToNegative
+ (mTmp, task.poly->getLeadMonomial(), task.multiply);
+ MATHICGB_ASSERT(ring().hashValid(task.multiply));
+ mTodo.push_back(task);
}
- return newCol;
+
+ // Create the new left or right column
+ const auto colCount = insertLeft ? mLeftColCount : mRightColCount;
+ if (colCount == std::numeric_limits<ColIndex>::max())
+ throw std::overflow_error("Too many columns in QuadMatrix");
+ const auto inserted = mMap.insert
+ (std::make_pair(mTmp, LeftRightColIndex(colCount, insertLeft)));
+ insertLeft ? ++mLeftColCount : ++mRightColCount;
+ MATHICGB_ASSERT(inserted.second);
+ MATHICGB_ASSERT(inserted.first.first != 0);
+ return std::make_pair(*inserted.first.first, inserted.first.second);
}
void F4MatrixBuilder::appendRowBottom(
@@ -371,16 +283,26 @@ void F4MatrixBuilder::appendRowBottom(
// todo: eliminate the code-duplication between here and appendRowTop.
MATHICGB_ASSERT(!multiple.isNull());
MATHICGB_ASSERT(&builder != 0);
- const ColReader colMap(builder.columnToIndexMap());
- for (auto it = begin; it != end; ++it) {
- const auto col = findOrCreateColumn(it.getMonomial(), multiple, colMap, builder);
+ auto it = begin;
+updateReader:
+ // Use an on-stack const reader to make it as obvious as possible to the
+ // optimizer's alias analysis that the pointer inside the reader never
+ // changes inside the loop.
+ const ColReader reader(mMap);
+ for (; it != end; ++it) {
+ const auto col = reader.findProduct(it.getMonomial(), multiple);
+ if (col.first == 0) {
+ createColumn(it.getMonomial(), multiple);
+ goto updateReader;
+ }
+
const auto origScalar = it.getCoefficient();
MATHICGB_ASSERT(origScalar != 0);
- const auto possiblyNegated =
+ const auto maybeNegated =
negate ? ring().coefficientNegateNonZero(origScalar) : origScalar;
- MATHICGB_ASSERT(possiblyNegated < std::numeric_limits<Scalar>::max());
- builder.appendEntryBottom(col, static_cast<Scalar>(possiblyNegated));
+ MATHICGB_ASSERT(maybeNegated < std::numeric_limits<Scalar>::max());
+ builder.appendEntryBottom(*col.first, static_cast<Scalar>(maybeNegated));
}
builder.rowDoneBottomLeftAndRight();
}
@@ -394,35 +316,41 @@ void F4MatrixBuilder::appendRowTop(
MATHICGB_ASSERT(&poly != 0);
MATHICGB_ASSERT(&builder != 0);
- ColReader colMap(builder.columnToIndexMap());
-
auto it = poly.begin();
const auto end = poly.end();
if ((std::distance(it, end) % 2) == 1) {
- const auto col = findOrCreateColumn(it.getMonomial(), multiple, colMap, builder);
+ ColReader reader(mMap);
+ const auto col = findOrCreateColumn(it.getMonomial(), multiple, reader);
MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
MATHICGB_ASSERT(it.getCoefficient());
- builder.appendEntryTop(col, static_cast<Scalar>(it.getCoefficient()));
+ builder.appendEntryTop(col.first, static_cast<Scalar>(it.getCoefficient()));
++it;
}
+updateReader:
+ ColReader colMap(mMap);
MATHICGB_ASSERT((std::distance(it, end) % 2) == 0);
while (it != end) {
MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
MATHICGB_ASSERT(it.getCoefficient() != 0);
const auto scalar1 = static_cast<Scalar>(it.getCoefficient());
const const_monomial mono1 = it.getMonomial();
- ++it;
- MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
- MATHICGB_ASSERT(it.getCoefficient() != 0);
- const auto scalar2 = static_cast<Scalar>(it.getCoefficient());
- const const_monomial mono2 = it.getMonomial();
- ++it;
+ auto it2 = it;
+ ++it2;
+ MATHICGB_ASSERT(it2.getCoefficient() < std::numeric_limits<Scalar>::max());
+ MATHICGB_ASSERT(it2.getCoefficient() != 0);
+ const auto scalar2 = static_cast<Scalar>(it2.getCoefficient());
+ const const_monomial mono2 = it2.getMonomial();
+
+ const auto colPair = colMap.findTwoProducts(mono1, mono2, multiple);
+ if (colPair.first == 0 || colPair.second == 0) {
+ createTwoColumns(mono1, mono2, multiple);
+ goto updateReader;
+ }
- const auto colPair =
- findOrCreateTwoColumns(mono1, mono2, multiple, colMap, builder);
- builder.appendEntryTop(colPair.first, scalar1);
- builder.appendEntryTop(colPair.second, scalar2);
+ builder.appendEntryTop(*colPair.first, scalar1);
+ builder.appendEntryTop(*colPair.second, scalar2);
+ it = ++it2;
}
builder.rowDoneTopLeftAndRight();
}
@@ -454,7 +382,7 @@ void F4MatrixBuilder::appendRowBottom(
++itA;
++itB;
- const ColReader colMap(builder.columnToIndexMap());
+ const ColReader colMap(mMap);
const const_monomial mulA = task.multiply;
const const_monomial mulB = task.sPairMultiply;
@@ -473,18 +401,19 @@ void F4MatrixBuilder::appendRowBottom(
coefficient coeff = 0;
LeftRightColIndex col;
- const auto colA = findOrCreateColumn(itA.getMonomial(), mulA, colMap, builder);
- const auto colB = findOrCreateColumn(itB.getMonomial(), mulB, colMap, builder);
- const auto cmp = ring().monomialCompare
- (builder.monomialOfCol(colA), builder.monomialOfCol(colB));
+ const auto colA = findOrCreateColumn(itA.getMonomial(), mulA, colMap);
+ const auto colB = findOrCreateColumn(itB.getMonomial(), mulB, colMap);
+ const auto cmp = ring().monomialCompare(colA.second, colB.second);
+ //const auto cmp = ring().monomialCompare
+ // (builder.monomialOfCol(colA), builder.monomialOfCol(colB));
if (cmp != LT) {
coeff = itA.getCoefficient();
- col = colA;
+ col = colA.first;
++itA;
}
if (cmp != GT) {
coeff = ring().coefficientSubtract(coeff, itB.getCoefficient());
- col = colB;
+ col = colB.first;
++itB;
}
MATHICGB_ASSERT(coeff < std::numeric_limits<Scalar>::max());
diff --git a/src/mathicgb/F4MatrixBuilder.hpp b/src/mathicgb/F4MatrixBuilder.hpp
index cbc197c..ab6d943 100755
--- a/src/mathicgb/F4MatrixBuilder.hpp
+++ b/src/mathicgb/F4MatrixBuilder.hpp
@@ -20,6 +20,8 @@ private:
typedef QuadMatrixBuilder::ColIndex ColIndex;
typedef QuadMatrixBuilder::LeftRightColIndex LeftRightColIndex;
typedef QuadMatrixBuilder::Scalar Scalar;
+ typedef QuadMatrixBuilder::Map Map;
+ typedef QuadMatrixBuilder::MonomialsType Monomials;
public:
/// memoryQuantum is how much to increase the memory size by each time the
@@ -66,13 +68,14 @@ public:
const PolyRing& ring() const {return mBuilder.ring();}
private:
- typedef const QuadMatrixBuilder::ColReader ColReader;
+ typedef const MonomialMap<LeftRightColIndex>::Reader ColReader;
/** Creates a column with monomial label x and schedules a new row to
reduce that column if possible. Here x is monoA if monoB is
null and otherwise x is the product of monoA and monoB. */
- MATHICGB_NO_INLINE LeftRightColIndex createColumn
- (QuadMatrixBuilder& builder, const_monomial monoA, const_monomial monoB);
+ MATHICGB_NO_INLINE
+ std::pair<LeftRightColIndex, ConstMonomial>
+ createColumn(const_monomial monoA, const_monomial monoB);
/// Represents the task of adding a row to the matrix. If sPairPoly is null
/// then the row to add is multiply * poly. Otherwise, the row to add is
@@ -99,35 +102,36 @@ private:
QuadMatrixBuilder& builder
);
- MATHICGB_NO_INLINE LeftRightColIndex findOrCreateColumn
- (const_monomial monoA, const_monomial monoB, QuadMatrixBuilder& builder);
- MATHICGB_INLINE LeftRightColIndex findOrCreateColumn(
+ MATHICGB_NO_INLINE
+ std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+ findOrCreateColumn(const_monomial monoA, const_monomial monoB);
+
+ MATHICGB_INLINE
+ std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+ findOrCreateColumn(
const_monomial monoA,
const_monomial monoB,
- const ColReader& colMap,
- QuadMatrixBuilder& builder);
+ const ColReader& colMap
+ );
- MATHICGB_NO_INLINE const std::pair<LeftRightColIndex, LeftRightColIndex>
- findOrCreateTwoColumns(
+ MATHICGB_NO_INLINE
+ void createTwoColumns(
const const_monomial monoA1,
const const_monomial monoA2,
- const const_monomial monoB,
- QuadMatrixBuilder& builder
- );
- MATHICGB_INLINE const std::pair<LeftRightColIndex, LeftRightColIndex>
- findOrCreateTwoColumns(
- const const_monomial monoA1,
- const const_monomial monoA2,
- const const_monomial monoB,
- const ColReader& colMap,
- QuadMatrixBuilder& builder
+ const const_monomial monoB
);
- const int mThreadCount;
+ std::mutex mCreateColumnLock;
+ ColIndex mLeftColCount;
+ ColIndex mRightColCount;
monomial mTmp;
- std::vector<RowTask> mTodo;
+ const int mThreadCount;
const PolyBasis& mBasis;
+ Monomials mMonomialsLeft;
+ Monomials mMonomialsRight;
QuadMatrixBuilder mBuilder;
+ Map mMap;
+ std::vector<RowTask> mTodo;
};
#endif
diff --git a/src/mathicgb/FixedSizeMonomialMap.h b/src/mathicgb/FixedSizeMonomialMap.h
index 1166aaa..0991708 100755
--- a/src/mathicgb/FixedSizeMonomialMap.h
+++ b/src/mathicgb/FixedSizeMonomialMap.h
@@ -25,6 +25,9 @@ public:
typedef T mapped_type;
typedef std::pair<const_monomial, mapped_type> value_type;
+ /// Iterates through entries in the hash table.
+ class const_iterator;
+
// Construct a hash table with at least requestedBucketCount buckets. There
// may be more buckets. Currently the number is rounded up to the next power
// of two.
@@ -75,7 +78,7 @@ public:
const size_t index = hashToIndex(mRing.monomialHashValue(node->mono));
Node* const next = node->next.load();
node->next.store(mBuckets[index].load());
- mBuckets[index].store(node);
+ mBuckets[index].store(node, std::memory_order_relaxed);
node = next;
}
}
@@ -88,8 +91,26 @@ public:
const PolyRing& ring() const {return mRing;}
- // Returns the value associated to mono or null if there is no such value.
- const mapped_type* find(const const_monomial mono) const {
+ /// The range [begin(), end()) contains all entries in the hash table.
+ /// Insertions invalidate all iterators. Beware that insertions can
+ /// happen concurrently.
+ const_iterator begin() const {
+ const auto bucketsBegin = mBuckets.get();
+ const auto bucketsEnd = bucketsBegin + bucketCount();
+ return const_iterator(bucketsBegin, bucketsEnd);
+ }
+
+ const_iterator end() const {
+ const auto bucketsBegin = mBuckets.get();
+ const auto bucketsEnd = bucketsBegin + bucketCount();
+ return const_iterator(bucketsEnd, bucketsEnd);
+ }
+
+ /// Returns the value associated to mono or null if there is no such value.
+ /// Also returns an internal monomial that equals mono of such a monomial
+ /// exists.
+ std::pair<const mapped_type*, ConstMonomial>
+ find(const const_monomial mono) const {
const HashValue monoHash = mRing.monomialHashValue(mono);
const Node* node = bucketAtIndex(hashToIndex(monoHash));
for (; node != 0; node = node->next.load(std::memory_order_consume)) {
@@ -98,13 +119,15 @@ public:
//if (monoHash != mRing.monomialHashValue(node->mono))
// continue;
if (mRing.monomialEqualHintTrue(mono, node->mono))
- return &node->value;
+ return std::make_pair(&node->value, node->mono);
}
- return 0;
+ return std::make_pair(static_cast<const mapped_type*>(0), ConstMonomial(0));
}
- // As find on the product a*b.
- MATHICGB_INLINE const mapped_type* findProduct(
+ // As find on the product a*b but also returns the monomial that is the
+ // product.
+ MATHICGB_INLINE
+ std::pair<const mapped_type*, ConstMonomial> findProduct(
const const_monomial a,
const const_monomial b
) const {
@@ -116,9 +139,9 @@ public:
//if (abHash != mRing.monomialHashValue(node->mono))
// continue;
if (mRing.monomialIsProductOfHintTrue(a, b, node->mono))
- return &node->value;
+ return std::make_pair(&node->value, node->mono);
}
- return 0;
+ return std::make_pair(static_cast<const mapped_type*>(0), ConstMonomial(0));
}
/// As findProduct but looks for a1*b and a2*b at one time.
@@ -138,16 +161,19 @@ public:
)
return std::make_pair(&node1->value, &node2->value);
else
- return std::make_pair(findProduct(a1, b), findProduct(a2, b));
+ return std::make_pair(findProduct(a1, b).first, findProduct(a2, b).first);
}
/// Makes value.first map to value.second unless value.first is already
/// present in the map - in that case nothing is done. If p is the returned
- /// pair then *p.first is the value that value.first maps to after the insert
- /// and p.second is true if an insertion was performed. *p.first will not
+ /// pair then *p.first.first is the value that value.first maps to after the insert
+ /// and p.second is true if an insertion was performed. *p.first.first will not
/// equal value.second if an insertion was not performed - unless the
/// inserted value equals the already present value.
- std::pair<const mapped_type*, bool> insert(const value_type& value) {
+ ///
+ /// p.first.second is a internal monomial that equals value.first.
+ std::pair<std::pair<const mapped_type*, ConstMonomial>, bool>
+ insert(const value_type& value) {
const std::lock_guard<std::mutex> lockGuard(mInsertionMutex);
// find() loads buckets with memory_order_consume, so it may seem like
// we need some extra synchronization to make sure that we have the
@@ -158,8 +184,8 @@ public:
// insertion mutex and by locking that mutex we have synchronized with
// all threads that previously did insertions.
{
- const mapped_type* const found = find(value.first);
- if (found != 0)
+ const auto found = find(value.first);
+ if (found.first != 0)
return std::make_pair(found, false); // key already present
}
@@ -176,7 +202,7 @@ public:
// lock only synchronizes with threads who later grab the lock - it does
// not synchronize with reading threads since they do not grab the lock.
mBuckets[index].store(node, std::memory_order_release);
- return std::make_pair(&node->value, true); // successful insertion
+ return std::make_pair(std::make_pair(&node->value, node->mono), true); // successful insertion
}
/// This operation removes all entries from the table. This operation
@@ -264,6 +290,81 @@ private:
const PolyRing& mRing;
memt::BufferPool mNodeAlloc; // nodes are allocated from here.
std::mutex mInsertionMutex;
+
+public:
+ class const_iterator {
+ public:
+ const_iterator(): mNode(0), mBucket(0), mBucketEnd(0) {}
+
+ const_iterator& operator++() {
+ MATHICGB_ASSERT(mNode != 0);
+ MATHICGB_ASSERT(mBucket != mBucketsEnd);
+ const Node* const node = mNode->next.load(std::memory_order_consume);
+ if (node != 0)
+ mNode = node;
+ else
+ advanceBucket();
+ return *this;
+ }
+
+ bool operator==(const const_iterator& it) const {
+ MATHICGB_ASSERT(fromSameMap(it));
+ return mNode == it.mNode;
+ }
+
+ bool operator!=(const const_iterator& it) const {
+ return !(*this == it);
+ }
+
+ const std::pair<mapped_type, ConstMonomial> operator*() const {
+ MATHICGB_ASSERT(mNode != 0);
+ return std::make_pair(mNode->value, mNode->mono);
+ }
+
+ private:
+ friend class FixedSizeMonomialMap<T>;
+ const_iterator(
+ const Atomic<Node*>* const bucketBegin,
+ const Atomic<Node*>* const bucketEnd
+ ):
+ mBucket(bucketBegin),
+ mBucketsEnd(bucketEnd)
+ {
+ if (bucketBegin == bucketEnd) {
+ mNode = 0;
+ return;
+ }
+ const Node* const node = bucketBegin->load(std::memory_order_consume);
+ if (node != 0)
+ mNode = node;
+ else
+ advanceBucket();
+ }
+
+ void advanceBucket() {
+ MATHICGB_ASSERT(mBucket != mBucketsEnd);
+ while (true) {
+ ++mBucket;
+ if (mBucket == mBucketsEnd) {
+ mNode = 0;
+ break;
+ }
+ const Node* const node = mBucket->load(std::memory_order_consume);
+ if (node != 0) {
+ mNode = node;
+ break;
+ }
+ }
+ }
+
+ bool fromSameMap(const const_iterator& it) const {
+ return mBucketsEnd == it.mBucketsEnd;
+ }
+
+ const Node* mNode;
+ const Atomic<Node*>* mBucket;
+ const Atomic<Node*>* mBucketsEnd;
+ };
};
#endif
diff --git a/src/mathicgb/MonomialMap.hpp b/src/mathicgb/MonomialMap.hpp
index c975733..e8b14ab 100755
--- a/src/mathicgb/MonomialMap.hpp
+++ b/src/mathicgb/MonomialMap.hpp
@@ -85,14 +85,17 @@ public:
}
/// Returns the value that mono maps to or null if no such key has been
- /// inserted. Misses can be spurious! Read the comments on the parent
+ /// inserted. Also returns the internal monomial that matches mono if such
+ /// a monomial exists. Misses can be spurious! Read the comments on the parent
/// class.
- const mapped_type* find(const_monomial mono) const {
+ std::pair<const mapped_type*, ConstMonomial>
+ find(const_monomial mono) const {
return mMap.find(mono);
}
- // As find but looks for the product of a and b.
- const mapped_type* findProduct(
+ // As find but looks for the product of a and b and also returns the
+ // monomal that is the product.
+ std::pair<const mapped_type*, ConstMonomial> findProduct(
const const_monomial a,
const const_monomial b
) const {
@@ -111,6 +114,14 @@ public:
return mMap.findTwoProducts(a1, a2, b);
}
+ typedef typename FixedSizeMonomialMap<T>::const_iterator const_iterator;
+
+ /// The range [begin(), end()) contains all entries in the hash table.
+ /// Insertions invalidate all iterators. Beware that insertions can
+ /// happen concurrently.
+ const_iterator begin() const {return mMap.begin();}
+ const_iterator end() const {return mMap.end();}
+
private:
const FixedSizeMonomialMap<T>& mMap;
};
@@ -121,11 +132,13 @@ public:
/// Makes value.first map to value.second unless value.first is already
/// present in the map - in that case nothing is done. If p is the returned
- /// pair then *p.first is the value that value.first maps to after the insert
- /// and p.second is true if an insertion was performed. *p.first will not
+ /// pair then *p.first.first is the value that value.first maps to after the insert
+ /// and p.second is true if an insertion was performed. *p.first.first will not
/// equal value.second if an insertion was not performed - unless the
- /// inserted value equals the already present value.
- std::pair<const mapped_type*, bool> insert(const value_type& value) {
+ /// inserted value equals the already present value. p.first.second is an
+ /// internal monomial that equals value.first.
+ std::pair<std::pair<const mapped_type*, ConstMonomial>, bool>
+ insert(const value_type& value) {
const std::lock_guard<std::mutex> lockGuard(mInsertionMutex);
// We can load mMap as std::memory_order_relaxed because we have already
diff --git a/src/mathicgb/QuadMatrix.cpp b/src/mathicgb/QuadMatrix.cpp
index bd57bb5..c049de2 100755
--- a/src/mathicgb/QuadMatrix.cpp
+++ b/src/mathicgb/QuadMatrix.cpp
@@ -9,13 +9,16 @@ bool QuadMatrix::debugAssertValid() const {
#ifndef MATHICGB_DEBUG
return true;
#else
- MATHICGB_ASSERT(topLeft.computeColCount() <= leftColumnMonomials.size());
- MATHICGB_ASSERT(bottomLeft.computeColCount() <= leftColumnMonomials.size());
- MATHICGB_ASSERT(topLeft.rowCount() == topRight.rowCount());
-
- MATHICGB_ASSERT(topRight.computeColCount() <= rightColumnMonomials.size());
- MATHICGB_ASSERT(bottomRight.computeColCount() <=
+ if (!leftColumnMonomials.empty() || !rightColumnMonomials.empty()) {
+ MATHICGB_ASSERT(topLeft.computeColCount() <= leftColumnMonomials.size());
+ MATHICGB_ASSERT(bottomLeft.computeColCount() <=
+ leftColumnMonomials.size());
+ MATHICGB_ASSERT(topRight.computeColCount() <= rightColumnMonomials.size());
+ MATHICGB_ASSERT(bottomRight.computeColCount() <=
rightColumnMonomials.size());
+ }
+
+ MATHICGB_ASSERT(topLeft.rowCount() == topRight.rowCount());
MATHICGB_ASSERT(bottomRight.rowCount() == bottomLeft.rowCount());
return true;
#endif
@@ -215,3 +218,87 @@ std::ostream& operator<<(std::ostream& out, const QuadMatrix& qm) {
qm.print(out);
return out;
}
+
+namespace {
+ class ColumnComparer {
+ public:
+ ColumnComparer(const PolyRing& ring): mRing(ring) {}
+
+ typedef SparseMatrix::ColIndex ColIndex;
+ typedef std::pair<monomial, ColIndex> Pair;
+ bool operator()(const Pair& a, const Pair b) const {
+ return mRing.monomialLT(b.first, a.first);
+ }
+
+ private:
+ const PolyRing& mRing;
+ };
+
+ std::vector<SparseMatrix::ColIndex> sortColumnMonomialsAndMakePermutation(
+ std::vector<monomial>& monomials,
+ const PolyRing& ring
+ ) {
+ typedef SparseMatrix::ColIndex ColIndex;
+ MATHICGB_ASSERT(monomials.size() <= std::numeric_limits<ColIndex>::max());
+ const ColIndex colCount = static_cast<ColIndex>(monomials.size());
+ // Monomial needs to be non-const as we are going to put these
+ // monomials back into the vector of monomials which is not const.
+ std::vector<std::pair<monomial, ColIndex>> columns;
+ columns.reserve(colCount);
+ for (ColIndex col = 0; col < colCount; ++col)
+ columns.push_back(std::make_pair(monomials[col], col));
+ std::sort(columns.begin(), columns.end(), ColumnComparer(ring));
+
+ // Apply sorting permutation to monomials. This is why it is necessary to
+ // copy the values in monomial out of there: in-place application of a
+ // permutation is messy.
+ MATHICGB_ASSERT(columns.size() == colCount);
+ MATHICGB_ASSERT(monomials.size() == colCount);
+ for (size_t col = 0; col < colCount; ++col) {
+ MATHICGB_ASSERT(col == 0 ||
+ ring.monomialLT(columns[col].first, columns[col - 1].first));
+ monomials[col] = columns[col].first;
+ }
+
+ // Construct permutation of indices to match permutation of monomials
+ std::vector<ColIndex> permutation(colCount);
+ for (ColIndex col = 0; col < colCount; ++col) {
+ // The monomial for column columns[col].second is now the
+ // monomial for col, so we need the inverse map for indices.
+ permutation[columns[col].second] = col;
+ }
+
+ return std::move(permutation);
+ }
+}
+
+void QuadMatrix::sortColumnsLeftRightParallel(const int threadCount) {
+ typedef SparseMatrix::ColIndex ColIndex;
+ std::vector<ColIndex> leftPermutation;
+ std::vector<ColIndex> rightPermutation;
+
+#pragma omp parallel for num_threads(threadCount) schedule(static)
+ for (OMPIndex i = 0; i < 2; ++i) {
+ if (i == 0)
+ leftPermutation =
+ sortColumnMonomialsAndMakePermutation(leftColumnMonomials, *ring);
+ else
+ rightPermutation =
+ sortColumnMonomialsAndMakePermutation(rightColumnMonomials, *ring);
+ }
+
+ // todo: parallelize per block instead of per matrix.
+#pragma omp parallel for num_threads(threadCount) schedule(dynamic)
+ for (OMPIndex i = 0; i < 4; ++i) {
+ if (i == 0)
+ topRight.applyColumnMap(rightPermutation);
+ else if (i == 1)
+ bottomRight.applyColumnMap(rightPermutation);
+ else if (i == 2)
+ topLeft.applyColumnMap(leftPermutation);
+ else {
+ MATHICGB_ASSERT(i == 3);
+ bottomLeft.applyColumnMap(leftPermutation);
+ }
+ }
+}
diff --git a/src/mathicgb/QuadMatrix.hpp b/src/mathicgb/QuadMatrix.hpp
index a508114..a66d316 100755
--- a/src/mathicgb/QuadMatrix.hpp
+++ b/src/mathicgb/QuadMatrix.hpp
@@ -55,6 +55,11 @@ public:
/// Shows whole matrix in a string. Useful for debugging.
std::string toString() const;
+ /// Sort the left columns to be in decreasing order according to the monomial
+ /// order from the ring. The operation is parallel using up to threadCount
+ /// threads.
+ void sortColumnsLeftRightParallel(int threadCount);
+
/// Makes a copy of this matrix whose rows are sorted in some canonical way.
/// TODO: Actually only coarsely sorts the top rows right now.
QuadMatrix toCanonical() const;
diff --git a/src/mathicgb/QuadMatrixBuilder.cpp b/src/mathicgb/QuadMatrixBuilder.cpp
index 5401221..93ddaa1 100755
--- a/src/mathicgb/QuadMatrixBuilder.cpp
+++ b/src/mathicgb/QuadMatrixBuilder.cpp
@@ -8,39 +8,23 @@
QuadMatrixBuilder::QuadMatrixBuilder(
const PolyRing& ring,
+ Map& map,
+ MonomialsType& monomialsLeft,
+ MonomialsType& monomialsRight,
const size_t memoryQuantum
):
- mMonomialToCol(ring),
+ mMonomialToCol(map),
mTopLeft(memoryQuantum),
mTopRight(memoryQuantum),
mBottomLeft(memoryQuantum),
- mBottomRight(memoryQuantum)
+ mBottomRight(memoryQuantum),
+ mMonomialsLeft(monomialsLeft),
+ mMonomialsRight(monomialsRight)
{}
void QuadMatrixBuilder::takeRowsFrom(QuadMatrix&& matrix) {
MATHICGB_ASSERT(&ring() == matrix.ring);
MATHICGB_ASSERT(matrix.debugAssertValid());
-#ifdef MATHICGB_DEBUG
- if (!matrix.leftColumnMonomials.empty() ||
- !matrix.rightColumnMonomials.empty()
- ) {
- // check left column monomials are the same
- for (ColIndex col = 0; col < matrix.leftColumnMonomials.size(); ++col) {
- MATHICGB_ASSERT(ring().monomialEQ
- (matrix.leftColumnMonomials[col], monomialOfLeftCol(col)));
- }
-
- // check right column monomials are the same
- for (ColIndex col = 0; col < matrix.rightColumnMonomials.size(); ++col) {
- MATHICGB_ASSERT(ring().monomialEQ
- (matrix.rightColumnMonomials[col], monomialOfRightCol(col)));
- }
- }
-#endif
-
- matrix.ring = 0;
- matrix.leftColumnMonomials.clear();
- matrix.rightColumnMonomials.clear();
mTopLeft.takeRowsFrom(std::move(matrix.topLeft));
mTopRight.takeRowsFrom(std::move(matrix.topRight));
@@ -56,16 +40,17 @@ namespace {
/// template in order to avoid referring to private types of
/// QuadMatrixBuilder.
template<class ToMono, class ToCol>
- QuadMatrixBuilder::LeftRightColIndex createCol
- (const_monomial mono,
- SparseMatrix& top,
- SparseMatrix& bottom,
- ToMono& toMonomial,
- ToCol& toCol,
- const PolyRing& ring,
- const bool left)
- {
- MATHICGB_ASSERT(typename ToCol::Reader(toCol).find(mono) == 0);
+ std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+ createCol(
+ const_monomial mono,
+ SparseMatrix& top,
+ SparseMatrix& bottom,
+ ToMono& toMonomial,
+ ToCol& toCol,
+ const PolyRing& ring,
+ const bool left
+ ) {
+ MATHICGB_ASSERT(typename ToCol::Reader(toCol).find(mono).first == 0);
const auto colCount =
static_cast<QuadMatrixBuilder::ColIndex>(toMonomial.size());
@@ -75,23 +60,31 @@ namespace {
toMonomial.push_back(0); // allocate memory now to avoid bad_alloc later
monomial copied = ring.allocMonomial();
ring.monomialCopy(mono, copied);
+ std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial> p;
try {
- toCol.insert(std::make_pair
- (copied,
- QuadMatrixBuilder::LeftRightColIndex(colCount, left)));
+ auto inserted = toCol.insert(std::make_pair(
+ copied, QuadMatrixBuilder::LeftRightColIndex(colCount, left))
+ );
+ MATHICGB_ASSERT(inserted.second);
+ MATHICGB_ASSERT(inserted.first.first != 0);
+ auto p(inserted.first);
+ toMonomial.back() = copied;
+
+ MATHICGB_ASSERT(ring.monomialEqualHintTrue(copied, p.second));
+ MATHICGB_ASSERT(*p.first == QuadMatrixBuilder::LeftRightColIndex(colCount, left));
+ return std::make_pair(*p.first, p.second);
} catch (...) {
toMonomial.pop_back();
ring.freeMonomial(copied);
throw;
}
- toMonomial.back() = copied;
-
- return QuadMatrixBuilder::LeftRightColIndex(colCount, left);
}
}
-QuadMatrixBuilder::LeftRightColIndex QuadMatrixBuilder::createColumnLeft
-(const_monomial monomialToBeCopied) {
+std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+QuadMatrixBuilder::createColumnLeft(
+ const_monomial monomialToBeCopied
+) {
return createCol
(monomialToBeCopied,
mTopLeft,
@@ -102,8 +95,10 @@ QuadMatrixBuilder::LeftRightColIndex QuadMatrixBuilder::createColumnLeft
true);
}
-QuadMatrixBuilder::LeftRightColIndex QuadMatrixBuilder::createColumnRight
-(const_monomial monomialToBeCopied) {
+std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+QuadMatrixBuilder::createColumnRight(
+ const_monomial monomialToBeCopied
+) {
return createCol
(monomialToBeCopied,
mTopRight,
@@ -114,183 +109,6 @@ QuadMatrixBuilder::LeftRightColIndex QuadMatrixBuilder::createColumnRight
false);
}
-namespace {
- class ColumnComparer {
- public:
- ColumnComparer(const FreeModuleOrder& order): mOrder(order) {}
-
- typedef QuadMatrixBuilder::ColIndex ColIndex;
- typedef std::pair<monomial, ColIndex> Pair;
- bool operator()(const Pair& a, const Pair b) const {
- return mOrder.signatureCompare(a.first, b.first) == GT;
- }
-
- private:
- const FreeModuleOrder& mOrder;
- };
-
- std::vector<SparseMatrix::ColIndex> sortColumnMonomialsAndMakePermutation(
- const FreeModuleOrder& order,
- std::vector<monomial>& monomials
- ) {
- typedef SparseMatrix::ColIndex ColIndex;
- MATHICGB_ASSERT(monomials.size() <= std::numeric_limits<ColIndex>::max());
- const ColIndex colCount = static_cast<ColIndex>(monomials.size());
- // Monomial needs to be non-const as we are going to put these
- // monomials back into the vector of monomials which is not const.
- std::vector<std::pair<monomial, ColIndex>> columns;
- columns.reserve(colCount);
- for (ColIndex col = 0; col < colCount; ++col)
- columns.push_back(std::make_pair(monomials[col], col));
- std::sort(columns.begin(), columns.end(), ColumnComparer(order));
-
- // Apply sorting permutation to monomials. This is why it is necessary to
- // copy the values in monomial out of there: in-place application of a
- // permutation is messy.
- MATHICGB_ASSERT(columns.size() == colCount);
- MATHICGB_ASSERT(monomials.size() == colCount);
- for (size_t col = 0; col < colCount; ++col) {
- MATHICGB_ASSERT(col == 0 || order.signatureCompare
- (columns[col - 1].first, columns[col].first) == GT);
- monomials[col] = columns[col].first;
- }
-
- // Construct permutation of indices to match permutation of monomials
- std::vector<ColIndex> permutation(colCount);
- for (ColIndex col = 0; col < colCount; ++col) {
- // The monomial for column columns[col].second is now the
- // monomial for col, so we need the inverse map for indices.
- permutation[columns[col].second] = col;
- }
-
- return std::move(permutation);
- }
-
-
-
- // The purpose of this function is to avoid code duplication for
- // left/right variants.
- void sortColumns(
- const FreeModuleOrder& order,
- std::vector<monomial>& monomials,
- SparseMatrix& topMatrix,
- SparseMatrix& bottomMatrix
- ) {
- typedef SparseMatrix::ColIndex ColIndex;
- const auto colCount = static_cast<ColIndex>(monomials.size());
-
- // Monomial needs to be non-const as we are going to put these
- // monomials back into the vector of monomials which is not const.
- std::vector<std::pair<monomial, ColIndex> > columns;
- columns.reserve(colCount);
- for (ColIndex col = 0; col < colCount; ++col)
- columns.push_back(std::make_pair(monomials[col], col));
- std::sort(columns.begin(), columns.end(), ColumnComparer(order));
-
- // Reorder monomials associated to each column to be in sorted order
- MATHICGB_ASSERT(columns.size() == colCount);
- MATHICGB_ASSERT(monomials.size() == colCount);
- for (size_t col = 0; col < colCount; ++col) {
- MATHICGB_ASSERT(col == 0 || order.signatureCompare
- (columns[col - 1].first, columns[col].first) == GT);
- monomials[col] = columns[col].first;
- }
-
- // Construct permutation of indices to match permutation of monomials
- std::vector<ColIndex> permutation(colCount);
- for (ColIndex col = 0; col < colCount; ++col) {
- // The monomial for column columns[col].second is now the
- // monomial for col, so we need the inverse map for indices.
- permutation[columns[col].second] = col;
- }
-
- // Change indices to match the permutation done on the columns
- topMatrix.applyColumnMap(permutation);
- bottomMatrix.applyColumnMap(permutation);
- }
-}
-
-void QuadMatrixBuilder::sortColumnsLeft(const FreeModuleOrder& order) {
- sortColumns(order, mMonomialsLeft, mTopLeft, mBottomLeft);
-}
-
-void QuadMatrixBuilder::sortColumnsRight(const FreeModuleOrder& order) {
- sortColumns(order, mMonomialsRight, mTopRight, mBottomRight);
-}
-
-void QuadMatrixBuilder::sortColumnsLeftRightParallel(
- const FreeModuleOrder& order,
- int threadCount
-) {
- std::vector<ColIndex> leftPermutation;
- std::vector<ColIndex> rightPermutation;
-
-#pragma omp parallel for num_threads(threadCount) schedule(static)
- for (OMPIndex i = 0; i < 2; ++i) {
- if (i == 0)
- leftPermutation =
- sortColumnMonomialsAndMakePermutation(order, mMonomialsLeft);
- else
- rightPermutation =
- sortColumnMonomialsAndMakePermutation(order, mMonomialsRight);
- }
-
-#pragma omp parallel for num_threads(threadCount) schedule(dynamic)
- for (OMPIndex i = 0; i < 4; ++i) {
- if (i == 0)
- mTopRight.applyColumnMap(rightPermutation);
- else if (i == 1)
- mBottomRight.applyColumnMap(rightPermutation);
- else if (i == 2)
- mTopLeft.applyColumnMap(leftPermutation);
- else {
- MATHICGB_ASSERT(i == 3);
- mBottomLeft.applyColumnMap(leftPermutation);
- }
- }
-}
-
-void QuadMatrixBuilder::print(std::ostream& out) const {
- mathic::ColumnPrinter printer;
- printer.addColumn(true, "", "");
- printer.addColumn(true, " | ", "");
-
- // column monomials
- out << "Left columns:";
- const auto leftColCount = static_cast<ColIndex>(mMonomialsLeft.size());
- for (ColIndex leftCol = 0; leftCol < leftColCount; ++leftCol) {
- out << ' ';
- ring().monomialDisplay(out, monomialOfLeftCol(leftCol), false, true);
- }
-
- out << "\nRight columns:";
- const auto rightColCount = static_cast<ColIndex>(mMonomialsRight.size());
- for (ColIndex rightCol = 0; rightCol < rightColCount; ++rightCol) {
- out << ' ';
- ring().monomialDisplay(out, monomialOfRightCol(rightCol), false, true);
- }
- out << '\n';
-
- // left side
- topLeft().print(printer[0]);
- printer[0] << '\n';
- bottomLeft().print(printer[0]);
-
- // right side
- topRight().print(printer[1]);
- printer[1] << '\n';
- bottomRight().print(printer[1]);
-
- out << printer;
-}
-
-// String representation intended for debugging.
-std::string QuadMatrixBuilder::toString() const {
- std::ostringstream out;
- print(out);
- return out.str();
-}
-
QuadMatrix QuadMatrixBuilder::buildMatrixAndClear() {
QuadMatrix out;
@@ -298,18 +116,12 @@ QuadMatrix QuadMatrixBuilder::buildMatrixAndClear() {
mTopRight.swap(out.topRight);
mBottomLeft.swap(out.bottomLeft);
mBottomRight.swap(out.bottomRight);
- mMonomialsLeft.swap(out.leftColumnMonomials);
- mMonomialsRight.swap(out.rightColumnMonomials);
out.ring = &ring();
mTopLeft.clear();
mTopRight.clear();
mBottomLeft.clear();
mBottomRight.clear();
- mMonomialsLeft.clear();
- mMonomialsRight.clear();
-
- mMonomialToCol.clearNonConcurrent();
MATHICGB_ASSERT(out.debugAssertValid());
return std::move(out);
diff --git a/src/mathicgb/QuadMatrixBuilder.hpp b/src/mathicgb/QuadMatrixBuilder.hpp
index 65451df..3794fe2 100755
--- a/src/mathicgb/QuadMatrixBuilder.hpp
+++ b/src/mathicgb/QuadMatrixBuilder.hpp
@@ -29,14 +29,6 @@ class QuadMatrixBuilder {
typedef SparseMatrix::ColIndex ColIndex;
typedef SparseMatrix::Scalar Scalar;
- QuadMatrixBuilder(const PolyRing& ring, size_t memoryQuantum = 0);
-
- /// Inserts the rows from builder. To avoid an assert either the matrix must
- /// have no column monomials specified or the monomials that are specified
- /// must match exactly to the column monomials for this object --- including
- /// the ordering of the monomials.
- void takeRowsFrom(QuadMatrix&& matrix);
-
/// The index of a column that can be either on the left or the
/// right side. The largest representable ColIndex is an invalid
/// index. This is the default value. The only allowed method to
@@ -93,14 +85,26 @@ class QuadMatrixBuilder {
bool mLeft;
};
- ColIndex leftColCount() const {
- return static_cast<ColIndex>(mMonomialsLeft.size());
- }
+ typedef MonomialMap<LeftRightColIndex> Map;
+ typedef std::vector<monomial> MonomialsType;
- ColIndex rightColCount() const {
- return static_cast<ColIndex>(mMonomialsRight.size());
- }
+ QuadMatrixBuilder(
+ const PolyRing& ring,
+ Map& map,
+ MonomialsType& monomialsLeft,
+ MonomialsType& monomialsRight,
+ size_t memoryQuantum = 0
+ );
+
+ /// Inserts the rows from builder. To avoid an assert either the matrix must
+ /// have no column monomials specified or the monomials that are specified
+ /// must match exactly to the column monomials for this object --- including
+ /// the ordering of the monomials.
+ void takeRowsFrom(QuadMatrix&& matrix);
+ size_t memoryQuantum() const {
+ return mTopLeft.memoryQuantum();
+ }
// **** Appending entries to top matrices.
// Same interface as SparseMatrix except with two matrices and here
@@ -163,102 +167,22 @@ class QuadMatrixBuilder {
/** Creates a new column associated to the monomial
monomialToBeCopied to the left matrices. There must not already
exist a column for this monomial on the left or on the right. */
- LeftRightColIndex createColumnLeft(const_monomial monomialToBeCopied);
+ std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+ createColumnLeft(const_monomial monomialToBeCopied);
/** Creates a new column associated to the monomial monomialToBeCopied
to the right matrices. There must not already exist a column for
this monomial on the left or on the right. */
- LeftRightColIndex createColumnRight(const_monomial monomialToBeCopied);
-
- /// As calling sortColumnsLeft() and sortColumnsRight(), but performs
- /// the operations in parallel using up to threadCount threads.
- void sortColumnsLeftRightParallel
- (const FreeModuleOrder& order, int threadCount);
-
- /** Sorts the left columns to be decreasing with respect to
- order. Also updates the column indices already in the matrix to
- reflect the new ordering. */
- void sortColumnsLeft(const FreeModuleOrder& order);
-
- /** Sorts the right columns to be decreasing with respect to
- order. Also updates the column indices already in the matrix to
- reflect the new ordering. */
- void sortColumnsRight(const FreeModuleOrder& order);
-
+ std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+ createColumnRight(const_monomial monomialToBeCopied);
// *** Querying columns
- typedef const MonomialMap<LeftRightColIndex>::Reader ColReader;
- const MonomialMap<LeftRightColIndex>& columnToIndexMap() const {
- return mMonomialToCol;
- }
-
- /** Returns a column for the findThis monomial. Searches on both the
- left and right side. Returns an invalid index if no such column
- exists. */
- LeftRightColIndex findColumn(const_monomial findThis) const {
- auto it = ColReader(mMonomialToCol).find(findThis);
- if (it != 0)
- return *it;
- else
- return LeftRightColIndex();
- }
-
- /// As findColumn(), but looks for a*b. This is faster than computing a*b
- /// and then looking that up.
- LeftRightColIndex findColumnProduct(
- const const_monomial a,
- const const_monomial b
- ) const {
- const auto it = ColReader(mMonomialToCol).findProduct(a, b);
- if (it != 0)
- return *it;
- else
- return LeftRightColIndex();
- }
-
- /// As findColumnProduct(), but looks for a1*b and a2*b at the same time.
- MATHICGB_INLINE std::pair<LeftRightColIndex, LeftRightColIndex>
- findTwoColumnsProduct(
- const const_monomial a1,
- const const_monomial a2,
- const const_monomial b
- ) const {
- const auto it = ColReader(mMonomialToCol).findTwoProducts(a1, a2, b);
- return std::make_pair(
- it.first != 0 ? *it.first : LeftRightColIndex(),
- it.second != 0 ? *it.second : LeftRightColIndex());
- }
-
- const_monomial monomialOfLeftCol(ColIndex col) const {
- MATHICGB_ASSERT(col < mMonomialsLeft.size());
- return mMonomialsLeft[col];
- }
-
- const_monomial monomialOfRightCol(ColIndex col) const {
- MATHICGB_ASSERT(col < mMonomialsRight.size());
- return mMonomialsRight[col];
- }
-
- const_monomial monomialOfCol(LeftRightColIndex col) const {
- MATHICGB_ASSERT(col.valid());
- if (col.left())
- return monomialOfLeftCol(col.leftIndex());
- else
- return monomialOfRightCol(col.rightIndex());
- }
-
const SparseMatrix& topLeft() const {return mTopLeft;}
const SparseMatrix& topRight() const {return mTopRight;}
const SparseMatrix& bottomLeft() const {return mBottomLeft;}
const SparseMatrix& bottomRight() const {return mBottomRight;}
- // String representation intended for debugging.
- void print(std::ostream& out) const;
-
- // String representation intended for debugging.
- std::string toString() const;
-
const PolyRing& ring() const {return mMonomialToCol.ring();}
/// Returns the built matrix and sets the builder to a state
@@ -266,12 +190,11 @@ class QuadMatrixBuilder {
QuadMatrix buildMatrixAndClear();
private:
- typedef std::vector<monomial> MonomialsType;
- MonomialsType mMonomialsLeft; /// stores one monomial per left column
- MonomialsType mMonomialsRight; /// stores one monomial per right column
+ MonomialsType& mMonomialsLeft; /// stores one monomial per left column
+ MonomialsType& mMonomialsRight; /// stores one monomial per right column
/// Used for fast determination of which column has a given monomial.
- MonomialMap<LeftRightColIndex> mMonomialToCol;
+ Map& mMonomialToCol;
SparseMatrix mTopLeft;
SparseMatrix mTopRight;
diff --git a/src/test/QuadMatrixBuilder.cpp b/src/test/QuadMatrixBuilder.cpp
index 04d31c4..ee4f828 100755
--- a/src/test/QuadMatrixBuilder.cpp
+++ b/src/test/QuadMatrixBuilder.cpp
@@ -26,20 +26,14 @@ namespace {
size_t colCount = 0;
for (Poly::iterator it = p.begin(); it != p.end(); ++it) {
QuadMatrixBuilder::LeftRightColIndex lrCol =
- b.createColumnLeft(it.getMonomial());
+ b.createColumnLeft(it.getMonomial()).first;
ASSERT_TRUE(lrCol.left());
ASSERT_FALSE(lrCol.right());
auto col = lrCol.leftIndex();
ASSERT_EQ(col, lrCol.index());
ASSERT_EQ(colCount, col);
++colCount;
- // not equal as pointers
- ASSERT_TRUE(it.getMonomial().unsafeGetRepresentation() !=
- b.monomialOfLeftCol(col).unsafeGetRepresentation());
- ASSERT_TRUE // equal as values
- (b.ring().monomialEQ(it.getMonomial(), b.monomialOfLeftCol(col)));
}
- ASSERT_EQ(colCount, b.leftColCount());
}
{
Poly p(b.ring());
@@ -48,39 +42,43 @@ namespace {
size_t colCount = 0;
for (Poly::iterator it = p.begin(); it != p.end(); ++it) {
QuadMatrixBuilder::LeftRightColIndex lrCol =
- b.createColumnRight(it.getMonomial());
+ b.createColumnRight(it.getMonomial()).first;
ASSERT_TRUE(lrCol.right());
ASSERT_FALSE(lrCol.left());
auto col = lrCol.rightIndex();
ASSERT_EQ(col, lrCol.index());
ASSERT_EQ(colCount, col);
++colCount;
- // not equal as pointers
- ASSERT_TRUE(it.getMonomial().unsafeGetRepresentation() !=
- b.monomialOfRightCol(col).unsafeGetRepresentation());
- ASSERT_TRUE // equal as values
- (b.ring().monomialEQ(it.getMonomial(), b.monomialOfRightCol(col)));
}
- ASSERT_EQ(colCount, b.rightColCount());
}
}
}
TEST(QuadMatrixBuilder, Empty) {
+ // test a builder with no rows and no columns
PolyRing ring(2, 0, 0);
- QuadMatrixBuilder b(ring); // test a builder with no rows and no columns
+ QuadMatrixBuilder::Map map(ring);
+ QuadMatrixBuilder::MonomialsType monoLeft;
+ QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder b(ring, map, monoLeft, monoRight);
const char* matrixStr =
"Left columns:\n"
"Right columns:\n"
"matrix with no rows | matrix with no rows\n"
" | \n"
"matrix with no rows | matrix with no rows\n";
- ASSERT_EQ(matrixStr, b.toString());
+ auto matrix = b.buildMatrixAndClear();
+ matrix.leftColumnMonomials = monoLeft;
+ matrix.rightColumnMonomials = monoRight;
+ ASSERT_EQ(matrixStr, matrix.toString());
}
TEST(QuadMatrixBuilder, Construction) {
std::unique_ptr<PolyRing> ring(ringFromString("32003 6 1\n1 1 1 1 1 1"));
- QuadMatrixBuilder b(*ring);
+ QuadMatrixBuilder::Map map(*ring);
+ QuadMatrixBuilder::MonomialsType monoLeft;
+ QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
createColumns("a<1>+<0>", "bc<0>+b<0>+c<0>", b);
// top row: nothing, nothing
@@ -112,12 +110,18 @@ TEST(QuadMatrixBuilder, Construction) {
"0: 1#4 | 0: \n"
"1: | 1: 0#5\n"
"2: | 2: \n";
- ASSERT_EQ(matrixStr, b.toString());
+ auto matrix = b.buildMatrixAndClear();
+ matrix.leftColumnMonomials = monoLeft;
+ matrix.rightColumnMonomials = monoRight;
+ ASSERT_EQ(matrixStr, matrix.toString());
}
TEST(QuadMatrixBuilder, ColumnQuery) {
std::unique_ptr<PolyRing> ring(ringFromString("32003 6 1\n1 1 1 1 1 1"));
- QuadMatrixBuilder b(*ring);
+ QuadMatrixBuilder::Map map(*ring);
+ QuadMatrixBuilder::MonomialsType monoLeft;
+ QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
createColumns("a<1>+<0>", "b<0>+c<0>+bc<0>", b);
Poly p(b.ring());
@@ -126,16 +130,17 @@ TEST(QuadMatrixBuilder, ColumnQuery) {
("10a<1>+11<0>+20b<0>+21c<0>+22bc<0>+30ab<0>+30e<0>+10a<1>");
p.parseDoNotOrder(in);
for (Poly::iterator it = p.begin(); it != p.end(); ++it) {
- QuadMatrixBuilder::LeftRightColIndex col = b.findColumn(it.getMonomial());
+ const QuadMatrixBuilder::LeftRightColIndex* col =
+ MonomialMap<QuadMatrixBuilder::LeftRightColIndex>::Reader(map).find(it.getMonomial()).first;
if (it.getCoefficient() / 10 == 3)
- ASSERT_FALSE(col.valid());
+ ASSERT_EQ(col, static_cast<void*>(0));
else {
- ASSERT_TRUE(col.valid());
- ASSERT_EQ(it.getCoefficient() % 10, col.index());
+ ASSERT_TRUE(col != static_cast<void*>(0));
+ ASSERT_EQ(it.getCoefficient() % 10, col->index());
if (it.getCoefficient() / 10 == 2)
- ASSERT_TRUE(col.right());
+ ASSERT_TRUE(col->right());
else
- ASSERT_TRUE(col.left());
+ ASSERT_TRUE(col->left());
}
}
}
@@ -148,21 +153,27 @@ TEST(QuadMatrixBuilder, SortColumns) {
// one row top, no rows bottom, no columns
{
- QuadMatrixBuilder b(*ring);
+ QuadMatrixBuilder::Map map(*ring);
+ QuadMatrixBuilder::MonomialsType monoLeft;
+ QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
b.rowDoneTopLeftAndRight();
- b.sortColumnsLeft(*order);
- b.sortColumnsRight(*order);
+ auto matrix = b.buildMatrixAndClear();
+ matrix.sortColumnsLeftRightParallel(1);
const char* matrixStr =
"Left columns:\n"
"Right columns:\n"
"0: | 0: \n"
" | \n"
"matrix with no rows | matrix with no rows\n";
- ASSERT_EQ(matrixStr, b.toString());
+ ASSERT_EQ(matrixStr, matrix.toString());
}
{
- QuadMatrixBuilder b(*ring);
+ QuadMatrixBuilder::Map map(*ring);
+ QuadMatrixBuilder::MonomialsType monoLeft;
+ QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
createColumns("<0>+a<0>", "b<0>+bcd<0>+bc<0>", b);
b.appendEntryTopLeft(0,1);
b.appendEntryTopLeft(1,2);
@@ -177,6 +188,9 @@ TEST(QuadMatrixBuilder, SortColumns) {
b.appendEntryBottomRight(1,9);
b.appendEntryBottomRight(2,10);
b.rowDoneBottomLeftAndRight();
+ auto matrix = b.buildMatrixAndClear();
+ matrix.leftColumnMonomials = monoLeft;
+ matrix.rightColumnMonomials = monoRight;
const char* matrixStr1 =
"Left columns: 1 a\n"
@@ -184,40 +198,28 @@ TEST(QuadMatrixBuilder, SortColumns) {
"0: 0#1 1#2 | 0: 0#3 1#4 2#5 \n"
" | \n"
"0: 0#6 1#7 | 0: 0#8 1#9 2#10\n";
- ASSERT_EQ(matrixStr1, b.toString());
+ ASSERT_EQ(matrixStr1, matrix.toString());
const char* matrixStr2 =
"Left columns: a 1\n"
- "Right columns: b bcd bc\n"
- "0: 1#1 0#2 | 0: 0#3 1#4 2#5 \n"
- " | \n"
- "0: 1#6 0#7 | 0: 0#8 1#9 2#10\n";
- b.sortColumnsLeft(*order);
- ASSERT_EQ(matrixStr2, b.toString());
-
- b.sortColumnsLeft(*order); // sort when already sorted
- ASSERT_EQ(matrixStr2, b.toString());
-
- const char* matrixStr3 =
- "Left columns: a 1\n"
"Right columns: bcd bc b\n"
"0: 1#1 0#2 | 0: 2#3 0#4 1#5 \n"
" | \n"
"0: 1#6 0#7 | 0: 2#8 0#9 1#10\n";
- b.sortColumnsRight(*order);
- ASSERT_EQ(matrixStr3, b.toString());
-
- b.sortColumnsRight(*order); // sort when already sorted
- ASSERT_EQ(matrixStr3, b.toString());
+ matrix.sortColumnsLeftRightParallel(1);
+ ASSERT_EQ(matrixStr2, matrix.toString());
- b.sortColumnsLeft(*order);
- ASSERT_EQ(matrixStr3, b.toString());
+ matrix.sortColumnsLeftRightParallel(1);
+ ASSERT_EQ(matrixStr2, matrix.toString());
}
}
TEST(QuadMatrixBuilder, BuildAndClear) {
std::unique_ptr<PolyRing> ring(ringFromString("32003 6 1\n1 1 1 1 1 1"));
- QuadMatrixBuilder b(*ring);
+ QuadMatrixBuilder::Map map(*ring);
+ QuadMatrixBuilder::MonomialsType monoLeft;
+ QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
createColumns("a<1>+<0>", "b<0>+c<0>+bc<0>", b);
b.appendEntryTopLeft(1, 1);
@@ -229,18 +231,8 @@ TEST(QuadMatrixBuilder, BuildAndClear) {
b.rowDoneBottomLeftAndRight();
QuadMatrix qm(b.buildMatrixAndClear());
-
- // test that the matrix was really cleared
- ASSERT_EQ(ring.get(), &b.ring()); // still same ring though
- const char* matrixStr =
- "Left columns:\n"
- "Right columns:\n"
- "matrix with no rows | matrix with no rows\n"
- " | \n"
- "matrix with no rows | matrix with no rows\n";
- ASSERT_EQ(matrixStr, b.toString());
- ASSERT_EQ(0, b.leftColCount());
- ASSERT_EQ(0, b.rightColCount());
+ qm.leftColumnMonomials = monoLeft;
+ qm.rightColumnMonomials = monoRight;
// test that the quad matrix is right
ASSERT_EQ("0: 1#1\n", qm.topLeft.toString());
diff --git a/src/test/SparseMatrix.cpp b/src/test/SparseMatrix.cpp
index 97ba272..ab3f299 100755
--- a/src/test/SparseMatrix.cpp
+++ b/src/test/SparseMatrix.cpp
@@ -51,7 +51,7 @@ TEST(SparseMatrix, Simple) {
mat.rowDone(); // add a row with two entries
ASSERT_EQ(3, mat.entryCount());
ASSERT_EQ(3, mat.rowCount());
- ASSERT_EQ(6, mat.computeColCount());
+ ASSERT_EQ(2002, mat.computeColCount());
ASSERT_EQ(5, mat.leadCol(2));
ASSERT_EQ(2, mat.entryCountInRow(2));
ASSERT_EQ("0: 5#101\n1:\n2: 5#102 2001#0\n", mat.toString());
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/mathicgb.git
More information about the debian-science-commits
mailing list