[mathicgb] 116/393: Fixed a matrix write bug and now the matrix action reduces the input file and either writes the reduced bottom right matrix out if a reference matrix does not exist or checks that the computed matrix matches the reference matrix if the reference matrix does exist.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:44 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 066c4f7b79664c9a3c45b3df612f441b68a46bda
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Tue Nov 20 16:26:57 2012 +0100
Fixed a matrix write bug and now the matrix action reduces the input file and either writes the reduced bottom right matrix out if a reference matrix does not exist or checks that the computed matrix matches the reference matrix if the reference matrix does exist.
---
.gitignore | 3 +-
Makefile.am | 13 +-
build/vs12/mathicgb-lib/mathicgb-lib.vcxproj | 2 +
.../vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters | 6 +
src/cli/MatrixAction.cpp | 93 +++-
src/mathicgb/CFile.cpp | 28 +
src/mathicgb/CFile.hpp | 37 ++
src/mathicgb/F4MatrixReducer.cpp | 12 +-
src/mathicgb/F4Reducer.cpp | 9 +-
src/mathicgb/PolyRing.hpp | 1 +
src/mathicgb/QuadMatrix.cpp | 616 ++++++++++-----------
src/mathicgb/QuadMatrix.hpp | 4 +-
src/mathicgb/SparseMatrix.cpp | 19 +-
src/mathicgb/SparseMatrix.hpp | 5 +
14 files changed, 505 insertions(+), 343 deletions(-)
diff --git a/.gitignore b/.gitignore
index 3e9d79e..aa5a783 100755
--- a/.gitignore
+++ b/.gitignore
@@ -1,8 +1,9 @@
# MathicGB generated files
*.gb
*.stats
-*.mat
*.qmat
+*.brmat
+*.rbrmat
# Compiled Object files
*.slo
diff --git a/Makefile.am b/Makefile.am
index 86cd0f9..386bde9 100755
--- a/Makefile.am
+++ b/Makefile.am
@@ -61,11 +61,8 @@ libmathicgb_la_SOURCES = src/mathicgb/BjarkeGeobucket2.cpp \
src/mathicgb/QuadMatrix.cpp src/mathicgb/F4MatrixReducer.cpp \
src/mathicgb/F4MatrixReducer.hpp src/mathicgb/MonomialMap.hpp \
src/mathicgb/RawVector.hpp src/mathicgb/Atomic.hpp \
- src/mathicgb/FixedSizeMonomialMap.hpp src/cli/CommonParams.hpp \
- src/cli/CommonParams.cpp src/cli/GBAction.hpp src/cli/GBAction.cpp \
- src/cli/GBCommonParams.hpp src/cli/GBCommonParams.cpp \
- src/cli/MatrixAction.cpp src/cli/MatrixAction.hpp \
- src/cli/SigGBAction.hpp src/cli/SigGBAction.cpp
+ src/mathicgb/FixedSizeMonomialMap.hpp src/mathicgb/CFile.hpp \
+ src/mathicgb/CFile.cpp
# When making a distribution file, Automake knows to include all files
# that are necessary to build the project. EXTRA_DIST specifies files
@@ -77,7 +74,11 @@ bin_PROGRAMS = mgb
# set up the tests. Listing the headers in sources ensure that those
# files are included in distributions.
mgb_CPPFLAGS = $(DEPS_CFLAGS)
-mgb_SOURCES = src/cli/GBMain.cpp
+mgb_SOURCES = src/cli/GBMain.cpp src/cli/CommonParams.hpp \
+ src/cli/CommonParams.cpp src/cli/GBAction.hpp src/cli/GBAction.cpp \
+ src/cli/GBCommonParams.hpp src/cli/GBCommonParams.cpp \
+ src/cli/MatrixAction.cpp src/cli/MatrixAction.hpp \
+ src/cli/SigGBAction.hpp src/cli/SigGBAction.cpp
mgb_LDADD = $(top_builddir)/libmathicgb.la
# set up tests to run on "make check"
diff --git a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj
index 3619f54..6b2adbb 100755
--- a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj
+++ b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj
@@ -442,6 +442,7 @@
<ClCompile Include="..\..\..\src\mathicgb\BjarkeGeobucket.cpp" />
<ClCompile Include="..\..\..\src\mathicgb\BjarkeGeobucket2.cpp" />
<ClCompile Include="..\..\..\src\mathicgb\BuchbergerAlg.cpp" />
+ <ClCompile Include="..\..\..\src\mathicgb\CFile.cpp" />
<ClCompile Include="..\..\..\src\mathicgb\ChainedHashTable.cpp" />
<ClCompile Include="..\..\..\src\mathicgb\DivisorLookup.cpp" />
<ClCompile Include="..\..\..\src\mathicgb\F4MatrixBuilder.cpp" />
@@ -480,6 +481,7 @@
<ClInclude Include="..\..\..\src\mathicgb\BjarkeGeobucket.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\BjarkeGeobucket2.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\BuchbergerAlg.hpp" />
+ <ClInclude Include="..\..\..\src\mathicgb\CFile.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\ChainedHashTable.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\DivisorLookup.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\DivLookup.hpp" />
diff --git a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters
index 4c2001b..a88b0f8 100755
--- a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters
+++ b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters
@@ -120,6 +120,9 @@
<ClCompile Include="..\..\..\src\mathicgb\TypicalReducer.cpp">
<Filter>Source Files</Filter>
</ClCompile>
+ <ClCompile Include="..\..\..\src\mathicgb\CFile.cpp">
+ <Filter>Source Files</Filter>
+ </ClCompile>
</ItemGroup>
<ItemGroup>
<ClInclude Include="..\..\..\src\mathicgb\BjarkeGeobucket.hpp">
@@ -275,5 +278,8 @@
<ClInclude Include="..\..\..\src\mathicgb\FixedSizeMonomialMap.h">
<Filter>Header Files</Filter>
</ClInclude>
+ <ClInclude Include="..\..\..\src\mathicgb\CFile.hpp">
+ <Filter>Header Files</Filter>
+ </ClInclude>
</ItemGroup>
</Project>
\ No newline at end of file
diff --git a/src/cli/MatrixAction.cpp b/src/cli/MatrixAction.cpp
index c5545e2..6b79cbd 100644
--- a/src/cli/MatrixAction.cpp
+++ b/src/cli/MatrixAction.cpp
@@ -3,11 +3,36 @@
#include "mathicgb/QuadMatrix.hpp"
#include "mathicgb/SparseMatrix.hpp"
+#include "mathicgb/F4MatrixReducer.hpp"
+#include "mathicgb/SparseMatrix.hpp"
+#include "mathicgb/CFile.hpp"
+#include <mathic.h>
#include <fstream>
#include <iostream>
+namespace {
+ static const char* QuadMatrixExtension = ".qmat";
+ static const char* LowerRightMatrixExtension = ".brmat";
+ static const char* ReducedLowerRightMatrixExtension = ".rbrmat";
+
+ /// Returns true if the file exists - or more precisely if it can be opened
+ /// for reading. Using this function to create a file only if it does not
+ /// exist implies a race condition in that the file could have been crated
+ /// after checking if it exists and before (re)creating it. Do not use this
+ /// approach if that is not acceptable. The advantage here is that this is
+ /// portable. There could be a solution with freopen, but unfortunately
+ /// freopen is allowed to fail on any change to the mode so it is not
+ /// a portable solution.
+ bool fileExists(const std::string fileName) {
+ return CFile(fileName, "r", CFile::NoThrowTag()).hasFile();
+ }
+}
+
MatrixAction::MatrixAction() {
- mParams.registerFileNameExtension(".mat");
+ mParams.registerFileNameExtension(QuadMatrixExtension);
+ mParams.registerFileNameExtension(LowerRightMatrixExtension);
+ mParams.registerFileNameExtension(ReducedLowerRightMatrixExtension);
+ mParams.registerFileNameExtension(".");
}
void MatrixAction::directOptions(
@@ -19,24 +44,66 @@ void MatrixAction::directOptions(
void MatrixAction::performAction() {
mParams.perform();
+ const auto fileNameStem = mParams.inputFileNameStem();
+ const auto extension = mParams.inputFileNameExtension();
+ const auto quadFileName = fileNameStem + QuadMatrixExtension;
+ const auto lowerRightFileName = fileNameStem + LowerRightMatrixExtension;
+ const auto reducedLowerRightFileName =
+ fileNameStem + ReducedLowerRightMatrixExtension;
+ std::string inputFileName;
- QuadMatrix matrix;
+ SparseMatrix lowerRightMatrix;
+ SparseMatrix::Scalar modulus;
+ if (
+ extension == QuadMatrixExtension ||
+ extension == "." ||
+ extension == ""
+ ) {
+ inputFileName = quadFileName;
+ CFile file(quadFileName, "rb");
+ QuadMatrix matrix;
+ modulus = matrix.read(file.handle());
+ fclose(file.handle());
+ // @todo: F4MatrixReducer should not take a PolyRing parameter.
+ PolyRing ring(modulus, 0, 0);
+ F4MatrixReducer reducer(ring);
+ // @todo: only reduce down to D, do not reduce D itself
+ lowerRightMatrix = reducer.reduce(matrix);
- // @todo: fix leak of file on exception
+ if (!fileExists(lowerRightFileName)) {
+ CFile file(lowerRightFileName, "wb");
+ lowerRightMatrix.write(modulus, file.handle());
+ }
+ } else if (extension == LowerRightMatrixExtension) {
+ inputFileName = lowerRightFileName;
+ CFile file(lowerRightFileName, "rb");
+ modulus = lowerRightMatrix.read(file.handle());
+ } else {
+ mathic::reportError
+ ("Unknown input file extension of " + mParams.mInputFile.value());
+ }
- // read matrix
- SparseMatrix::Scalar modulus;
{
- FILE* file = fopen((mParams.inputFileNameStem() + ".mat").c_str(), "rb");
- modulus = matrix.read(file);
- fclose(file);
+ // @todo: expose D -> reduced D code and call it here
+ //PolyRing ring(modulus, 0, 0);
+ //F4MatrixReducer reducer(ring);
}
+ lowerRightMatrix.sortRowsByIncreasingPivots();
- // write matrix
- {
- FILE* file = fopen((mParams.inputFileNameStem() + ".out").c_str(), "wb");
- matrix.write(modulus, file);
- fclose(file);
+ if (!fileExists(reducedLowerRightFileName)) {
+ CFile file(reducedLowerRightFileName.c_str(), "wb");
+ lowerRightMatrix.write(modulus, file.handle());
+ } else {
+ SparseMatrix referenceMatrix;
+ CFile file(reducedLowerRightFileName.c_str(), "rb");
+ referenceMatrix.read(file.handle());
+ if (lowerRightMatrix != referenceMatrix) {
+ std::cerr << "Reducing " << inputFileName
+ << " does not yield the matrix " << reducedLowerRightFileName << ".\n";
+ } else if (tracingLevel > 0) {
+ std::cerr << "Match for " << inputFileName
+ << " -> " << ReducedLowerRightMatrixExtension << ".\n";
+ }
}
}
diff --git a/src/mathicgb/CFile.cpp b/src/mathicgb/CFile.cpp
new file mode 100755
index 0000000..b882aa0
--- /dev/null
+++ b/src/mathicgb/CFile.cpp
@@ -0,0 +1,28 @@
+#include "stdinc.h"
+#include "CFile.hpp"
+
+#include <mathic.h>
+#include <sstream>
+
+CFile::CFile(const std::string& fileName, const char* mode, NoThrowTag):
+ mFile(fopen(fileName.c_str(), mode)
+) {}
+
+CFile::CFile(const std::string& fileName, const char* mode):
+ mFile(fopen(fileName.c_str(), mode)
+) {
+ if (mFile == 0) {
+ std::ostringstream error;
+ error << "Could not open file " << fileName << " in mode " << mode << '.';
+ mathic::reportError(error.str());
+ }
+}
+
+CFile::~CFile() {
+ close();
+}
+
+void CFile::close() {
+ if (mFile != 0)
+ fclose(mFile);
+}
\ No newline at end of file
diff --git a/src/mathicgb/CFile.hpp b/src/mathicgb/CFile.hpp
new file mode 100755
index 0000000..1474443
--- /dev/null
+++ b/src/mathicgb/CFile.hpp
@@ -0,0 +1,37 @@
+#ifndef MATHICGB_C_FILE_GUARD
+#define MATHICGB_C_FILE_GUARD
+
+#include <string>
+#include <cstdio>
+
+/// RAII handle for a C FILE*.
+///
+/// The purpose of using the C IO interface instead of iostreams is that the
+/// former is faster to a ridiculous degree. This class wraps the C IO
+/// interface to be more useful in a C++ context. For example the file is
+/// automatically closed in the destructor and if the file cannot be opened
+/// then an exception is thrown instead of returning a null pointer.
+class CFile {
+public:
+ struct NoThrowTag {};
+
+ /// Sets the handle to null if the file cannot be opened - does not
+ /// throw an exception. The purpose of the NoTrowTag parameter is only
+ /// to indicate that no exception should be thrown on error.
+ CFile(const std::string& fileName, const char* mode, NoThrowTag);
+
+ /// Opens the file and throws an exception if the file cannot be opened.
+ CFile(const std::string& fileName, const char* mode);
+
+ ~CFile();
+
+ bool hasFile() const {return mFile != 0;}
+
+ FILE* handle() {return mFile;}
+ void close();
+
+private:
+ FILE* mFile;
+};
+
+#endif
diff --git a/src/mathicgb/F4MatrixReducer.cpp b/src/mathicgb/F4MatrixReducer.cpp
index 28376df..7caffb0 100755
--- a/src/mathicgb/F4MatrixReducer.cpp
+++ b/src/mathicgb/F4MatrixReducer.cpp
@@ -132,10 +132,10 @@ namespace {
const SparseMatrix& reduceByLeft = qm.topLeft;
const SparseMatrix& reduceByRight = qm.topRight;
- const auto leftColCount =
- static_cast<SparseMatrix::ColIndex>(qm.leftColumnMonomials.size());
- const auto rightColCount =
- static_cast<SparseMatrix::ColIndex>(qm.rightColumnMonomials.size());
+ const auto leftColCount = qm.computeLeftColCount();
+// static_cast<SparseMatrix::ColIndex>(qm.leftColumnMonomials.size());
+ const auto rightColCount = static_cast<SparseMatrix::ColIndex>(qm.computeRightColCount());
+// static_cast<SparseMatrix::ColIndex>(qm.rightColumnMonomials.size());
MATHICGB_ASSERT(leftColCount == reduceByLeft.rowCount());
const auto pivotCount = leftColCount;
const auto rowCount = toReduceLeft.rowCount();
@@ -373,8 +373,8 @@ SparseMatrix F4MatrixReducer::reduce(const QuadMatrix& matrix) {
if (tracingLevel >= 3)
matrix.printSizes(std::cerr);
- const auto rightColCount =
- static_cast<SparseMatrix::ColIndex>(matrix.rightColumnMonomials.size());
+ const auto rightColCount = matrix.computeRightColCount();
+ //static_cast<SparseMatrix::ColIndex>(matrix.rightColumnMonomials.size());
SparseMatrix newPivots(::reduce(matrix, mModulus));
::reduceToEchelonForm(newPivots, rightColCount, mModulus);
return std::move(newPivots);
diff --git a/src/mathicgb/F4Reducer.cpp b/src/mathicgb/F4Reducer.cpp
index 5941f48..4fa2b3c 100755
--- a/src/mathicgb/F4Reducer.cpp
+++ b/src/mathicgb/F4Reducer.cpp
@@ -201,12 +201,9 @@ void F4Reducer::saveMatrix(const QuadMatrix& matrix) {
return;
++mMatrixSaveCount;
std::ostringstream fileName;
- fileName << mStoreToFile << '-' << mMatrixSaveCount << ".mat";
- if (tracingLevel > 2) {
- std::cerr << "F4Reducer: Saving "
- << mathic::ColumnPrinter::commafy(entryCount)
- << " entry matrix to " << fileName.str() << '\n';
- }
+ fileName << mStoreToFile << '-' << mMatrixSaveCount << ".qmat";
+ if (tracingLevel > 2)
+ std::cerr << "F4Reducer: Saving matrix to " << fileName.str() << '\n';
FILE* file = fopen(fileName.str().c_str(), "wb");
// @todo: fix leak of file on exception
matrix.write(static_cast<SparseMatrix::Scalar>(mRing.charac()), file);
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index 90bd347..ca0e3ba 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -11,6 +11,7 @@
#include <memtailor.h>
#include <cstdio>
#include <cstring>
+#include <limits>
#define LT (-1)
#define EQ 0
diff --git a/src/mathicgb/QuadMatrix.cpp b/src/mathicgb/QuadMatrix.cpp
index b2315cd..3ec8bd7 100755
--- a/src/mathicgb/QuadMatrix.cpp
+++ b/src/mathicgb/QuadMatrix.cpp
@@ -1,321 +1,321 @@
-#include "stdinc.h"
-#include "QuadMatrix.hpp"
-
-#include <mathic.h>
-#include <tbb/tbb.h>
-#include <ostream>
-#include <sstream>
-
-bool QuadMatrix::debugAssertValid() const {
-#ifndef MATHICGB_DEBUG
- return true;
-#else
- 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
-}
-
-void QuadMatrix::print(std::ostream& out) const {
- MATHICGB_ASSERT(debugAssertValid());
-
- typedef SparseMatrix::ColIndex ColIndex;
- mathic::ColumnPrinter printer;
- printer.addColumn(true, "", "");
- printer.addColumn(true, " | ", "");
-
- // column monomials
- out << "Left columns:";
- const auto leftColCount = leftColumnMonomials.size();
- for (ColIndex leftCol = 0; leftCol < leftColCount; ++leftCol) {
- out << ' ';
- ring->monomialDisplay(out, leftColumnMonomials[leftCol], false, true);
- }
-
- out << "\nRight columns:";
- const auto rightColCount = rightColumnMonomials.size();
- for (ColIndex rightCol = 0; rightCol < rightColCount; ++rightCol) {
- out << ' ';
- ring->monomialDisplay(out, rightColumnMonomials[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;
-}
-
+#include "stdinc.h"
+#include "QuadMatrix.hpp"
+
+#include <mathic.h>
+#include <tbb/tbb.h>
+#include <ostream>
+#include <sstream>
+
+bool QuadMatrix::debugAssertValid() const {
+#ifndef MATHICGB_DEBUG
+ return true;
+#else
+ 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
+}
+
+void QuadMatrix::print(std::ostream& out) const {
+ MATHICGB_ASSERT(debugAssertValid());
+
+ typedef SparseMatrix::ColIndex ColIndex;
+ mathic::ColumnPrinter printer;
+ printer.addColumn(true, "", "");
+ printer.addColumn(true, " | ", "");
+
+ // column monomials
+ out << "Left columns:";
+ const auto leftColCount = leftColumnMonomials.size();
+ for (ColIndex leftCol = 0; leftCol < leftColCount; ++leftCol) {
+ out << ' ';
+ ring->monomialDisplay(out, leftColumnMonomials[leftCol], false, true);
+ }
+
+ out << "\nRight columns:";
+ const auto rightColCount = rightColumnMonomials.size();
+ for (ColIndex rightCol = 0; rightCol < rightColCount; ++rightCol) {
+ out << ' ';
+ ring->monomialDisplay(out, rightColumnMonomials[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;
+}
+
size_t QuadMatrix::rowCount() const {
return topLeft.rowCount() + bottomLeft.rowCount();
}
-size_t QuadMatrix::computeLeftColCount() const {
+SparseMatrix::ColIndex QuadMatrix::computeLeftColCount() const {
return std::max(topLeft.computeColCount(), bottomLeft.computeColCount());
}
-size_t QuadMatrix::computeRightColCount() const {
+SparseMatrix::ColIndex QuadMatrix::computeRightColCount() const {
return std::max(topRight.computeColCount(), bottomRight.computeColCount());
}
-size_t QuadMatrix::entryCount() const {
- return
- topLeft.entryCount() + topRight.entryCount() +
- bottomLeft.entryCount() + bottomRight.entryCount();
-}
-
-std::string QuadMatrix::toString() const {
- std::ostringstream out;
- print(out);
- return out.str();
-}
-
-size_t QuadMatrix::memoryUse() const {
- return topLeft.memoryUse() + topRight.memoryUse() +
- bottomLeft.memoryUse() + bottomRight.memoryUse();
-}
-
-size_t QuadMatrix::memoryUseTrimmed() const {
- return topLeft.memoryUseTrimmed() + topRight.memoryUseTrimmed() +
- bottomLeft.memoryUseTrimmed() + bottomRight.memoryUseTrimmed();
-}
-
-void QuadMatrix::printSizes(std::ostream& out) const {
- typedef mathic::ColumnPrinter ColPr;
-
- ColPr pr;
- pr.addColumn(false, " ", "");
- pr.addColumn(false, "", "");
- pr.addColumn(false, "", "");
- const char* const line = "----------";
-
- pr[0] << '\n';
- pr[1] << ColPr::commafy(leftColumnMonomials.size()) << " \n";
- pr[2] << ColPr::commafy(rightColumnMonomials.size()) << " \n";
-
- pr[0] << "/\n";
- pr[1] << line << "|\n";
- pr[2] << line << "\\\n";
-
- pr[0] << ColPr::commafy(topLeft.rowCount()) << " |\n";
- pr[1] << ColPr::commafy(topLeft.entryCount()) << " |\n";
- pr[2] << ColPr::commafy(topRight.entryCount()) << " |\n";
-
- pr[0] << "|\n";
- pr[1] << ColPr::bytesInUnit(topLeft.memoryUse()) << " |\n";
- pr[2] << ColPr::bytesInUnit(topRight.memoryUse()) << " |\n";
-
- pr[0] << "|\n";
- pr[1] << ColPr::percent(topLeft.memoryUse(), memoryUse()) << " |\n";
- pr[2] << ColPr::percent(topRight.memoryUse(), memoryUse()) << " |\n";
-
- pr[0] << "|\n";
- pr[1] << line << "|\n";
- pr[2] << line << "|\n";
-
- pr[0] << ColPr::commafy(bottomLeft.rowCount()) << " |\n";
- pr[1] << ColPr::commafy(bottomLeft.entryCount()) << " |\n";
- pr[2] << ColPr::commafy(bottomRight.entryCount()) << " |\n";
-
- pr[0] << "|\n";
- pr[1] << ColPr::bytesInUnit(bottomLeft.memoryUse()) << " |\n";
- pr[2] << ColPr::bytesInUnit(bottomRight.memoryUse()) << " |\n";
-
- pr[0] << "|\n";
- pr[1] << ColPr::percent(bottomLeft.memoryUse(), memoryUse()) << " |\n";
- pr[2] << ColPr::percent(bottomRight.memoryUse(), memoryUse()) << " |\n";
-
- pr[0] << "\\\n";
- pr[1] << line << "|\n";
- pr[2] << line << "/\n";
-
- out << '\n' << pr
- << "Total memory: " << ColPr::bytesInUnit(memoryUse()) << " ("
- << ColPr::percent(memoryUseTrimmed(), memoryUse())
- << " in use)\n";
-}
-
-QuadMatrix QuadMatrix::toCanonical() const {
- class RowComparer {
- public:
- RowComparer(const SparseMatrix& matrix): mMatrix(matrix) {}
- bool operator()(size_t a, size_t b) const {
- // if you need this to work for empty rows or identical leading columns
- // then update this code.
- MATHICGB_ASSERT(!mMatrix.emptyRow(a));
- MATHICGB_ASSERT(!mMatrix.emptyRow(b));
- auto itA = mMatrix.rowBegin(a);
- const auto endA = mMatrix.rowEnd(a);
- auto itB = mMatrix.rowBegin(b);
- const auto endB = mMatrix.rowEnd(b);
- for (; itA != endA; ++itA, ++itB) {
- if (itB == endB)
- return true;
-
- if (itA.index() > itB.index())
- return true;
- if (itA.index() < itB.index())
- return false;
-
- if (itA.scalar() > itB.scalar())
- return false;
- if (itA.scalar() < itB.scalar())
- return true;
- }
- return false;
- }
-
- private:
- const SparseMatrix& mMatrix;
- };
-
- const auto leftColCount = leftColumnMonomials.size();
- const auto rightColCount = rightColumnMonomials.size();
-
- // 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)
- rows.push_back(row);
- {
- RowComparer comparer(topLeft);
- std::sort(rows.begin(), rows.end(), comparer);
- }
-
- matrix.topLeft.clear();
- matrix.topRight.clear();
- for (size_t 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)
- rows.push_back(row);
- {
- RowComparer comparer(bottomLeft);
- std::sort(rows.begin(), rows.end(), comparer);
- }
-
- matrix.bottomLeft.clear();
- matrix.bottomRight.clear();
- for (size_t i = 0; i < rows.size(); ++i) {
- matrix.bottomLeft.appendRow(bottomLeft, rows[i]);
- matrix.bottomRight.appendRow(bottomRight, rows[i]);
- }
- }
-
- matrix.leftColumnMonomials = leftColumnMonomials;
- matrix.rightColumnMonomials = rightColumnMonomials;
- matrix.ring = ring;
-
- return std::move(matrix);
-}
-
-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() {
- typedef SparseMatrix::ColIndex ColIndex;
- std::vector<ColIndex> leftPermutation;
- std::vector<ColIndex> rightPermutation;
-
- tbb::parallel_for(0, 2, 1, [&](int i) {
- if (i == 0)
- leftPermutation =
- sortColumnMonomialsAndMakePermutation(leftColumnMonomials, *ring);
- else
- rightPermutation =
- sortColumnMonomialsAndMakePermutation(rightColumnMonomials, *ring);
- });
-
- tbb::parallel_for(0, 4, 1, [&](int 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);
- }
- });
-}
-
+size_t QuadMatrix::entryCount() const {
+ return
+ topLeft.entryCount() + topRight.entryCount() +
+ bottomLeft.entryCount() + bottomRight.entryCount();
+}
+
+std::string QuadMatrix::toString() const {
+ std::ostringstream out;
+ print(out);
+ return out.str();
+}
+
+size_t QuadMatrix::memoryUse() const {
+ return topLeft.memoryUse() + topRight.memoryUse() +
+ bottomLeft.memoryUse() + bottomRight.memoryUse();
+}
+
+size_t QuadMatrix::memoryUseTrimmed() const {
+ return topLeft.memoryUseTrimmed() + topRight.memoryUseTrimmed() +
+ bottomLeft.memoryUseTrimmed() + bottomRight.memoryUseTrimmed();
+}
+
+void QuadMatrix::printSizes(std::ostream& out) const {
+ typedef mathic::ColumnPrinter ColPr;
+
+ ColPr pr;
+ pr.addColumn(false, " ", "");
+ pr.addColumn(false, "", "");
+ pr.addColumn(false, "", "");
+ const char* const line = "----------";
+
+ pr[0] << '\n';
+ pr[1] << ColPr::commafy(leftColumnMonomials.size()) << " \n";
+ pr[2] << ColPr::commafy(rightColumnMonomials.size()) << " \n";
+
+ pr[0] << "/\n";
+ pr[1] << line << "|\n";
+ pr[2] << line << "\\\n";
+
+ pr[0] << ColPr::commafy(topLeft.rowCount()) << " |\n";
+ pr[1] << ColPr::commafy(topLeft.entryCount()) << " |\n";
+ pr[2] << ColPr::commafy(topRight.entryCount()) << " |\n";
+
+ pr[0] << "|\n";
+ pr[1] << ColPr::bytesInUnit(topLeft.memoryUse()) << " |\n";
+ pr[2] << ColPr::bytesInUnit(topRight.memoryUse()) << " |\n";
+
+ pr[0] << "|\n";
+ pr[1] << ColPr::percent(topLeft.memoryUse(), memoryUse()) << " |\n";
+ pr[2] << ColPr::percent(topRight.memoryUse(), memoryUse()) << " |\n";
+
+ pr[0] << "|\n";
+ pr[1] << line << "|\n";
+ pr[2] << line << "|\n";
+
+ pr[0] << ColPr::commafy(bottomLeft.rowCount()) << " |\n";
+ pr[1] << ColPr::commafy(bottomLeft.entryCount()) << " |\n";
+ pr[2] << ColPr::commafy(bottomRight.entryCount()) << " |\n";
+
+ pr[0] << "|\n";
+ pr[1] << ColPr::bytesInUnit(bottomLeft.memoryUse()) << " |\n";
+ pr[2] << ColPr::bytesInUnit(bottomRight.memoryUse()) << " |\n";
+
+ pr[0] << "|\n";
+ pr[1] << ColPr::percent(bottomLeft.memoryUse(), memoryUse()) << " |\n";
+ pr[2] << ColPr::percent(bottomRight.memoryUse(), memoryUse()) << " |\n";
+
+ pr[0] << "\\\n";
+ pr[1] << line << "|\n";
+ pr[2] << line << "/\n";
+
+ out << '\n' << pr
+ << "Total memory: " << ColPr::bytesInUnit(memoryUse()) << " ("
+ << ColPr::percent(memoryUseTrimmed(), memoryUse())
+ << " in use)\n";
+}
+
+QuadMatrix QuadMatrix::toCanonical() const {
+ class RowComparer {
+ public:
+ RowComparer(const SparseMatrix& matrix): mMatrix(matrix) {}
+ bool operator()(size_t a, size_t b) const {
+ // if you need this to work for empty rows or identical leading columns
+ // then update this code.
+ MATHICGB_ASSERT(!mMatrix.emptyRow(a));
+ MATHICGB_ASSERT(!mMatrix.emptyRow(b));
+ auto itA = mMatrix.rowBegin(a);
+ const auto endA = mMatrix.rowEnd(a);
+ auto itB = mMatrix.rowBegin(b);
+ const auto endB = mMatrix.rowEnd(b);
+ for (; itA != endA; ++itA, ++itB) {
+ if (itB == endB)
+ return true;
+
+ if (itA.index() > itB.index())
+ return true;
+ if (itA.index() < itB.index())
+ return false;
+
+ if (itA.scalar() > itB.scalar())
+ return false;
+ if (itA.scalar() < itB.scalar())
+ return true;
+ }
+ return false;
+ }
+
+ private:
+ const SparseMatrix& mMatrix;
+ };
+
+ const auto leftColCount = leftColumnMonomials.size();
+ const auto rightColCount = rightColumnMonomials.size();
+
+ // 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)
+ rows.push_back(row);
+ {
+ RowComparer comparer(topLeft);
+ std::sort(rows.begin(), rows.end(), comparer);
+ }
+
+ matrix.topLeft.clear();
+ matrix.topRight.clear();
+ for (size_t 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)
+ rows.push_back(row);
+ {
+ RowComparer comparer(bottomLeft);
+ std::sort(rows.begin(), rows.end(), comparer);
+ }
+
+ matrix.bottomLeft.clear();
+ matrix.bottomRight.clear();
+ for (size_t i = 0; i < rows.size(); ++i) {
+ matrix.bottomLeft.appendRow(bottomLeft, rows[i]);
+ matrix.bottomRight.appendRow(bottomRight, rows[i]);
+ }
+ }
+
+ matrix.leftColumnMonomials = leftColumnMonomials;
+ matrix.rightColumnMonomials = rightColumnMonomials;
+ matrix.ring = ring;
+
+ return std::move(matrix);
+}
+
+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() {
+ typedef SparseMatrix::ColIndex ColIndex;
+ std::vector<ColIndex> leftPermutation;
+ std::vector<ColIndex> rightPermutation;
+
+ tbb::parallel_for(0, 2, 1, [&](int i) {
+ if (i == 0)
+ leftPermutation =
+ sortColumnMonomialsAndMakePermutation(leftColumnMonomials, *ring);
+ else
+ rightPermutation =
+ sortColumnMonomialsAndMakePermutation(rightColumnMonomials, *ring);
+ });
+
+ tbb::parallel_for(0, 4, 1, [&](int 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);
+ }
+ });
+}
+
void QuadMatrix::write(
const SparseMatrix::Scalar modulus,
FILE* file
diff --git a/src/mathicgb/QuadMatrix.hpp b/src/mathicgb/QuadMatrix.hpp
index 4e3d84f..8d27911 100755
--- a/src/mathicgb/QuadMatrix.hpp
+++ b/src/mathicgb/QuadMatrix.hpp
@@ -62,10 +62,10 @@ public:
size_t rowCount() const;
/// Return the number of left columns.
- size_t computeLeftColCount() const;
+ SparseMatrix::ColIndex computeLeftColCount() const;
/// Return the number of right columns.
- size_t computeRightColCount() const;
+ SparseMatrix::ColIndex computeRightColCount() const;
void write(SparseMatrix::Scalar modulus, FILE* file) const;
diff --git a/src/mathicgb/SparseMatrix.cpp b/src/mathicgb/SparseMatrix.cpp
index 9a24622..87f0e59 100755
--- a/src/mathicgb/SparseMatrix.cpp
+++ b/src/mathicgb/SparseMatrix.cpp
@@ -153,6 +153,23 @@ void SparseMatrix::swap(SparseMatrix& matrix) {
swap(mMemoryQuantum, matrix.mMemoryQuantum);
}
+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) {
+ if (entryCountInRow(row) != matrix.entryCountInRow(row))
+ return false;
+ const auto end = rowEnd(row);
+ auto it = rowBegin(row);
+ auto matrixIt = matrix.rowBegin(row);
+ for (auto it = rowBegin(row); it != end; ++it, ++matrixIt)
+ if (*it != *matrixIt)
+ return false;
+ }
+ return true;
+}
+
SparseMatrix::ColIndex SparseMatrix::computeColCount() const {
// Obviously this can be done faster, but there has not been a need for that
// so far.
@@ -400,7 +417,7 @@ void SparseMatrix::write(const Scalar modulus, FILE* file) const {
// write scalars
for (SparseMatrix::RowIndex row = 0; row < storedRowCount; ++row)
- fwrite(&rowBegin(row).index(), sizeof(uint16), entryCountInRow(row), file);
+ fwrite(&rowBegin(row).scalar(), sizeof(uint16), entryCountInRow(row), file);
// write indices
for (SparseMatrix::RowIndex row = 0; row < storedRowCount; ++row)
diff --git a/src/mathicgb/SparseMatrix.hpp b/src/mathicgb/SparseMatrix.hpp
index 9980ba8..af19b30 100755
--- a/src/mathicgb/SparseMatrix.hpp
+++ b/src/mathicgb/SparseMatrix.hpp
@@ -74,6 +74,11 @@ public:
SparseMatrix& operator=(const SparseMatrix&);
void swap(SparseMatrix& matrix);
+ bool operator==(const SparseMatrix& matrix) const;
+ bool operator!=(const SparseMatrix& matrix) const {
+ return !(*this == matrix);
+ }
+
// Removes all rows from *this.
void clear();
--
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