[mathicgb] 144/393: Moved reducer 26 to using preblocks (no scalars) for the first version of the F4 matrix.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:50 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 a801c6e0b60929ee8bf6afe3ca3a8425841582a7
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Fri Feb 1 18:41:20 2013 +0100
Moved reducer 26 to using preblocks (no scalars) for the first version of the F4 matrix.
---
src/mathicgb/F4MatrixBuilder2.cpp | 409 ++++++++++++++++++++++++++------------
src/mathicgb/F4MatrixBuilder2.hpp | 9 +-
src/mathicgb/F4MatrixReducer.cpp | 35 ++--
src/mathicgb/Poly.hpp | 6 +-
src/mathicgb/QuadMatrix.cpp | 14 +-
src/mathicgb/SparseMatrix.cpp | 21 +-
src/mathicgb/SparseMatrix.hpp | 6 +-
7 files changed, 335 insertions(+), 165 deletions(-)
diff --git a/src/mathicgb/F4MatrixBuilder2.cpp b/src/mathicgb/F4MatrixBuilder2.cpp
index 339869c..a47677b 100755
--- a/src/mathicgb/F4MatrixBuilder2.cpp
+++ b/src/mathicgb/F4MatrixBuilder2.cpp
@@ -9,6 +9,121 @@ MATHICGB_DEFINE_LOG_DOMAIN(
"Displays statistics about F4 matrix construction."
);
+class F4MatrixBuilder2::F4PreBlock {
+public:
+ typedef uint32 RowIndex;
+ typedef uint32 ColIndex;
+ typedef coefficient ExternalScalar;
+ typedef SparseMatrix::Scalar Scalar;
+
+ struct Row {
+ const ColIndex* indices;
+ const Scalar* scalars;
+ const ExternalScalar* externalScalars;
+ ColIndex entryCount;
+ };
+
+ RowIndex rowCount() const {return static_cast<RowIndex>(mRows.size());}
+
+ Row row(const RowIndex row) const {
+ MATHICGB_ASSERT(row < mRows.size());
+ const auto& r = mRows[row];
+ Row rr;
+ rr.indices = mIndices.data() + r.indicesBegin;
+ rr.entryCount = r.entryCount;
+ if (r.externalScalars == 0) {
+ rr.scalars = mScalars.data() + r.scalarsBegin;
+ rr.externalScalars = 0;
+ } else {
+ rr.scalars = 0;
+ rr.externalScalars = r.externalScalars;
+ }
+ return rr;
+ }
+
+ ColIndex* makeRowWithTheseScalars(const Poly& scalars) {
+ MATHICGB_ASSERT(rowCount() < std::numeric_limits<RowIndex>::max());
+ MATHICGB_ASSERT
+ (scalars.termCount() < std::numeric_limits<ColIndex>::max());
+
+ InternalRow row;
+ row.indicesBegin = mIndices.size();
+ row.scalarsBegin = std::numeric_limits<decltype(row.scalarsBegin)>::max();
+ row.entryCount = static_cast<ColIndex>(scalars.termCount());
+ row.externalScalars = scalars.coefficientBegin();
+ mRows.push_back(row);
+
+ mIndices.resize(mIndices.size() + row.entryCount);
+ return mIndices.data() + row.indicesBegin;
+ }
+
+ std::pair<ColIndex*, Scalar*> makeRow(ColIndex entryCount) {
+ MATHICGB_ASSERT(rowCount() < std::numeric_limits<RowIndex>::max());
+
+ InternalRow row;
+ row.indicesBegin = mIndices.size();
+ row.scalarsBegin = mScalars.size();
+ row.entryCount = entryCount;
+ row.externalScalars = 0;
+ mRows.push_back(row);
+
+ mIndices.resize(mIndices.size() + entryCount);
+ mScalars.resize(mScalars.size() + entryCount);
+ return std::make_pair(
+ mIndices.data() + row.indicesBegin,
+ mScalars.data() + row.scalarsBegin
+ );
+ }
+
+ void removeLastEntries(const RowIndex row, const ColIndex count) {
+ MATHICGB_ASSERT(row < rowCount());
+ MATHICGB_ASSERT(mRows[row].entryCount >= count);
+ mRows[row].entryCount -= count;
+ if (row != rowCount() - 1)
+ return;
+ mIndices.resize(mIndices.size() - count);
+ if (mRows[row].externalScalars == 0)
+ mScalars.resize(mScalars.size() - count);
+ }
+
+private:
+ struct InternalRow {
+ size_t indicesBegin;
+ size_t scalarsBegin;
+ ColIndex entryCount;
+ const ExternalScalar* externalScalars;
+ };
+
+ std::vector<ColIndex> mIndices;
+ std::vector<Scalar> mScalars;
+ std::vector<InternalRow> mRows;
+};
+
+void toSparseMatrix(const F4MatrixBuilder2::F4PreBlock& block, SparseMatrix& matrix) {
+ typedef F4MatrixBuilder2::F4PreBlock F4PreBlock;
+ const auto rowCount = block.rowCount();
+ for (F4PreBlock::RowIndex r = 0; r < rowCount; ++r) {
+ const auto row = block.row(r);
+ const auto entryCount = row.entryCount;
+ const F4PreBlock::ColIndex* const indices = row.indices;
+ MATHICGB_ASSERT(row.scalars == 0 || row.externalScalars == 0);
+ if (row.scalars != 0) {
+ const F4PreBlock::Scalar* const scalars = row.scalars;
+ for (F4PreBlock::ColIndex col = 0; col < entryCount; ++col)
+ matrix.appendEntry(indices[col], scalars[col]);
+ } else if (row.externalScalars != 0) {
+ const F4PreBlock::ExternalScalar* const scalars = row.externalScalars;
+ for (F4PreBlock::ColIndex col = 0; col < entryCount; ++col) {
+ const auto scalar = static_cast<F4PreBlock::Scalar>(scalars[col]);
+ MATHICGB_ASSERT
+ (static_cast<F4PreBlock::ExternalScalar>(scalar) == scalar);
+ matrix.appendEntry(indices[col], scalar);
+ }
+ }
+ matrix.rowDone();
+ }
+}
+
MATHICGB_NO_INLINE
std::pair<F4MatrixBuilder2::ColIndex, ConstMonomial>
F4MatrixBuilder2::findOrCreateColumn(
@@ -126,13 +241,13 @@ void F4MatrixBuilder2::buildMatrixAndClear(QuadMatrix& quadMatrix) {
// we are calling here can add more pending items.
struct ThreadData {
- SparseMatrix matrix;
+ F4PreBlock block;
monomial tmp1;
monomial tmp2;
};
tbb::enumerable_thread_specific<ThreadData> threadData([&](){
- ThreadData data = {SparseMatrix(mMemoryQuantum)};
+ ThreadData data;
{
tbb::mutex::scoped_lock guard(mCreateColumnLock);
data.tmp1 = ring().allocMonomial();
@@ -155,7 +270,7 @@ void F4MatrixBuilder2::buildMatrixAndClear(QuadMatrix& quadMatrix) {
data.tmp1
);
appendRowSPair
- (&poly, data.tmp1, task.sPairPoly, data.tmp2, data.matrix, feeder);
+ (&poly, data.tmp1, task.sPairPoly, data.tmp2, data.block, feeder);
return;
}
if (task.desiredLead.isNull())
@@ -164,26 +279,10 @@ void F4MatrixBuilder2::buildMatrixAndClear(QuadMatrix& quadMatrix) {
ring().monomialDivide
(task.desiredLead, poly.getLeadMonomial(), data.tmp1);
MATHICGB_ASSERT(ring().hashValid(data.tmp1));
- appendRow(
- data.tmp1,
- task.poly->begin(),
- task.poly->end(),
- data.matrix,
- feeder
- );
+ appendRow(data.tmp1, *task.poly, data.block, feeder);
});
MATHICGB_ASSERT(!threadData.empty()); // as mTodo empty causes early return
- auto matrix = std::move(threadData.begin()->matrix);
- const auto end = threadData.end();
- for (auto it = threadData.begin(); it != end; ++it) {
- if (it != threadData.begin())
- matrix.takeRowsFrom(std::move(it->matrix));
- ring().freeMonomial(it->tmp1);
- ring().freeMonomial(it->tmp2);
- }
- threadData.clear();
-
// Free the monomials from all the tasks
const auto todoEnd = mTodo.end();
for (auto it = mTodo.begin(); it != todoEnd; ++it)
@@ -211,120 +310,136 @@ void F4MatrixBuilder2::buildMatrixAndClear(QuadMatrix& quadMatrix) {
}
// Decide which rows are reducers (top) and which are reducees (bottom).
- const auto rowCount = matrix.rowCount();
- std::vector<RowIndex> reducerRows(mLeftColCount, rowCount);
- std::vector<RowIndex> reduceeRows;
- for (RowIndex row = 0; row < rowCount; ++row) {
- if (matrix.emptyRow(row))
- continue;
-
- MATHICGB_ASSERT(!matrix.emptyRow(row));
- // Determine leading (minimum index) left entry.
- const auto lead = [&] {
- MATHICGB_ASSERT(mTranslate.size() <= std::numeric_limits<ColIndex>::max());
- const auto end = matrix.rowEnd(row);
- for (auto it = matrix.rowBegin(row); it != end; ++it) {
- MATHICGB_ASSERT(it.index() < mTranslate.size());
- if (mTranslate[it.index()].left)
- return mTranslate[it.index()].index;
+ const auto noReducer = std::numeric_limits<RowIndex>::max();
+ F4PreBlock::Row noRow = {};
+ noRow.indices = 0;
+ std::vector<F4PreBlock::Row> reducerRows(mLeftColCount, noRow);
+ std::vector<F4PreBlock::Row> reduceeRows;
+
+ SparseMatrix matrix(mMemoryQuantum);
+ const auto end = threadData.end();
+ for (auto it = threadData.begin(); it != end; ++it) {
+ auto& block = it->block;
+ const auto rowCount = block.rowCount();
+ for (RowIndex r = 0; r < rowCount; ++r) {
+ const auto row = block.row(r);
+ if (row.entryCount == 0)
+ continue;
+
+ // Determine leading (minimum index) left entry.
+ const auto lead = [&] {
+ MATHICGB_ASSERT(mTranslate.size() <=
+ std::numeric_limits<ColIndex>::max());
+ const auto end = row.indices + row.entryCount;
+ for (auto it = row.indices; it != end; ++it) {
+ MATHICGB_ASSERT(*it < mTranslate.size());
+ if (mTranslate[*it].left)
+ return mTranslate[*it].index;
+ }
+ return mLeftColCount; // No left entries at all.
+ }();
+ if (!mTranslate[*row.indices].left) {
+ reduceeRows.push_back(row); // no left entries
+ continue; // todo: for now
}
- return mLeftColCount; // no left entries at all
- }();
- if (!mTranslate[matrix.leadCol(row)].left) {
- reduceeRows.push_back(row); // no left entries
- continue;
+ // Decide if this should be a reducer or reducee row.
+ if (lead == mLeftColCount) {
+ reduceeRows.push_back(row); // no left entries
+ continue;
+ }
+ auto& reducer = reducerRows[lead];
+ if (reducer.entryCount == 0)
+ reducer = row; // row is first reducer with this lead
+ else if (reducer.entryCount > row.entryCount) {
+ reduceeRows.push_back(reducer); // row sparser (=better) reducer
+ reducer = row;
+ } else
+ reduceeRows.push_back(row);
}
- // decide if this is a reducer or reducee row
- if (lead == mLeftColCount) {
- reduceeRows.push_back(row); // no left entries
- continue;
- }
- auto& reducer = reducerRows[lead];
- if (reducer == rowCount)
- reducer = row; // row is first reducer with this lead
- else if (matrix.entryCountInRow(reducer) > matrix.entryCountInRow(row)) {
- reduceeRows.push_back(reducer); // row sparser (=better) reducer
- reducer = row;
- } else
- reduceeRows.push_back(row);
+ ring().freeMonomial(it->tmp1);
+ ring().freeMonomial(it->tmp2);
}
MATHICGB_ASSERT(reducerRows.size() == mLeftColCount);
#ifdef MATHICGB_DEBUG
for (size_t i = 0; i < reducerRows.size(); ++i) {
const auto row = reducerRows[i];
- MATHICGB_ASSERT(row < matrix.rowCount());
- MATHICGB_ASSERT(!matrix.emptyRow(row));
- MATHICGB_ASSERT(mTranslate[matrix.leadCol(row)].left);
- MATHICGB_ASSERT(mTranslate[matrix.leadCol(row)].index == i);
+ MATHICGB_ASSERT(row.entryCount > 0);
+ MATHICGB_ASSERT(mTranslate[*row.indices].left);
+ MATHICGB_ASSERT(mTranslate[*row.indices].index == i);
}
#endif
quadMatrix.ring = &ring();
auto splitLeftRight = [this](
- SparseMatrix& from,
- const std::vector<RowIndex>& fromRows,
+ const std::vector<F4PreBlock::Row>& from,
const bool makeLeftUnitary,
SparseMatrix& left,
SparseMatrix& right
) {
left.clear();
right.clear();
- const auto fromRowsEnd = fromRows.end();
- for (
- auto fromRowsIt = fromRows.begin();
- fromRowsIt != fromRowsEnd;
- ++fromRowsIt
- ) {
- const auto row = *fromRowsIt;
- const auto fromEnd = from.rowEnd(row);
- auto fromIt = from.rowBegin(row);
- MATHICGB_ASSERT(!from.emptyRow(row));
- if (
- makeLeftUnitary &&
- (!mTranslate[fromIt.index()].left || fromIt.scalar() != 1)
- ) {
- auto firstLeft = fromIt;
- while (!mTranslate[firstLeft.index()].left) {
- // We cannot make the left part unitary if the left part is a zero
- // row, so makeLeftUnitary implies no zero left rows.
- MATHICGB_ASSERT(firstLeft != fromEnd);
- ++firstLeft;
+ const auto fromEnd = from.end();
+ for (auto fromIt = from.begin(); fromIt != fromEnd; ++fromIt) {
+ const auto row = *fromIt;
+ MATHICGB_ASSERT(row.entryCount != 0);
+ MATHICGB_ASSERT(row.scalars == 0 || row.externalScalars == 0);
+
+ if (row.externalScalars != 0) {
+ auto indices = row.indices;
+ auto indicesEnd = row.indices + row.entryCount;
+ auto scalars = row.externalScalars;
+ for (; indices != indicesEnd; ++indices, ++scalars) {
+ const auto scalar = static_cast<Scalar>(*scalars);
+ const auto index = *indices;
+ const auto translated = mTranslate[index];
+ if (translated.left)
+ left.appendEntry(translated.index, scalar);
+ else
+ right.appendEntry(translated.index, scalar);
}
- if (firstLeft.scalar() != 1) {
- const auto modulus = static_cast<Scalar>(ring().charac());
- const auto inverse = modularInverse(firstLeft.scalar(), modulus);
- for (auto it = fromIt; it != fromEnd; ++it)
- it.setScalar(modularProduct(it.scalar(), inverse, modulus));
- MATHICGB_ASSERT(firstLeft.scalar() == 1);
+ } else {
+ auto indices = row.indices;
+ auto indicesEnd = row.indices + row.entryCount;
+ auto scalars = row.scalars;
+ for (; indices != indicesEnd; ++indices, ++scalars) {
+ const auto index = *indices;
+ const auto translated = mTranslate[index];
+ if (translated.left)
+ left.appendEntry(translated.index, *scalars);
+ else
+ right.appendEntry(translated.index, *scalars);
}
}
- for (; fromIt != fromEnd; ++fromIt) {
- MATHICGB_ASSERT(fromIt.index() < mTranslate.size());
- const auto translated = mTranslate[fromIt.index()];
- if (translated.left)
- left.appendEntry(translated.index, fromIt.scalar());
- else
- right.appendEntry(translated.index, fromIt.scalar());
- }
+ const auto rowIndex = left.rowCount();
+ MATHICGB_ASSERT(rowIndex == right.rowCount());
left.rowDone();
right.rowDone();
+
+ if (
+ makeLeftUnitary &&
+ !left.emptyRow(rowIndex) &&
+ left.rowBegin(rowIndex).scalar() != 1
+ ) {
+ const auto modulus = static_cast<Scalar>(ring().charac());
+ const auto inverse =
+ modularInverse(left.rowBegin(rowIndex).scalar(), modulus);
+ left.multiplyRow(rowIndex, inverse, modulus);
+ right.multiplyRow(rowIndex, inverse, modulus);
+ MATHICGB_ASSERT(left.rowBegin(rowIndex).scalar() == 1);
+ }
+
MATHICGB_ASSERT(left.rowCount() == right.rowCount());
MATHICGB_ASSERT(!makeLeftUnitary || !left.emptyRow(left.rowCount() - 1));
MATHICGB_ASSERT
(!makeLeftUnitary || left.rowBegin(left.rowCount() - 1).scalar() == 1);
}
};
+ splitLeftRight(reducerRows, true, quadMatrix.topLeft, quadMatrix.topRight);
splitLeftRight
- (matrix, reducerRows, true, quadMatrix.topLeft, quadMatrix.topRight);
- splitLeftRight(
- matrix,
- reduceeRows,
- false,
- quadMatrix.bottomLeft,
- quadMatrix.bottomRight
- );
+ (reduceeRows, false, quadMatrix.bottomLeft, quadMatrix.bottomRight);
+ threadData.clear();
#ifdef MATHICGB_DEBUG
for (size_t side = 0; side < 2; ++side) {
@@ -334,11 +449,12 @@ void F4MatrixBuilder2::buildMatrixAndClear(QuadMatrix& quadMatrix) {
MATHICGB_ASSERT(!it->isNull());
}
}
- for (size_t row = 0; row < quadMatrix.topLeft.rowCount(); ++row) {
+ for (RowIndex row = 0; row < quadMatrix.topLeft.rowCount(); ++row) {
MATHICGB_ASSERT(quadMatrix.topLeft.leadCol(row) == row);
}
#endif
+ // todo: do this together with left/right split
quadMatrix.sortColumnsLeftRightParallel();
#ifdef MATHICGB_DEBUG
MATHICGB_ASSERT(quadMatrix.debugAssertValid());
@@ -401,21 +517,30 @@ F4MatrixBuilder2::createColumn(
void F4MatrixBuilder2::appendRow(
const const_monomial multiple,
- const Poly::const_iterator begin,
- const Poly::const_iterator end,
- SparseMatrix& matrix,
+ const Poly& poly,
+ F4PreBlock& block,
TaskFeeder& feeder
) {
MATHICGB_ASSERT(!multiple.isNull());
+ const auto begin = poly.begin();
+ const auto end = poly.end();
+ const auto count = std::distance(begin, end);
+ MATHICGB_ASSERT(count < std::numeric_limits<ColIndex>::max());
+ auto indices = block.makeRowWithTheseScalars(poly);
+
auto it = begin;
- if ((std::distance(it, end) % 2) == 1) {
+ if ((count % 2) == 1) {
ColReader reader(mMap);
const auto col = findOrCreateColumn
(it.getMonomial(), multiple, reader, feeder);
MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
MATHICGB_ASSERT(it.getCoefficient());
- matrix.appendEntry(col.first, static_cast<Scalar>(it.getCoefficient()));
+ //matrix.appendEntry(col.first, static_cast<Scalar>(it.getCoefficient()));
+ *indices = col.first;
+ ++indices;
+ //*row.first++ = col.first;
+ //*row.second++ = static_cast<Scalar>(it.getCoefficient());
++it;
}
updateReader:
@@ -440,11 +565,22 @@ updateReader:
goto updateReader;
}
- matrix.appendEntry(*colPair.first, scalar1);
- matrix.appendEntry(*colPair.second, scalar2);
+ //matrix.appendEntry(*colPair.first, scalar1);
+ //matrix.appendEntry(*colPair.second, scalar2);
+
+ *indices = *colPair.first;
+ ++indices;
+ *indices = *colPair.second;
+ ++indices;
+
+ //*row.first++ = *colPair.first;
+ //*row.second++ = scalar1;
+ //*row.first++ = *colPair.second;
+ //*row.second++ = scalar2;
+
it = ++it2;
}
- matrix.rowDone();
+ //matrix.rowDone();
}
void F4MatrixBuilder2::appendRowSPair(
@@ -452,7 +588,7 @@ void F4MatrixBuilder2::appendRowSPair(
monomial multiply,
const Poly* sPairPoly,
monomial sPairMultiply,
- SparseMatrix& matrix,
+ F4PreBlock& block,
TaskFeeder& feeder
) {
MATHICGB_ASSERT(!poly->isZero());
@@ -473,6 +609,14 @@ void F4MatrixBuilder2::appendRowSPair(
++itA;
++itB;
+ // @todo: handle overflow of termCount addition here
+ MATHICGB_ASSERT(poly->termCount() + sPairPoly->termCount() - 2 <=
+ std::numeric_limits<ColIndex>::max());
+ const auto maxCols =
+ static_cast<ColIndex>(poly->termCount() + sPairPoly->termCount() - 2);
+ auto row = block.makeRow(maxCols);
+ const auto indicesBegin = row.first;
+
const ColReader colMap(mMap);
const const_monomial mulA = multiply;
@@ -497,24 +641,31 @@ void F4MatrixBuilder2::appendRowSPair(
++itB;
}
MATHICGB_ASSERT(coeff < std::numeric_limits<Scalar>::max());
- if (coeff != 0)
- matrix.appendEntry(col, static_cast<Scalar>(coeff));
+ if (coeff != 0) {
+ //matrix.appendEntry(col, static_cast<Scalar>(coeff));
+ *row.first++ = col;
+ *row.second++ = static_cast<Scalar>(coeff);
+ }
}
- // these calls also end the row.
- if (itA != endA)
- appendRow(mulA, itA, endA, matrix, feeder);
- else {
- const auto toNegateCount = std::distance(itB, endB);
- appendRow(mulB, itB, endB, matrix, feeder);
- const auto row = matrix.rowCount() - 1;
- const auto end = matrix.rowEnd(row);
- auto it = matrix.rowBegin(row);
- it += matrix.entryCountInRow(row) - toNegateCount;
- for (; it != end; ++it) {
- const auto negative = ring().coefficientNegate(it.scalar());
- MATHICGB_ASSERT(negative < std::numeric_limits<Scalar>::max());
- it.setScalar(static_cast<Scalar>(negative));
- }
+ for (; itA != endA; ++itA) {
+ const auto colA = findOrCreateColumn
+ (itA.getMonomial(), mulA, colMap, feeder);
+ //matrix.appendEntry(colA.first, static_cast<Scalar>(itA.getCoefficient()));
+ *row.first++ = colA.first;
+ *row.second++ = static_cast<Scalar>(itA.getCoefficient());
}
+
+ for (; itB != endB; ++itB) {
+ const auto colB = findOrCreateColumn
+ (itB.getMonomial(), mulB, colMap, feeder);
+ const auto negative = ring().coefficientNegate(itB.getCoefficient());
+ //matrix.appendEntry(colB.first, static_cast<Scalar>(negative));
+ *row.first++ = colB.first;
+ *row.second++ = static_cast<Scalar>(negative);
+ }
+
+ const auto toRemove =
+ maxCols - static_cast<ColIndex>(row.first - indicesBegin);
+ block.removeLastEntries(block.rowCount() - 1, toRemove);
}
diff --git a/src/mathicgb/F4MatrixBuilder2.hpp b/src/mathicgb/F4MatrixBuilder2.hpp
index 6274943..04aef01 100755
--- a/src/mathicgb/F4MatrixBuilder2.hpp
+++ b/src/mathicgb/F4MatrixBuilder2.hpp
@@ -65,6 +65,8 @@ public:
const PolyRing& ring() const {return mBasis.ring();}
+ class F4PreBlock;
+
private:
typedef const Map::Reader ColReader;
typedef std::vector<monomial> Monomials;
@@ -93,9 +95,8 @@ private:
void appendRow(
const_monomial multiple,
- const Poly::const_iterator begin,
- const Poly::const_iterator end,
- SparseMatrix& matrix,
+ const Poly& poly,
+ F4PreBlock& block,
TaskFeeder& feeder
);
void appendRowSPair(
@@ -103,7 +104,7 @@ private:
monomial multiply,
const Poly* sPairPoly,
monomial sPairMultiply,
- SparseMatrix& matrix,
+ F4PreBlock& block,
TaskFeeder& feeder
);
diff --git a/src/mathicgb/F4MatrixReducer.cpp b/src/mathicgb/F4MatrixReducer.cpp
index 730d1aa..bab9c49 100644
--- a/src/mathicgb/F4MatrixReducer.cpp
+++ b/src/mathicgb/F4MatrixReducer.cpp
@@ -153,7 +153,7 @@ namespace {
}
void rowReduceByUnitary(
- const size_t pivotRow,
+ const SparseMatrix::RowIndex pivotRow,
const SparseMatrix& matrix,
const SparseMatrix::Scalar modulus
) {
@@ -229,11 +229,11 @@ namespace {
std::vector<SparseMatrix::RowIndex> rowOrder(rowCount);
tbb::mutex lock;
- tbb::parallel_for(tbb::blocked_range<size_t>(0, rowCount),
- [&](const tbb::blocked_range<size_t>& range)
+ tbb::parallel_for(tbb::blocked_range<SparseMatrix::RowIndex>(0, rowCount),
+ [&](const tbb::blocked_range<SparseMatrix::RowIndex>& range)
{for (auto it = range.begin(); it != range.end(); ++it)
{
- const size_t row = it;
+ const auto row = it;
auto& denseRow = denseRowPerThread.local();
denseRow.clear(leftColCount);
@@ -272,12 +272,12 @@ namespace {
rowOrder[tmp.rowCount() - 1] = row;
}});
- tbb::parallel_for(tbb::blocked_range<size_t>(0, rowCount),
- [&](const tbb::blocked_range<size_t>& range)
+ tbb::parallel_for(tbb::blocked_range<SparseMatrix::RowIndex>(0, rowCount),
+ [&](const tbb::blocked_range<SparseMatrix::RowIndex>& range)
{for (auto iter = range.begin(); iter != range.end(); ++iter)
{
- const size_t i = iter;
- const size_t row = rowOrder[i];
+ const auto i = iter;
+ const auto row = rowOrder[i];
auto& denseRow = denseRowPerThread.local();
denseRow.clear(rightColCount);
@@ -315,11 +315,11 @@ namespace {
// convert to dense representation
std::vector<DenseRow> dense(rowCount);
- tbb::parallel_for(tbb::blocked_range<size_t>(0, rowCount),
- [&](const tbb::blocked_range<size_t>& range)
+ tbb::parallel_for(tbb::blocked_range<SparseMatrix::RowIndex>(0, rowCount),
+ [&](const tbb::blocked_range<SparseMatrix::RowIndex>& range)
{for (auto it = range.begin(); it != range.end(); ++it)
{
- const size_t row = it;
+ const auto row = it;
if (toReduce.emptyRow(row))
continue;
dense[row].clear(colCount);
@@ -350,19 +350,19 @@ namespace {
//std::cout << "reducing " << reduced.rowCount() << " out of " << toReduce.rowCount() << std::endl;
tbb::mutex lock;
- tbb::parallel_for(tbb::blocked_range<size_t>(0, rowCount),
- [&](const tbb::blocked_range<size_t>& range)
+ tbb::parallel_for(tbb::blocked_range<SparseMatrix::RowIndex>(0, rowCount),
+ [&](const tbb::blocked_range<SparseMatrix::RowIndex>& range)
{for (auto it = range.begin(); it != range.end(); ++it)
{
- const size_t row = it;
+ const auto row = it;
MATHICGB_ASSERT(leadCols[row] <= colCount);
DenseRow& denseRow = dense[row];
if (denseRow.empty())
continue;
// reduce by each row of reduced.
- for (size_t reducerRow = 0; reducerRow < reducerCount; ++reducerRow) {
- size_t const col = reduced.rowBegin(reducerRow).index();
+ for (SparseMatrix::RowIndex reducerRow = 0; reducerRow < reducerCount; ++reducerRow) {
+ const auto col = reduced.rowBegin(reducerRow).index();
if (denseRow[col] == 0 || (isPivotRow[row] && col == leadCols[row]))
continue;
denseRow.rowReduceByUnitary(reducerRow, reduced, modulus);
@@ -556,7 +556,8 @@ void rowReducedEchelonMatrix(
const SparseMatrix::Scalar modulus
) {
assert(matrix.empty() || matrix[0].size() == colCount);
- const SparseMatrix::RowIndex rowCount=matrix.size();
+ const SparseMatrix::RowIndex rowCount =
+ static_cast<SparseMatrix::RowIndex>(matrix.size());
// pivotRowOfCol[i] is the pivot in column i or rowCount
// if we have not identified such a pivot so far.
std::vector<SparseMatrix::RowIndex> pivotRowOfCol(colCount, rowCount);
diff --git a/src/mathicgb/Poly.hpp b/src/mathicgb/Poly.hpp
index 61a2c5b..657c148 100755
--- a/src/mathicgb/Poly.hpp
+++ b/src/mathicgb/Poly.hpp
@@ -105,6 +105,8 @@ public:
iterator begin() { return iterator(*this); }
iterator end() { return iterator(*this,1); }
+ const coefficient* coefficientBegin() const {return coeffs.data();}
+
static Poly * add(const PolyRing *R,
iterator first1,
@@ -124,7 +126,9 @@ public:
const_coefficient getLeadCoefficient() const { return coeffs[0]; }
exponent getLeadComponent() const { return R->monomialGetComponent(&(monoms[0])); }
bool isZero() const { return coeffs.empty(); }
- size_t nTerms() const { return coeffs.size(); }
+
+ size_t nTerms() const { return coeffs.size(); } /// @todo: deprecated
+ size_t termCount() const {return coeffs.size();}
size_t getMemoryUse() const;
diff --git a/src/mathicgb/QuadMatrix.cpp b/src/mathicgb/QuadMatrix.cpp
index 29ce67a..a39395f 100755
--- a/src/mathicgb/QuadMatrix.cpp
+++ b/src/mathicgb/QuadMatrix.cpp
@@ -176,7 +176,7 @@ QuadMatrix QuadMatrix::toCanonical() const {
class RowComparer {
public:
RowComparer(const SparseMatrix& matrix): mMatrix(matrix) {}
- bool operator()(size_t a, size_t b) const {
+ bool operator()(SparseMatrix::RowIndex a, SparseMatrix::RowIndex b) const {
// if you need this to work for empty rows or identical leading columns
// then update this code.
MATHICGB_ASSERT(!mMatrix.emptyRow(a));
@@ -212,8 +212,8 @@ QuadMatrix QuadMatrix::toCanonical() const {
// todo: eliminate left/right code duplication here
QuadMatrix matrix;
{ // left side
- std::vector<size_t> rows;
- for (size_t row = 0; row < topLeft.rowCount(); ++row)
+ std::vector<SparseMatrix::RowIndex> rows;
+ for (SparseMatrix::RowIndex row = 0; row < topLeft.rowCount(); ++row)
rows.push_back(row);
{
RowComparer comparer(topLeft);
@@ -222,14 +222,14 @@ QuadMatrix QuadMatrix::toCanonical() const {
matrix.topLeft.clear();
matrix.topRight.clear();
- for (size_t i = 0; i < rows.size(); ++i) {
+ for (SparseMatrix::RowIndex i = 0; i < rows.size(); ++i) {
matrix.topLeft.appendRow(topLeft, rows[i]);
matrix.topRight.appendRow(topRight, rows[i]);
}
}
{ // right side
- std::vector<size_t> rows;
- for (size_t row = 0; row < bottomLeft.rowCount(); ++row)
+ std::vector<SparseMatrix::RowIndex> rows;
+ for (SparseMatrix::RowIndex row = 0; row < bottomLeft.rowCount(); ++row)
rows.push_back(row);
{
RowComparer comparer(bottomLeft);
@@ -238,7 +238,7 @@ QuadMatrix QuadMatrix::toCanonical() const {
matrix.bottomLeft.clear();
matrix.bottomRight.clear();
- for (size_t i = 0; i < rows.size(); ++i) {
+ for (SparseMatrix::RowIndex i = 0; i < rows.size(); ++i) {
matrix.bottomLeft.appendRow(bottomLeft, rows[i]);
matrix.bottomRight.appendRow(bottomRight, rows[i]);
}
diff --git a/src/mathicgb/SparseMatrix.cpp b/src/mathicgb/SparseMatrix.cpp
index 7e5e4e6..9957500 100755
--- a/src/mathicgb/SparseMatrix.cpp
+++ b/src/mathicgb/SparseMatrix.cpp
@@ -48,10 +48,10 @@ void SparseMatrix::sortRowsByIncreasingPivots() {
const auto rowCount = this->rowCount();
std::vector<RowIndex> rows(rowCount);
- for (size_t row = 0; row < rowCount; ++row)
+ for (RowIndex row = 0; row < rowCount; ++row)
rows[row] = row;
- const auto lexLess = [&](const size_t a, const size_t b) -> bool {
+ const auto lexLess = [&](const RowIndex a, const RowIndex b) -> bool {
auto aIt = rowBegin(a);
auto bIt = rowBegin(b);
const auto aEnd = rowEnd(a);
@@ -85,6 +85,17 @@ void SparseMatrix::applyColumnMap(const std::vector<ColIndex>& colMap) {
}
}
+void SparseMatrix::multiplyRow(
+ const RowIndex row,
+ const Scalar multiplier,
+ const Scalar modulus
+) {
+ MATHICGB_ASSERT(row < rowCount());
+ const auto end = rowEnd(row);
+ for (auto it = rowBegin(row); it != end; ++it)
+ it.setScalar(modularProduct(it.scalar(), multiplier, modulus));
+}
+
void SparseMatrix::print(std::ostream& out) const {
if (rowCount() == 0)
out << "matrix with no rows\n";
@@ -175,7 +186,7 @@ SparseMatrix& SparseMatrix::operator=(const SparseMatrix& matrix) {
// A version that works on each block would be faster, but this is not
// used anywhere time-critical right now. Improve this if it turns
// up in profiling at some point.
- for (size_t row = 0; row < matrix.rowCount(); ++row)
+ for (RowIndex row = 0; row < matrix.rowCount(); ++row)
appendRow(matrix, row);
return *this;
}
@@ -191,7 +202,7 @@ bool SparseMatrix::operator==(const SparseMatrix& matrix) const {
const auto count = rowCount();
if (count != matrix.rowCount())
return false;
- for (size_t row = 0; row < count; ++row) {
+ for (RowIndex row = 0; row < count; ++row) {
if (entryCountInRow(row) != matrix.entryCountInRow(row))
return false;
const auto end = rowEnd(row);
@@ -208,7 +219,7 @@ SparseMatrix::ColIndex SparseMatrix::computeColCount() const {
// Obviously this can be done faster, but there has not been a need for that
// so far.
ColIndex colCount = 0;
- for (size_t row = 0; row < rowCount(); ++row) {
+ for (RowIndex row = 0; row < rowCount(); ++row) {
const auto end = rowEnd(row);
for (auto it = rowBegin(row); it != end; ++it)
colCount = std::max(colCount, it.index() + 1);
diff --git a/src/mathicgb/SparseMatrix.hpp b/src/mathicgb/SparseMatrix.hpp
index 09ed6d0..bbeef08 100755
--- a/src/mathicgb/SparseMatrix.hpp
+++ b/src/mathicgb/SparseMatrix.hpp
@@ -42,7 +42,7 @@ support a wider range of types of matrices in future.
*/
class SparseMatrix {
public:
- typedef size_t RowIndex;
+ typedef uint32 RowIndex;
typedef uint32 ColIndex;
typedef uint16 Scalar;
class ConstRowIterator;
@@ -88,7 +88,7 @@ public:
/// the memory out of matrix.
void takeRowsFrom(SparseMatrix&& matrix);
- RowIndex rowCount() const {return mRows.size();}
+ RowIndex rowCount() const {return static_cast<RowIndex>(mRows.size());}
ColIndex computeColCount() const;
size_t memoryQuantum() const {return mMemoryQuantum;}
@@ -231,6 +231,8 @@ public:
/// Replaces all column indices i with colMap[i].
void applyColumnMap(const std::vector<ColIndex>& colMap);
+ void multiplyRow(RowIndex row, Scalar multiplier, Scalar modulus);
+
/// Let poly be the dot product of colMonomials and the given row.
void rowToPolynomial(
RowIndex row,
--
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