[mathicgb] 173/393: Improved code in F4MatrixProjection and added ScopeExit functionality.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:55 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 c0732a94c75a43dbd34d3f5d988e590992e91284
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Mon Feb 18 16:39:36 2013 +0100
Improved code in F4MatrixProjection and added ScopeExit functionality.
---
Makefile.am | 2 +-
build/vs12/mathicgb-lib/mathicgb-lib.vcxproj | 1 +
.../vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters | 3 +
src/mathicgb/F4MatrixBuilder2.cpp | 26 +-
src/mathicgb/F4MatrixProjection.cpp | 706 ++++++++++-----------
src/mathicgb/F4MatrixProjection.hpp | 110 +++-
src/mathicgb/MonomialMap.hpp | 2 +-
src/mathicgb/ScopeExit.hpp | 98 +++
8 files changed, 534 insertions(+), 414 deletions(-)
diff --git a/Makefile.am b/Makefile.am
index 9d903bb..9a2af76 100755
--- a/Makefile.am
+++ b/Makefile.am
@@ -68,7 +68,7 @@ libmathicgb_la_SOURCES = src/mathicgb/BjarkeGeobucket2.cpp \
src/mathicgb/F4MatrixBuilder2.hpp src/mathicgb/F4MatrixBuilder2.cpp \
src/mathicgb/LogDomainSet.cpp src/mathicgb/ProtoMatrix.hpp \
src/mathicgb/ProtoMatrix.cpp src/mathicgb/F4MatrixProject.hpp \
- src/mathicgb/F4MatrixProjection.cpp
+ src/mathicgb/F4MatrixProjection.cpp src/mathicgb/ScopeExit.hpp
# When making a distribution file, Automake knows to include all files
# that are necessary to build the project. EXTRA_DIST specifies files
diff --git a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj
index 754053e..373b9b7 100755
--- a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj
+++ b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj
@@ -531,6 +531,7 @@
<ClInclude Include="..\..\..\src\mathicgb\ReducerNoDedup.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\ReducerPack.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\ReducerPackDedup.hpp" />
+ <ClInclude Include="..\..\..\src\mathicgb\ScopeExit.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\SignatureGB.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\SigSPairQueue.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\SigSPairs.hpp" />
diff --git a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters
index fde5175..3390b61 100755
--- a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters
+++ b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters
@@ -311,5 +311,8 @@
<ClInclude Include="..\..\..\src\mathicgb\F4MatrixProjection.hpp">
<Filter>Header Files</Filter>
</ClInclude>
+ <ClInclude Include="..\..\..\src\mathicgb\ScopeExit.hpp">
+ <Filter>Header Files</Filter>
+ </ClInclude>
</ItemGroup>
</Project>
\ No newline at end of file
diff --git a/src/mathicgb/F4MatrixBuilder2.cpp b/src/mathicgb/F4MatrixBuilder2.cpp
index f21a526..0551718 100755
--- a/src/mathicgb/F4MatrixBuilder2.cpp
+++ b/src/mathicgb/F4MatrixBuilder2.cpp
@@ -173,18 +173,32 @@ void F4MatrixBuilder2::buildMatrixAndClear(QuadMatrix& quadMatrix) {
ring().freeMonomial(it->desiredLead);
mTodo.clear();
- // Collect pre-blocks from each thread
- std::vector<F4ProtoMatrix*> blocks;
- blocks.reserve(threadData.size());
+ // Move the proto-matrices across all threads into the projection.
+ F4MatrixProjection projection
+ (ring(), static_cast<ColIndex>(mMap.entryCount()));
const auto end = threadData.end();
for (auto it = threadData.begin(); it != end; ++it) {
- blocks.emplace_back(&it->block);
+ projection.addProtoMatrix(std::move(it->block));
ring().freeMonomial(it->tmp1);
ring().freeMonomial(it->tmp2);
- continue;
}
- quadMatrix = F4MatrixProjection().project(std::move(mIsColumnToLeft), std::move(mMap), std::move(blocks));
+ // Sort columns by monomial and tell the projection of the resulting order
+ MonomialMap<ColIndex>::Reader reader(mMap);
+ typedef std::pair<ColIndex, ConstMonomial> IndexMono;
+ std::vector<IndexMono> columns(reader.begin(), reader.end());
+ const auto cmp = [&](const IndexMono& a, const IndexMono b) {
+ return ring().monomialLT(b.second, a.second);
+ };
+ tbb::parallel_sort(columns.begin(), columns.end(), cmp);
+
+ const auto colEnd = columns.end();
+ for (auto it = columns.begin(); it != colEnd; ++it) {
+ const auto p = *it;
+ projection.addColumn(p.first, p.second, mIsColumnToLeft[p.first]);
+ }
+
+ quadMatrix = projection.makeProjectionAndClear();
threadData.clear();
#ifdef MATHICGB_DEBUG
diff --git a/src/mathicgb/F4MatrixProjection.cpp b/src/mathicgb/F4MatrixProjection.cpp
index d5f6b89..ffed5c4 100644
--- a/src/mathicgb/F4MatrixProjection.cpp
+++ b/src/mathicgb/F4MatrixProjection.cpp
@@ -1,417 +1,361 @@
-#include "stdinc.h"
-#include "F4MatrixProjection.hpp"
-
-namespace {
- class LeftRightProjection {
- public:
- typedef SparseMatrix::ColIndex ColIndex;
-
- LeftRightProjection(
- const std::vector<char>& isColToLeft,
- const MonomialMap<ColIndex>& map
- ) {
- const auto& ring = map.ring();
-
- // Sort columns by monomial while keeping track of original index.
- MonomialMap<ColIndex>::Reader reader(map);
- typedef std::pair<ColIndex, ConstMonomial> IndexMono;
- std::vector<IndexMono> columns(reader.begin(), reader.end());
- const auto cmp = [&ring](const IndexMono& a, const IndexMono b) {
- return ring.monomialLT(b.second, a.second);
- };
- tbb::parallel_sort(columns.begin(), columns.end(), cmp);
-
- // Copy monomials and construct projection mapping.
- MATHICGB_ASSERT
- (isColToLeft.size() <= std::numeric_limits<ColIndex>::max());
- ColIndex colCount = static_cast<ColIndex>(isColToLeft.size());
- mProject.resize(isColToLeft.size());
- for (size_t i = 0; i < colCount; ++i) {
- const auto indexMono = columns[i];
- monomial mono = ring.allocMonomial();
- ring.monomialCopy(indexMono.second, mono);
-
- auto& projected = mProject[indexMono.first];
- projected.left = isColToLeft[indexMono.first];
- if (projected.left) {
- projected.index = static_cast<ColIndex>(mLeftMonomials.size());
- mLeftMonomials.push_back(mono);
- } else {
- projected.index = static_cast<ColIndex>(mRightMonomials.size());
- mRightMonomials.push_back(mono);
- }
- }
- MATHICGB_ASSERT
- (mLeftMonomials.size() + mRightMonomials.size() == isColToLeft.size());
- }
-
- struct Projected {
- ColIndex index;
- bool left;
- };
+#include "stdinc.h"
+#include "F4MatrixProjection.hpp"
+
+#include "ScopeExit.hpp"
+
+F4MatrixProjection::F4MatrixProjection(
+ const PolyRing& ring,
+ ColIndex colCount
+):
+ mRing(ring),
+ mColProjectTo(colCount)
+{}
+
+void F4MatrixProjection::addColumn(
+ const ColIndex projectFrom,
+ const const_monomial mono,
+ const bool isLeft
+) {
+ MATHICGB_ASSERT(projectFrom < mColProjectTo.size());
+ MATHICGB_ASSERT
+ (mLeftMonomials.size() + mRightMonomials.size() < mColProjectTo.size());
+
+ auto monoCopy = mRing.allocMonomial();
+ MATHICGB_SCOPE_EXIT(monoGuard) {mRing.freeMonomial(monoCopy);};
+ mRing.monomialCopy(mono, monoCopy);
+
+ auto& projected = mColProjectTo[projectFrom];
+ if (isLeft) {
+ projected.isLeft = true;
+ projected.index = static_cast<ColIndex>(mLeftMonomials.size());
+ mLeftMonomials.push_back(monoCopy);
+ } else {
+ projected.isLeft = false;
+ projected.index = static_cast<ColIndex>(mRightMonomials.size());
+ mRightMonomials.push_back(monoCopy);
+ }
- Projected project(const ColIndex index) const {
- MATHICGB_ASSERT(index < mProject.size());
- return mProject[index];
- }
+ monoGuard.dismiss();
+}
- void project(
- const std::vector<F4ProtoMatrix*>& preBlocks,
- SparseMatrix& left,
- SparseMatrix& right,
- const PolyRing& ring
- ) const {
- left.clear();
- right.clear();
- const auto modulus = static_cast<SparseMatrix::Scalar>(ring.charac());
-
- const auto end = preBlocks.end();
- for (auto it = preBlocks.begin(); it != end; ++it) {
- auto& block = **it;
- const auto rowCount = block.rowCount();
- for (SparseMatrix::RowIndex r = 0; r < rowCount; ++r) {
- const auto row = block.row(r);
- if (row.entryCount == 0)
- continue;
- 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<SparseMatrix::Scalar>(*scalars);
- const auto index = *indices;
- const auto translated = project(index);
- if (translated.left)
- left.appendEntry(translated.index, scalar);
- else
- right.appendEntry(translated.index, scalar);
- }
- } 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 = project(index);
- if (translated.left)
- left.appendEntry(translated.index, *scalars);
- else
- right.appendEntry(translated.index, *scalars);
- }
- }
- MATHICGB_ASSERT(left.rowCount() == right.rowCount());
- left.rowDone();
- right.rowDone();
+void F4MatrixProjection::setupRowProjection() {
+ mTopRowProjectFrom.resize(mLeftMonomials.size());
+
+ const auto noReducer = std::numeric_limits<RowIndex>::max();
+ F4ProtoMatrix::Row noRow = {};
+ noRow.indices = 0;
+ const auto noCol = std::numeric_limits<SparseMatrix::ColIndex>::max();
+ const auto modulus = static_cast<SparseMatrix::Scalar>(ring().charac());
+
+ const auto end = mMatrices.end();
+ for (auto it = mMatrices.begin(); it != end; ++it) {
+ auto& matrix = **it;
+ const auto rowCount = matrix.rowCount();
+ for (RowIndex r = 0; r < rowCount; ++r) {
+ const auto row = matrix.row(r);
+ if (row.entryCount == 0)
+ continue;
+
+ // Determine leading (minimum index) left entry.
+ const auto lead = [&] {
+ for (SparseMatrix::ColIndex col = 0; col < row.entryCount; ++col) {
+ auto const translated = mColProjectTo[row.indices[col]];
+ if (translated.isLeft)
+ return std::make_pair(col, translated.index);
}
+ return std::make_pair(noCol, noCol);
+ }();
+ // Decide if this should be a reducer or reducee row.
+ if (lead.second == noCol) {
+ const RowProjectFrom p = {1, row};
+ mBottomRowProjectFrom.push_back(p); // no left entries
+ continue;
}
- }
-
- void project(
- const std::vector<std::pair<SparseMatrix::Scalar, F4ProtoMatrix::Row>>& from,
- SparseMatrix& left,
- SparseMatrix& right,
- const PolyRing& ring
- ) const {
- left.clear();
- right.clear();
- const auto fromEnd = from.end();
- for (auto fromIt = from.begin(); fromIt != fromEnd; ++fromIt) {
- const auto row = fromIt->second;
- MATHICGB_ASSERT(row.entryCount != 0);
- MATHICGB_ASSERT(row.scalars == 0 || row.externalScalars == 0);
- const auto modulus = static_cast<SparseMatrix::Scalar>(ring.charac());
-
- 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<SparseMatrix::Scalar>(*scalars);
- const auto index = *indices;
- const auto translated = project(index);
- if (translated.left)
- left.appendEntry(translated.index, scalar);
- else
- right.appendEntry(translated.index, scalar);
- }
- } 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 = project(index);
- if (translated.left)
- left.appendEntry(translated.index, *scalars);
- else
- right.appendEntry(translated.index, *scalars);
- }
+ MATHICGB_ASSERT(row.scalars != 0 || row.externalScalars != 0);
+
+ const auto reducer = mTopRowProjectFrom[lead.second].row;
+ if (
+ reducer.entryCount != 0 && // already have a reducer and...
+ reducer.entryCount < row.entryCount // ...it is sparser/better
+ ) {
+ const RowProjectFrom p = {1, row};
+ mBottomRowProjectFrom.push_back(p);
+ } else {
+ if (reducer.entryCount != 0) {
+ const RowProjectFrom p = {1, reducer};
+ mBottomRowProjectFrom.push_back(p);
}
- const auto rowIndex = left.rowCount();
- MATHICGB_ASSERT(rowIndex == right.rowCount());
- left.rowDone();
- right.rowDone();
-
- if (fromIt->first != 1) {
- MATHICGB_ASSERT(fromIt->first != 0);
- left.multiplyRow(rowIndex, fromIt->first, modulus);
- right.multiplyRow(rowIndex, fromIt->first, modulus);
- MATHICGB_ASSERT(left.rowBegin(rowIndex).scalar() == 1);
- }
-
- MATHICGB_ASSERT(left.rowCount() == right.rowCount());
+ const auto leadScalar = row.scalars != 0 ? row.scalars[lead.first] :
+ static_cast<SparseMatrix::Scalar>
+ (row.externalScalars[lead.first]);
+ MATHICGB_ASSERT(leadScalar != 0);
+ const auto inverse = leadScalar == 1 ?
+ 1 : modularInverse(leadScalar, modulus);
+ const RowProjectFrom p = {inverse, row};
+ mTopRowProjectFrom[lead.second] = p;
}
}
-
- const std::vector<monomial>& leftMonomials() const {
- return mLeftMonomials;
- }
-
- std::vector<monomial> moveLeftMonomials() {
- return std::move(mLeftMonomials);
+ }
+#ifdef MATHICGB_DEBUG
+ for (size_t i = 0; i < mTopRowProjectFrom.size(); ++i) {
+ const auto p = mTopRowProjectFrom[i];
+ MATHICGB_ASSERT(p.row.entryCount > 0);
+ MATHICGB_ASSERT(p.multiplyBy != 0);
+
+ // Find leading left entry.
+ ColIndex col = 0;
+ while (!mColProjectTo[p.row.indices[col]].isLeft) {
+ ++col;
+ MATHICGB_ASSERT(col < p.row.entryCount);
}
+ const auto projected = mColProjectTo[p.row.indices[col]];
- std::vector<monomial> moveRightMonomials() {
- return std::move(mRightMonomials);
- }
+ // The leading left entry of row i should be in column i.
+ MATHICGB_ASSERT(projected.index == i);
- private:
- std::vector<Projected> mProject;
- std::vector<monomial> mLeftMonomials;
- std::vector<monomial> mRightMonomials;
- };
-
- class TopBottomProjectionLate {
- public:
- TopBottomProjectionLate(
- const SparseMatrix& left,
- const SparseMatrix& right,
- const SparseMatrix::ColIndex leftColCount,
- const PolyRing& ring
- ):
- mModulus(static_cast<SparseMatrix::Scalar>(ring.charac()))
- {
- const auto noRow = std::numeric_limits<SparseMatrix::ColIndex>::max();
- mTopRows.resize(leftColCount, std::make_pair(0, noRow));
-
- MATHICGB_ASSERT(left.computeColCount() == leftColCount);
- MATHICGB_ASSERT(left.rowCount() >= leftColCount);
- MATHICGB_ASSERT(left.rowCount() == right.rowCount());
+ // After multiplication, the leading left scalar should be 1.
+ const auto leadScalar = p.row.scalars != 0 ?
+ p.row.scalars[col] : static_cast<Scalar>(p.row.externalScalars[col]);
+ MATHICGB_ASSERT(modularProduct(leadScalar, p.multiplyBy, modulus) == 1);
+ }
+ for (size_t i = 0; i < mBottomRowProjectFrom.size(); ++i) {
+ const auto p = mBottomRowProjectFrom[i];
+ MATHICGB_ASSERT(p.row.entryCount > 0);
+ MATHICGB_ASSERT(p.multiplyBy != 0);
+ }
+#endif
+}
- std::vector<SparseMatrix::ColIndex> topEntryCounts(leftColCount);
-
- const auto rowCount = left.rowCount();
- for (SparseMatrix::RowIndex row = 0; row < rowCount; ++row) {
- const auto leftEntryCount = left.entryCountInRow(row);
- const auto entryCount = leftEntryCount + right.entryCountInRow(row);
- MATHICGB_ASSERT(entryCount >= leftEntryCount); // no overflow
- if (entryCount == 0)
- continue; // ignore zero rows
- if (leftEntryCount == 0) {
- mBottomRows.push_back(std::make_pair(1, row)); //can't be top/reducer
- continue;
- }
- const auto lead = left.rowBegin(row).index();
- if (mTopRows[lead].second != noRow && topEntryCounts[lead]<entryCount)
- mBottomRows.push_back(std::make_pair(1, row)); //other reducer better
- else {
- if (mTopRows[lead].second != noRow)
- mBottomRows.push_back(std::make_pair(1, mTopRows[lead].second));
- topEntryCounts[lead] = entryCount;
- mTopRows[lead].second = row; // do scalar .first later
- }
- }
+void F4MatrixProjection::projectRows(
+ SparseMatrix&& in,
+ SparseMatrix& top,
+ SparseMatrix& bottom
+) {
+ const auto modulus = static_cast<Scalar>(ring().charac());
+
+ top.clear();
+ bottom.clear();
+
+ const auto rowCountTop =
+ static_cast<SparseMatrix::RowIndex>(mTopRows.size());
+ for (SparseMatrix::RowIndex toRow = 0; toRow < rowCountTop; ++toRow) {
+ top.appendRow(in, mTopRows[toRow].row);
+ if (mTopRows[toRow].multiplyBy != 1)
+ top.multiplyRow(toRow, mTopRows[toRow].multiplyBy, modulus);
+ }
- const auto modulus = static_cast<SparseMatrix::Scalar>(ring.charac());
- for (SparseMatrix::RowIndex r = 0; r < leftColCount; ++r) {
- const auto row = mTopRows[r].second;
- MATHICGB_ASSERT(row != noRow);
- MATHICGB_ASSERT(left.entryCountInRow(row) > 0);
- MATHICGB_ASSERT(left.rowBegin(row).index() == r);
- MATHICGB_ASSERT(left.rowBegin(row).scalar() != 0);
- MATHICGB_ASSERT(topEntryCounts[r] ==
- left.entryCountInRow(row) + right.entryCountInRow(row));
-
- const auto leadScalar = left.rowBegin(row).scalar();
- mTopRows[r].first = leadScalar == 1 ? 1 : // 1 is the common case
- modularInverse(leadScalar, modulus);
- }
+ const auto rowCountBottom =
+ static_cast<SparseMatrix::RowIndex>(mBottomRows.size());
+ for (SparseMatrix::RowIndex toRow = 0; toRow < rowCountBottom; ++toRow) {
+ bottom.appendRow(in, mBottomRows[toRow].row);
+ if (mBottomRows[toRow].multiplyBy != 1)
+ bottom.multiplyRow(toRow, mBottomRows[toRow].multiplyBy, modulus);
+ }
-#ifdef MATHICGB_DEBUG
- for (SparseMatrix::RowIndex r = 0; r < mBottomRows.size(); ++r) {
- const auto row = mBottomRows[r].second;
- MATHICGB_ASSERT(
- left.entryCountInRow(row) + right.entryCountInRow(row) > 0);
- MATHICGB_ASSERT(mBottomRows[r].first == 1);
- }
-#endif
- }
+ in.clear();
+}
- void project(
- SparseMatrix&& in,
- SparseMatrix& top,
- SparseMatrix& bottom
- ) {
- top.clear();
- bottom.clear();
-
- const auto rowCountTop =
- static_cast<SparseMatrix::RowIndex>(mTopRows.size());
- for (SparseMatrix::RowIndex toRow = 0; toRow < rowCountTop; ++toRow) {
- top.appendRow(in, mTopRows[toRow].second);
- if (mTopRows[toRow].first != 1)
- top.multiplyRow(toRow, mTopRows[toRow].first, mModulus);
- }
+void F4MatrixProjection::setupRowProjectionLate(
+ const SparseMatrix& left,
+ const SparseMatrix& right
+) {
+ const auto leftColCount = static_cast<ColIndex>(mLeftMonomials.size());
+ const auto modulus = static_cast<Scalar>(ring().charac());
+ const auto noRow = std::numeric_limits<ColIndex>::max();
+ {
+ const RowProjectFromLate init = {0, noRow};
+ mTopRows.resize(leftColCount, init);
+ }
- const auto rowCountBottom =
- static_cast<SparseMatrix::RowIndex>(mBottomRows.size());
- for (SparseMatrix::RowIndex toRow = 0; toRow < rowCountBottom; ++toRow) {
- bottom.appendRow(in, mBottomRows[toRow].second);
- if (mBottomRows[toRow].first != 1)
- bottom.multiplyRow(toRow, mBottomRows[toRow].first, mModulus);
+ MATHICGB_ASSERT(left.computeColCount() == leftColCount);
+ MATHICGB_ASSERT(left.rowCount() >= leftColCount);
+ MATHICGB_ASSERT(left.rowCount() == right.rowCount());
+
+ std::vector<SparseMatrix::ColIndex> topEntryCounts(leftColCount);
+
+ const auto rowCount = left.rowCount();
+ for (SparseMatrix::RowIndex row = 0; row < rowCount; ++row) {
+ const auto leftEntryCount = left.entryCountInRow(row);
+ const auto entryCount = leftEntryCount + right.entryCountInRow(row);
+ MATHICGB_ASSERT(entryCount >= leftEntryCount); // no overflow
+ if (entryCount == 0)
+ continue; // ignore zero rows
+ if (leftEntryCount == 0) {
+ RowProjectFromLate p = {1, row};
+ mBottomRows.push_back(p); // can't be top/reducer
+ continue;
+ }
+ const auto lead = left.rowBegin(row).index();
+ if (mTopRows[lead].row != noRow && topEntryCounts[lead]<entryCount) {
+ RowProjectFromLate p = {1, row};
+ mBottomRows.push_back(p); // other reducer better
+ } else {
+ if (mTopRows[lead].row != noRow) {
+ RowProjectFromLate p = {1, mTopRows[lead].row};
+ mBottomRows.push_back(p);
}
-
- in.clear();
+ topEntryCounts[lead] = entryCount;
+ mTopRows[lead].row = row; // do scalar .first later
}
+ }
- private:
- const SparseMatrix::Scalar mModulus;
- std::vector<std::pair<SparseMatrix::Scalar, SparseMatrix::RowIndex>> mTopRows;
- std::vector<std::pair<SparseMatrix::Scalar, SparseMatrix::RowIndex>> mBottomRows;
- };
-
- class TopBottomProjection {
- public:
- TopBottomProjection(
- const std::vector<F4ProtoMatrix*>& blocks,
- const LeftRightProjection& leftRight,
- const PolyRing& ring
- ):
- mReducerRows(leftRight.leftMonomials().size())
- {
- typedef SparseMatrix::RowIndex RowIndex;
- const auto noReducer = std::numeric_limits<RowIndex>::max();
- F4ProtoMatrix::Row noRow = {};
- noRow.indices = 0;
- const auto noCol = std::numeric_limits<SparseMatrix::ColIndex>::max();
-
- const auto modulus = static_cast<SparseMatrix::Scalar>(ring.charac());
-
- const auto end = blocks.end();
- for (auto it = blocks.begin(); it != end; ++it) {
- auto& block = **it;
- 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 = [&] {
- for (SparseMatrix::ColIndex col = 0; col < row.entryCount; ++col) {
- auto const translated = leftRight.project(row.indices[col]);
- if (translated.left)
- return std::make_pair(col, translated.index);
- }
- return std::make_pair(noCol, noCol);
- }();
- // Decide if this should be a reducer or reducee row.
- if (lead.second == noCol) {
- mReduceeRows.push_back(std::make_pair(1, row)); // no left entries
- continue;
- }
- MATHICGB_ASSERT(row.scalars != 0 || row.externalScalars != 0);
-
- const auto reducer = mReducerRows[lead.second].second;
- if (
- reducer.entryCount != 0 && // already have a reducer and...
- reducer.entryCount < row.entryCount // ...it is sparser/better
- )
- mReduceeRows.push_back(std::make_pair(1, row));
- else {
- if (reducer.entryCount != 0)
- mReduceeRows.push_back(std::make_pair(1, reducer));
- const auto leadScalar = row.scalars != 0 ? row.scalars[lead.first] :
- static_cast<SparseMatrix::Scalar>
- (row.externalScalars[lead.first]);
- MATHICGB_ASSERT(leadScalar != 0);
- const auto inverse = leadScalar == 1 ?
- 1 : modularInverse(leadScalar, modulus);
- mReducerRows[lead.second] = std::make_pair(inverse, row);
- }
- }
- }
+ for (SparseMatrix::RowIndex r = 0; r < leftColCount; ++r) {
+ const auto row = mTopRows[r].row;
+ MATHICGB_ASSERT(row != noRow);
+ MATHICGB_ASSERT(left.entryCountInRow(row) > 0);
+ MATHICGB_ASSERT(left.rowBegin(row).index() == r);
+ MATHICGB_ASSERT(left.rowBegin(row).scalar() != 0);
+ MATHICGB_ASSERT(topEntryCounts[r] ==
+ left.entryCountInRow(row) + right.entryCountInRow(row));
+
+ const auto leadScalar = left.rowBegin(row).scalar();
+ mTopRows[r].multiplyBy = leadScalar == 1 ? 1 : // 1 is the common case
+ modularInverse(leadScalar, modulus);
+ }
#ifdef MATHICGB_DEBUG
- for (size_t i = 0; i < mReducerRows.size(); ++i) {
- const auto row = mReducerRows[i];
- MATHICGB_ASSERT(row.second.entryCount > 0);
- for (SparseMatrix::ColIndex col = 0; ; ++col) {
- MATHICGB_ASSERT(col < row.second.entryCount);
- const auto projected = leftRight.project(row.second.indices[col]);
- if (projected.left) {
- MATHICGB_ASSERT(projected.index == i);
- const auto leadScalar = row.second.scalars != 0 ?
- row.second.scalars[col] :
- static_cast<SparseMatrix::Scalar>(row.second.externalScalars[col]);
- MATHICGB_ASSERT(modularProduct(leadScalar, row.first, modulus) == 1);
- break;
- }
- }
- }
- for (size_t i = 0; i < mReduceeRows.size(); ++i) {
- const auto row = mReduceeRows[i];
- MATHICGB_ASSERT(row.second.entryCount > 0);
- MATHICGB_ASSERT(row.first == 1);
+ for (SparseMatrix::RowIndex r = 0; r < mBottomRows.size(); ++r) {
+ const auto row = mBottomRows[r].row;
+ MATHICGB_ASSERT(
+ left.entryCountInRow(row) + right.entryCountInRow(row) > 0);
+ MATHICGB_ASSERT(mBottomRows[r].multiplyBy == 1);
}
#endif
- }
+}
- const std::vector<std::pair<SparseMatrix::Scalar, F4ProtoMatrix::Row>>& reducerRows() const {
- return mReducerRows;
+void F4MatrixProjection::projectRowsLeftRight(
+ const std::vector<RowProjectFrom>& from,
+ SparseMatrix& left,
+ SparseMatrix& right
+) const {
+ left.clear();
+ right.clear();
+ const auto fromEnd = from.end();
+ for (auto fromIt = from.begin(); fromIt != fromEnd; ++fromIt) {
+ const auto row = fromIt->row;
+ MATHICGB_ASSERT(row.entryCount != 0);
+ MATHICGB_ASSERT(row.scalars == 0 || row.externalScalars == 0);
+ const auto modulus = static_cast<SparseMatrix::Scalar>(ring().charac());
+
+ 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<SparseMatrix::Scalar>(*scalars);
+ const auto index = *indices;
+ MATHICGB_ASSERT(index < mColProjectTo.size());
+ const auto translated = mColProjectTo[index];
+ if (translated.isLeft)
+ left.appendEntry(translated.index, scalar);
+ else
+ right.appendEntry(translated.index, scalar);
+ }
+ } else {
+ auto indices = row.indices;
+ auto indicesEnd = row.indices + row.entryCount;
+ auto scalars = row.scalars;
+ for (; indices != indicesEnd; ++indices, ++scalars) {
+ const auto index = *indices;
+ MATHICGB_ASSERT(index < mColProjectTo.size());
+ const auto translated = mColProjectTo[index];
+ if (translated.isLeft)
+ left.appendEntry(translated.index, *scalars);
+ else
+ right.appendEntry(translated.index, *scalars);
+ }
}
-
- const std::vector<std::pair<SparseMatrix::Scalar, F4ProtoMatrix::Row>>& reduceeRows() const {
- return mReduceeRows;
+ const auto rowIndex = left.rowCount();
+ MATHICGB_ASSERT(rowIndex == right.rowCount());
+ left.rowDone();
+ right.rowDone();
+
+ if (fromIt->multiplyBy != 1) {
+ MATHICGB_ASSERT(fromIt->multiplyBy != 0);
+ left.multiplyRow(rowIndex, fromIt->multiplyBy, modulus);
+ right.multiplyRow(rowIndex, fromIt->multiplyBy, modulus);
+ MATHICGB_ASSERT(left.rowBegin(rowIndex).scalar() == 1);
}
- private:
- std::vector<std::pair<SparseMatrix::Scalar, F4ProtoMatrix::Row>> mReducerRows;
- std::vector<std::pair<SparseMatrix::Scalar, F4ProtoMatrix::Row>> mReduceeRows;
- };
+ MATHICGB_ASSERT(left.rowCount() == right.rowCount());
+ }
}
-
-
-QuadMatrix F4MatrixProjection::project(
- std::vector<char>&& isColumnToLeft,
- MonomialMap<ColIndex>&& map,
- std::vector<F4ProtoMatrix*>&& matrix
-) {
- const auto& ring = map.ring();
+
+void F4MatrixProjection::projectLeftRight(
+ const std::vector<F4ProtoMatrix*>& preBlocks,
+ SparseMatrix& left,
+ SparseMatrix& right
+) const {
+ left.clear();
+ right.clear();
+ const auto modulus = static_cast<SparseMatrix::Scalar>(ring().charac());
+
+ const auto end = preBlocks.end();
+ for (auto it = preBlocks.begin(); it != end; ++it) {
+ auto& block = **it;
+ const auto rowCount = block.rowCount();
+ for (SparseMatrix::RowIndex r = 0; r < rowCount; ++r) {
+ const auto row = block.row(r);
+ if (row.entryCount == 0)
+ continue;
+ 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<SparseMatrix::Scalar>(*scalars);
+ const auto index = *indices;
+ MATHICGB_ASSERT(index < mColProjectTo.size());
+ const auto translated = mColProjectTo[index];
+ if (translated.isLeft)
+ left.appendEntry(translated.index, scalar);
+ else
+ right.appendEntry(translated.index, scalar);
+ }
+ } else {
+ auto indices = row.indices;
+ auto indicesEnd = row.indices + row.entryCount;
+ auto scalars = row.scalars;
+ for (; indices != indicesEnd; ++indices, ++scalars) {
+ const auto index = *indices;
+ MATHICGB_ASSERT(index < mColProjectTo.size());
+ const auto translated = mColProjectTo[index];
+ if (translated.isLeft)
+ left.appendEntry(translated.index, *scalars);
+ else
+ right.appendEntry(translated.index, *scalars);
+ }
+ }
+ MATHICGB_ASSERT(left.rowCount() == right.rowCount());
+ left.rowDone();
+ right.rowDone();
+ }
+ }
+}
+
+QuadMatrix F4MatrixProjection::makeProjectionAndClear() {
QuadMatrix quadMatrix;
- // Create projections
- LeftRightProjection projection(isColumnToLeft, map);
if (true) {
- TopBottomProjection topBottom(matrix, projection, ring);
-
- // Project the pre-blocks into the matrix
- projection.project(topBottom.reducerRows(), quadMatrix.topLeft, quadMatrix.topRight, ring);
- projection.project(topBottom.reduceeRows(), quadMatrix.bottomLeft, quadMatrix.bottomRight, ring);
+ setupRowProjection();
+ projectRowsLeftRight
+ (mTopRowProjectFrom, quadMatrix.topLeft, quadMatrix.topRight);
+ projectRowsLeftRight
+ (mBottomRowProjectFrom, quadMatrix.bottomLeft, quadMatrix.bottomRight);
} else {
SparseMatrix left;
SparseMatrix right;
- projection.project(matrix, left, right, ring);
- TopBottomProjectionLate topBottom(left, right, left.computeColCount(), ring);
- topBottom.project(std::move(left), quadMatrix.topLeft, quadMatrix.bottomLeft);
- topBottom.project(std::move(right), quadMatrix.topRight, quadMatrix.bottomRight);
+ projectLeftRight(mMatrices, left, right);
+ setupRowProjectionLate(left, right);
+ projectRows(std::move(left), quadMatrix.topLeft, quadMatrix.bottomLeft);
+ projectRows(std::move(right), quadMatrix.topRight, quadMatrix.bottomRight);
}
- quadMatrix.ring = ˚
- quadMatrix.leftColumnMonomials = projection.moveLeftMonomials();
- quadMatrix.rightColumnMonomials = projection.moveRightMonomials();
+ quadMatrix.ring = &ring();
+ quadMatrix.leftColumnMonomials = std::move(mLeftMonomials);
+ quadMatrix.rightColumnMonomials = std::move(mRightMonomials);
return std::move(quadMatrix);
-}
+}
diff --git a/src/mathicgb/F4MatrixProjection.hpp b/src/mathicgb/F4MatrixProjection.hpp
index 3dcd2d9..ff9bcac 100644
--- a/src/mathicgb/F4MatrixProjection.hpp
+++ b/src/mathicgb/F4MatrixProjection.hpp
@@ -1,25 +1,85 @@
-#ifndef MATHICGB_F4_MATRIX_PROJECTION_GUARD
-#define MATHICGB_F4_MATRIX_PROJECTION_GUARD
-
-#include "QuadMatrix.hpp"
-#include "SparseMatrix.hpp"
-#include "F4ProtoMatrix.hpp"
-#include "MonomialMap.hpp"
-#include <vector>
-
-class F4MatrixProjection {
-public:
- typedef SparseMatrix::ColIndex ColIndex;
-
- F4MatrixProjection() {}
-
- QuadMatrix project(
- std::vector<char>&& isColumnToLeft,
- MonomialMap<ColIndex>&& map,
- std::vector<F4ProtoMatrix*>&& matrix
- );
-
-private:
-};
-
-#endif
+#ifndef MATHICGB_F4_MATRIX_PROJECTION_GUARD
+#define MATHICGB_F4_MATRIX_PROJECTION_GUARD
+
+#include "QuadMatrix.hpp"
+#include "SparseMatrix.hpp"
+#include "F4ProtoMatrix.hpp"
+#include "MonomialMap.hpp"
+#include "ScopeExit.hpp"
+#include <vector>
+
+class F4MatrixProjection {
+public:
+ typedef SparseMatrix::RowIndex RowIndex;
+ typedef SparseMatrix::ColIndex ColIndex;
+ typedef SparseMatrix::Scalar Scalar;
+
+ F4MatrixProjection(const PolyRing& ring, ColIndex colCount);
+
+ void addProtoMatrix(F4ProtoMatrix&& matrix) {mMatrices.push_back(&matrix);}
+
+ // No reference to mono is retained.
+ void addColumn(ColIndex index, const_monomial mono, const bool isLeft);
+
+ QuadMatrix makeProjectionAndClear();
+
+ const PolyRing& ring() const {return mRing;}
+
+private:
+ std::vector<F4ProtoMatrix*> mMatrices;
+
+ // *** Projection of columns
+ struct ColProjectTo {
+ ColIndex index;
+ bool isLeft;
+ };
+ std::vector<ColProjectTo> mColProjectTo;
+
+ // *** Projection: Simultaneously do let/right and row permutation
+ void setupRowProjection();
+ struct RowProjectFrom {
+ Scalar multiplyBy;
+ F4ProtoMatrix::Row row;
+ };
+ void projectRowsLeftRight(
+ const std::vector<RowProjectFrom>& from,
+ SparseMatrix& left,
+ SparseMatrix& right
+ ) const;
+ std::vector<RowProjectFrom> mTopRowProjectFrom;
+ std::vector<RowProjectFrom> mBottomRowProjectFrom;
+
+ // *** Projection: Do left/right, then permute rows
+ void setupRowProjectionLate(
+ const SparseMatrix& left,
+ const SparseMatrix& right
+ );
+ struct RowProjectFromLate {
+ Scalar multiplyBy;
+ RowIndex row;
+ };
+ void projectLeftRight(
+ const std::vector<F4ProtoMatrix*>& preblocks,
+ SparseMatrix& left,
+ SparseMatrix& right
+ ) const;
+ void projectRows(
+ SparseMatrix&& in,
+ SparseMatrix& top,
+ SparseMatrix& bottom
+ );
+ std::vector<RowProjectFromLate> mTopRows;
+ std::vector<RowProjectFromLate> mBottomRows;
+
+
+
+
+
+
+ std::vector<monomial> mLeftMonomials;
+ std::vector<monomial> mRightMonomials;
+
+ const PolyRing& mRing;
+};
+
+#endif
diff --git a/src/mathicgb/MonomialMap.hpp b/src/mathicgb/MonomialMap.hpp
index f137012..9b0c535 100755
--- a/src/mathicgb/MonomialMap.hpp
+++ b/src/mathicgb/MonomialMap.hpp
@@ -5,7 +5,7 @@
#include "Atomic.hpp"
#include "PolyRing.hpp"
#include <memtailor.h>
-#include <tbb/tbb.h>
+#include <tbb/tbb.h>
#include <limits>
#include <vector>
#include <algorithm>
diff --git a/src/mathicgb/ScopeExit.hpp b/src/mathicgb/ScopeExit.hpp
new file mode 100755
index 0000000..5ca876a
--- /dev/null
+++ b/src/mathicgb/ScopeExit.hpp
@@ -0,0 +1,98 @@
+#ifndef MATHICGB_SCOPE_EXIT_GUARD
+#define MATHICGB_SCOPE_EXIT_GUARD
+
+// Guard holds an action to call and calls it unless it has been released.
+template<class T>
+class Guard {
+public:
+ ~Guard() {
+ if (mOwning && mActive)
+ mAction();
+ }
+
+private:
+ friend struct GuardMaker;
+ Guard(T&& action, const bool& active):
+ mAction(std::move(action)), mOwning(true), mActive(active) {}
+
+ // Most compilers should elide the call to this construtor, but it must be
+ // here anyway and we should support even a crazy compiler that decides to
+ // call it.
+ Guard(Guard<T>&& guard):
+ mAction(std::move(guard.mAction)), mOwning(true), mActive(guard.mActive)
+ {
+ assert(guard.mActive);
+ guard.mOwning = false; // to avoid calling mAction twice
+ }
+
+ bool mOwning;
+ const bool& mActive;
+ const T mAction;
+};
+
+// The class user code interacts with to dismiss an action.
+class Dismisser {
+public:
+ Dismisser(bool& active): mActive(active) {}
+ void dismiss() {mActive = false;}
+
+private:
+ bool& mActive;
+};
+
+// Helper class that allows convenient syntax for the macro by overloading
+// operator+.
+struct GuardMaker {
+public:
+ GuardMaker(const bool& active): mActive(active) {}
+
+ template<class T>
+ Guard<T> operator+(T&& t) {return Guard<T>(std::forward<T>(t), mActive);}
+
+private:
+ const bool& mActive;
+};
+
+#define MYLIB__CAT_HELPER(A, B) A##B
+#define MYLIB__CAT(A, B) MYLIB__CAT_HELPER(A, B)
+#define MYLIB__UNIQUE(NAME) MYLIB__CAT(MyLib_,MYLIB__CAT(NAME,__LINE__))
+
+// Example, with no need to dismiss:
+// FILE* file = fopen("file.txt", "r");
+// MATHICGB_SCOPE_EXIT() {
+// fclose(file);
+// std::cout << "file closed";
+// };
+// // ...
+// return; // the file is closed
+//
+// Example, with need to dismiss:
+// v.push_back(5);
+// MATHICGB_SCOPE_EXIT(name) {v.pop_back();};
+// // ...
+// if (error)
+// return; // the pop_back is done
+// name.dismiss();
+// return; // the pop_back is not done
+//
+// The middle line is a no-op if the name parameter expands to nothing.
+// When NAME expands to nothing, we need the static_cast to prevent
+// the compiler from parsing the middle line as a redeclaration of
+// myBool. In the final line we use that a const reference keeps
+// temporary objects alive until the end of the scope. It would be correct
+// to copy, but this way we can keep the copy constructor private
+// and help out any compiler that has issue with eliding that copy.
+#define MATHICGB_SCOPE_EXIT(NAME) \
+ bool MYLIB__UNIQUE(active) = true; \
+ ::Dismisser NAME(static_cast<bool&>(MYLIB__UNIQUE(active))); \
+ const auto& MYLIB__UNIQUE(guard) = \
+ ::GuardMaker(MYLIB__UNIQUE(active)) + [&]
+
+// Without this pragma, MSVC will say
+// warning C4003: not enough actual parameters for macro 'MYLIB_SCOPE_EXIT'
+// when using MYLIB_SCOPE_EXIT without a name parameter. Not very happy about
+// turning off the warning. I wonder if there is a way to avoid the warning in
+// this case without turning it off.
+#pragma warning (disable: 4003)
+
+#endif
--
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