[mathicgb] 69/393: Added a RawVector class that works like std::vector but has manual memory management. This allowed to optimize appendEntry on SparseMatrix which reduced hcyc8 from 7.7s to 7.4s for a 3.8% performance increase.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:34 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 f5efe9ba0fcd38f3cd323965782ccf5365a0ae26
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Fri Oct 19 18:33:57 2012 +0200
Added a RawVector class that works like std::vector but has manual memory management. This allowed to optimize appendEntry on SparseMatrix which reduced hcyc8 from 7.7s to 7.4s for a 3.8% performance increase.
---
Makefile.am | 3 +-
build/vs12/mathicgb-lib/mathicgb-lib.vcxproj | 1 +
.../vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters | 3 +
src/mathicgb/F4Reducer.cpp | 1 +
src/mathicgb/MonomialMap.hpp | 545 +++++++++++----------
src/mathicgb/QuadMatrix.cpp | 2 +-
src/mathicgb/QuadMatrix.hpp | 8 +-
src/mathicgb/RawVector.hpp | 286 +++++++++++
src/mathicgb/SparseMatrix.cpp | 6 +-
src/mathicgb/SparseMatrix.hpp | 98 +++-
10 files changed, 652 insertions(+), 301 deletions(-)
diff --git a/Makefile.am b/Makefile.am
index 1e0ed08..546b0a5 100755
--- a/Makefile.am
+++ b/Makefile.am
@@ -61,7 +61,8 @@ libmathicgb_la_SOURCES = src/mathicgb/BjarkeGeobucket2.cpp \
src/mathicgb/F4Reducer.cpp src/mathicgb/F4MatrixBuilder.hpp \
src/mathicgb/F4MatrixBuilder.cpp src/mathicgb/QuadMatrix.hpp \
src/mathicgb/QuadMatrix.cpp src/mathicgb/F4MatrixReducer.cpp \
- src/mathicgb/F4MatrixReducer.hpp src/mathicgb/MonomialMap.hpp
+ src/mathicgb/F4MatrixReducer.hpp src/mathicgb/MonomialMap.hpp \
+ src/mathicgb/RawVector.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 d853efd..433b5cb 100755
--- a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj
+++ b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj
@@ -468,6 +468,7 @@
<ClInclude Include="..\..\..\src\mathicgb\PolyRing.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\QuadMatrix.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\QuadMatrixBuilder.hpp" />
+ <ClInclude Include="..\..\..\src\mathicgb\RawVector.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\Reducer.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\ReducerDedup.hpp" />
<ClInclude Include="..\..\..\src\mathicgb\ReducerHash.hpp" />
diff --git a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters
index 84c8e87..18837e9 100755
--- a/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters
+++ b/build/vs12/mathicgb-lib/mathicgb-lib.vcxproj.filters
@@ -278,5 +278,8 @@
<ClInclude Include="..\..\..\src\mathicgb\MonomialMap.hpp">
<Filter>Header Files</Filter>
</ClInclude>
+ <ClInclude Include="..\..\..\src\mathicgb\RawVector.hpp">
+ <Filter>Source Files</Filter>
+ </ClInclude>
</ItemGroup>
</Project>
\ No newline at end of file
diff --git a/src/mathicgb/F4Reducer.cpp b/src/mathicgb/F4Reducer.cpp
index bfe9b36..7046ddf 100755
--- a/src/mathicgb/F4Reducer.cpp
+++ b/src/mathicgb/F4Reducer.cpp
@@ -45,6 +45,7 @@ std::unique_ptr<Poly> F4Reducer::classicReduceSPoly
if (tracingLevel >= 2)
std::cerr << "F4Reducer: Reducing single S-pair.\n";
+
QuadMatrix qm;
{
F4MatrixBuilder builder(basis);
diff --git a/src/mathicgb/MonomialMap.hpp b/src/mathicgb/MonomialMap.hpp
index 7971587..cd43d16 100755
--- a/src/mathicgb/MonomialMap.hpp
+++ b/src/mathicgb/MonomialMap.hpp
@@ -1,47 +1,47 @@
-#ifndef MATHICGB_MONOMIAL_MAP_GUARD
-#define MATHICGB_MONOMIAL_MAP_GUARD
-
-#include "PolyRing.hpp"
-#include <limits>
-
-#if defined(MATHICGB_USE_STD_HASH) && defined(MATHICGB_USE_CUSTOM_HASH)
-#error Only select one kind of hash table
-#endif
-
-// set default hash table type if nothing has been specified
-#if !defined(MATHICGB_USE_STD_HASH) && !defined(MATHICGB_USE_CUSTOM_HASH)
-#define MATHICGB_USE_CUSTOM_HASH
-#endif
-
-#ifdef MATHICGB_USE_CUSTOM_HASH
+#ifndef MATHICGB_MONOMIAL_MAP_GUARD
+#define MATHICGB_MONOMIAL_MAP_GUARD
+
+#include "PolyRing.hpp"
+#include <limits>
+
+#if defined(MATHICGB_USE_STD_HASH) && defined(MATHICGB_USE_CUSTOM_HASH)
+#error Only select one kind of hash table
+#endif
+
+// set default hash table type if nothing has been specified
+#if !defined(MATHICGB_USE_STD_HASH) && !defined(MATHICGB_USE_CUSTOM_HASH)
+#define MATHICGB_USE_CUSTOM_HASH
+#endif
+
+#ifdef MATHICGB_USE_CUSTOM_HASH
#include <memtailor.h>
-#include <vector>
-#include <algorithm>
-#elif defined(MATHICGB_USE_STD_HASH)
-#include <unordered_map>
-#else
-#include <map>
-#endif
-
-/** A mapping from monomials to MapToType. The interface is the same as
- for STL maps. The interface is not complete. Add methods you need when
- you need it.
-*/
-template<class MapToType>
-class MonomialMap;
-
-namespace MonomialMapInternal {
- // The map type MapClass is defined here. This is here in
- // this namespace to avoid cluttering the class definition with
- // a lot of ifdef's.
-
-#ifdef MATHICGB_USE_CUSTOM_HASH
+#include <vector>
+#include <algorithm>
+#elif defined(MATHICGB_USE_STD_HASH)
+#include <unordered_map>
+#else
+#include <map>
+#endif
+
+/** A mapping from monomials to MapToType. The interface is the same as
+ for STL maps. The interface is not complete. Add methods you need when
+ you need it.
+*/
+template<class MapToType>
+class MonomialMap;
+
+namespace MonomialMapInternal {
+ // The map type MapClass is defined here. This is here in
+ // this namespace to avoid cluttering the class definition with
+ // a lot of ifdef's.
+
+#ifdef MATHICGB_USE_CUSTOM_HASH
template<class MTT>
class MapClass {
public:
typedef MapClass Map;
- typedef std::pair<const_monomial, MTT> value_type;
- typedef MTT mapped_type;
+ typedef std::pair<const_monomial, MTT> value_type;
+ typedef MTT mapped_type;
private:
typedef exponent HashValue;
@@ -51,25 +51,25 @@ namespace MonomialMapInternal {
exponent mono[1];
};
- public:
- MapClass(const PolyRing& ring):
- mRing(ring),
- mTable(),
- mNodeAlloc(sizeof(Node) - sizeof(exponent) + ring.maxMonomialByteSize())
- {
- mGrowWhenThisManyEntries = 0;
- mTableSize = mTable.size();
- mEntryCount = 0;
- growTable();
- }
-
- Map& map() {return *this;}
- const Map& map() const {return *this;}
- const PolyRing& ring() const {return mRing;}
-
- mapped_type* find(const_monomial mono) {
- const HashValue monoHash = mRing.monomialHashValue(mono);
- Node* node = entry(hashToIndex(monoHash));
+ public:
+ MapClass(const PolyRing& ring):
+ mRing(ring),
+ mTable(),
+ mNodeAlloc(sizeof(Node) - sizeof(exponent) + ring.maxMonomialByteSize())
+ {
+ mGrowWhenThisManyEntries = 0;
+ mTableSize = mTable.size();
+ mEntryCount = 0;
+ growTable();
+ }
+
+ Map& map() {return *this;}
+ const Map& map() const {return *this;}
+ const PolyRing& ring() const {return mRing;}
+
+ mapped_type* find(const_monomial mono) {
+ const HashValue monoHash = mRing.monomialHashValue(mono);
+ Node* node = entry(hashToIndex(monoHash));
for (; node != 0; node = node->next) {
if (monoHash == mRing.monomialHashValue(node->mono) &&
mRing.monomialEqualHintTrue(mono, node->mono)
@@ -78,11 +78,11 @@ namespace MonomialMapInternal {
}
}
return 0;
- }
-
- mapped_type* findProduct(const_monomial a, const_monomial b) {
- const HashValue abHash = mRing.monomialHashOfProduct(a, b);
- Node* node = entry(hashToIndex(abHash));
+ }
+
+ mapped_type* findProduct(const_monomial a, const_monomial b) {
+ const HashValue abHash = mRing.monomialHashOfProduct(a, b);
+ Node* node = entry(hashToIndex(abHash));
for (; node != 0; node = node->next) {
if (abHash == mRing.monomialHashValue(node->mono) &&
mRing.monomialIsProductOfHintTrue(a, b, node->mono)
@@ -91,29 +91,32 @@ namespace MonomialMapInternal {
}
}
return 0;
- }
-
- void insert(const value_type& value) {
- Node* node = static_cast<Node*>(mNodeAlloc.alloc());
- const size_t index = hashToIndex(mRing.monomialHashValue(value.first));
- ring().monomialCopy(value.first, Monomial(node->mono));
- new (&node->value) mapped_type(value.second);
- node->next = entry(index);
+ }
+
+ void insert(const value_type& value) {
+ Node* node = static_cast<Node*>(mNodeAlloc.alloc());
+ const size_t index = hashToIndex(mRing.monomialHashValue(value.first));
+ {
+ Monomial nodeTmp(node->mono);
+ ring().monomialCopy(value.first, nodeTmp);
+ }
+ new (&node->value) mapped_type(value.second);
+ node->next = entry(index);
mTable[index] = node;
++mEntryCount;
MATHICGB_ASSERT(mEntryCount <= mGrowWhenThisManyEntries);
if (mEntryCount == mGrowWhenThisManyEntries)
growTable();
- }
-
- void clear() {
- std::fill(mTable.begin(), mTable.end(), static_cast<Node*>(0));
- mNodeAlloc.freeAllBuffers();
- mEntryCount = 0;
- }
-
- size_t size() const {return mEntryCount;}
-
+ }
+
+ void clear() {
+ std::fill(mTable.begin(), mTable.end(), static_cast<Node*>(0));
+ mNodeAlloc.freeAllBuffers();
+ mEntryCount = 0;
+ }
+
+ size_t size() const {return mEntryCount;}
+
private:
size_t hashToIndex(HashValue hash) const {
const auto index = hash & mHashToIndexMask;
@@ -147,11 +150,11 @@ namespace MonomialMapInternal {
for (auto tableIt = mTable.begin(); tableIt != tableEnd; ++tableIt) {
Node* node = *tableIt;
while (node != 0) {
- const size_t index =
- mRing.monomialHashValue(node->mono) & newHashToIndexMask;
- MATHICGB_ASSERT(index < newTable.size());
- Node* const next = node->next;
- node->next = newTable[index];
+ const size_t index =
+ mRing.monomialHashValue(node->mono) & newHashToIndexMask;
+ MATHICGB_ASSERT(index < newTable.size());
+ Node* const next = node->next;
+ node->next = newTable[index];
newTable[index] = node;
node = next;
}
@@ -175,188 +178,188 @@ namespace MonomialMapInternal {
memt::BufferPool mNodeAlloc; // nodes are allocated from here.
};
-#elif defined(MATHICGB_USE_STD_HASH)
- class Hash {
- public:
- Hash(const PolyRing& ring): mRing(ring) {}
- size_t operator()(const_monomial m) const {
- return mRing.monomialHashValue(m);
- }
- const PolyRing& ring() const {return mRing;}
-
- private:
- const PolyRing& mRing;
- };
-
- class Equal {
- public:
- Equal(const PolyRing& ring): mRing(ring) {}
- size_t operator()(const_monomial a, const_monomial b) const {
- return mRing.monomialEQ(a, b);
- }
- const PolyRing& ring() const {return mRing;}
-
- private:
- const PolyRing& mRing;
- };
-
- // This has to be a template to support rebind().
- template<class T>
- class TemplateHashAllocator {
- public:
- typedef T value_type;
- typedef T* pointer;
- typedef T& reference;
- typedef const T* const_pointer;
- typedef const T& const_reference;
- typedef ::size_t size_type;
- typedef ::ptrdiff_t difference_type;
-
- template<class T2>
- struct rebind {
- typedef TemplateHashAllocator<T2> other;
- };
-
- TemplateHashAllocator() {}
- TemplateHashAllocator(memt::Arena& arena): mArena(arena) {}
- template<class X>
- TemplateHashAllocator(const TemplateHashAllocator<X>& a):
- mArena(a.arena()) {}
- TemplateHashAllocator(const TemplateHashAllocator<T>& a):
- mArena(a.arena()) {}
-
- pointer address(reference x) {return &x;}
- const_pointer address(const_reference x) const {return &x;}
- pointer allocate(size_type n, void* hint = 0) {
- return static_cast<pointer>(mArena.alloc(sizeof(T) * n));
- }
- void deallocate(pointer p, size_t n) {}
- size_type max_size() const {return std::numeric_limits<size_type>::max();}
- void construct(pointer p, const_reference val) {new (p) T(val);}
- void destroy(pointer p) {p->~T();}
- memt::Arena& arena() const {return mArena;}
-
- private:
- mutable memt::Arena& mArena;
- };
-
- template<class MTT>
- class MapClass {
- public:
- typedef TemplateHashAllocator<std::pair<const const_monomial, MTT>>
- HashAllocator;
- typedef std::unordered_map<const_monomial, MTT, Hash, Equal, HashAllocator>
- Map;
- typedef typename Map::iterator iterator;
- typedef typename Map::const_iterator const_iterator;
- typedef typename Map::value_type value_type;
-
- MapClass(const PolyRing& ring):
- mArena(),
- mMap(100, Hash(ring), Equal(ring), HashAllocator(mArena)),
- mTmp(ring.allocMonomial())
- {
- mMap.max_load_factor(0.3f);
- }
-
- ~MapClass() {
- ring().freeMonomial(mTmp);
- }
-
- Map& map() {return mMap;}
- const Map& map() const {return mMap;}
- const PolyRing& ring() const {return mMap.key_eq().ring();}
-
- value_type* findProduct(const_monomial a, const_monomial b) {
- ring().setGivenHash(mTmp, ring().monomialHashOfProduct(a, b));
- size_t bucket = mMap.bucket(mTmp);
- auto stop = mMap.end(bucket);
- for (auto it = mMap.begin(bucket); it != stop; ++it)
- if (ring().monomialIsProductOf(a, b, it->first))
- return &*it;
- return 0;
- }
-
- private:
- memt::Arena mArena;
- Map mMap;
- monomial mTmp;
- };
-#elif defined(MATHICGB_USE_CUSTOM_HASH)
-#else
- /// We need SOME ordering to make std::map work.
- class ArbitraryOrdering {
- public:
- ArbitraryOrdering(const PolyRing& ring): mRing(ring) {}
- bool operator()(const_monomial a, const_monomial b) const {
- return mRing.monomialLT(a, b);
- }
- const PolyRing& ring() const {return mRing;}
-
- private:
- const PolyRing& mRing;
- };
-
- template<class MTT>
- class MapClass {
- public:
- typedef std::map<const_monomial, MTT, ArbitraryOrdering> Map;
-
- MapClass(const PolyRing& ring): mMap(ArbitraryOrdering(ring)) {}
- Map& map() {return mMap;}
- const Map& map() const {return mMap;}
- const PolyRing& ring() const {return mMap.key_cmp().ring();}
-
- private:
- Map mMap;
- };
-#endif
-}
-
-template<class MTT>
-class MonomialMap {
-public:
- typedef MTT mapped_type;
- typedef MonomialMapInternal::MapClass<MTT> MapType;
-
- //typedef typename MapType::Map::iterator iterator;
- //typedef typename MapType::Map::const_iterator const_iterator;
- typedef typename MapType::Map::value_type value_type;
-
- MonomialMap(const PolyRing& ring): mMap(ring) {}
-
-/* iterator begin() {return map().begin();}
- const_iterator begin() const {return map().begin();}
- iterator end() {return map().end();}
- const_iterator end() const {return map().end();}*/
-
- mapped_type* find(const_monomial m) {return mMap.find(m);}
-
- mapped_type* findProduct(const_monomial a, const_monomial b) {
- return mMap.findProduct(a, b);
- }
-
- const mapped_type* find(const_monomial m) const {
- return const_cast<MonomialMap&>(*this).find(m);
- }
-
- const mapped_type* findProduct(const_monomial a, const_monomial b) const {
- return const_cast<MonomialMap&>(*this).findProduct(a, b);
- }
-
- size_t size() const {return map().size();}
- void clear() {map().clear();}
- void insert(const value_type& val) {
- map().insert(val);
- }
-
- inline const PolyRing& ring() const {return mMap.ring();}
-
-private:
- typename MapType::Map& map() {return mMap.map();}
- const typename MapType::Map& map() const {return mMap.map();}
-
- MapType mMap;
-};
-
-#endif
+#elif defined(MATHICGB_USE_STD_HASH)
+ class Hash {
+ public:
+ Hash(const PolyRing& ring): mRing(ring) {}
+ size_t operator()(const_monomial m) const {
+ return mRing.monomialHashValue(m);
+ }
+ const PolyRing& ring() const {return mRing;}
+
+ private:
+ const PolyRing& mRing;
+ };
+
+ class Equal {
+ public:
+ Equal(const PolyRing& ring): mRing(ring) {}
+ size_t operator()(const_monomial a, const_monomial b) const {
+ return mRing.monomialEQ(a, b);
+ }
+ const PolyRing& ring() const {return mRing;}
+
+ private:
+ const PolyRing& mRing;
+ };
+
+ // This has to be a template to support rebind().
+ template<class T>
+ class TemplateHashAllocator {
+ public:
+ typedef T value_type;
+ typedef T* pointer;
+ typedef T& reference;
+ typedef const T* const_pointer;
+ typedef const T& const_reference;
+ typedef ::size_t size_type;
+ typedef ::ptrdiff_t difference_type;
+
+ template<class T2>
+ struct rebind {
+ typedef TemplateHashAllocator<T2> other;
+ };
+
+ TemplateHashAllocator() {}
+ TemplateHashAllocator(memt::Arena& arena): mArena(arena) {}
+ template<class X>
+ TemplateHashAllocator(const TemplateHashAllocator<X>& a):
+ mArena(a.arena()) {}
+ TemplateHashAllocator(const TemplateHashAllocator<T>& a):
+ mArena(a.arena()) {}
+
+ pointer address(reference x) {return &x;}
+ const_pointer address(const_reference x) const {return &x;}
+ pointer allocate(size_type n, void* hint = 0) {
+ return static_cast<pointer>(mArena.alloc(sizeof(T) * n));
+ }
+ void deallocate(pointer p, size_t n) {}
+ size_type max_size() const {return std::numeric_limits<size_type>::max();}
+ void construct(pointer p, const_reference val) {new (p) T(val);}
+ void destroy(pointer p) {p->~T();}
+ memt::Arena& arena() const {return mArena;}
+
+ private:
+ mutable memt::Arena& mArena;
+ };
+
+ template<class MTT>
+ class MapClass {
+ public:
+ typedef TemplateHashAllocator<std::pair<const const_monomial, MTT>>
+ HashAllocator;
+ typedef std::unordered_map<const_monomial, MTT, Hash, Equal, HashAllocator>
+ Map;
+ typedef typename Map::iterator iterator;
+ typedef typename Map::const_iterator const_iterator;
+ typedef typename Map::value_type value_type;
+
+ MapClass(const PolyRing& ring):
+ mArena(),
+ mMap(100, Hash(ring), Equal(ring), HashAllocator(mArena)),
+ mTmp(ring.allocMonomial())
+ {
+ mMap.max_load_factor(0.3f);
+ }
+
+ ~MapClass() {
+ ring().freeMonomial(mTmp);
+ }
+
+ Map& map() {return mMap;}
+ const Map& map() const {return mMap;}
+ const PolyRing& ring() const {return mMap.key_eq().ring();}
+
+ value_type* findProduct(const_monomial a, const_monomial b) {
+ ring().setGivenHash(mTmp, ring().monomialHashOfProduct(a, b));
+ size_t bucket = mMap.bucket(mTmp);
+ auto stop = mMap.end(bucket);
+ for (auto it = mMap.begin(bucket); it != stop; ++it)
+ if (ring().monomialIsProductOf(a, b, it->first))
+ return &*it;
+ return 0;
+ }
+
+ private:
+ memt::Arena mArena;
+ Map mMap;
+ monomial mTmp;
+ };
+#elif defined(MATHICGB_USE_CUSTOM_HASH)
+#else
+ /// We need SOME ordering to make std::map work.
+ class ArbitraryOrdering {
+ public:
+ ArbitraryOrdering(const PolyRing& ring): mRing(ring) {}
+ bool operator()(const_monomial a, const_monomial b) const {
+ return mRing.monomialLT(a, b);
+ }
+ const PolyRing& ring() const {return mRing;}
+
+ private:
+ const PolyRing& mRing;
+ };
+
+ template<class MTT>
+ class MapClass {
+ public:
+ typedef std::map<const_monomial, MTT, ArbitraryOrdering> Map;
+
+ MapClass(const PolyRing& ring): mMap(ArbitraryOrdering(ring)) {}
+ Map& map() {return mMap;}
+ const Map& map() const {return mMap;}
+ const PolyRing& ring() const {return mMap.key_cmp().ring();}
+
+ private:
+ Map mMap;
+ };
+#endif
+}
+
+template<class MTT>
+class MonomialMap {
+public:
+ typedef MTT mapped_type;
+ typedef MonomialMapInternal::MapClass<MTT> MapType;
+
+ //typedef typename MapType::Map::iterator iterator;
+ //typedef typename MapType::Map::const_iterator const_iterator;
+ typedef typename MapType::Map::value_type value_type;
+
+ MonomialMap(const PolyRing& ring): mMap(ring) {}
+
+/* iterator begin() {return map().begin();}
+ const_iterator begin() const {return map().begin();}
+ iterator end() {return map().end();}
+ const_iterator end() const {return map().end();}*/
+
+ mapped_type* find(const_monomial m) {return mMap.find(m);}
+
+ mapped_type* findProduct(const_monomial a, const_monomial b) {
+ return mMap.findProduct(a, b);
+ }
+
+ const mapped_type* find(const_monomial m) const {
+ return const_cast<MonomialMap&>(*this).find(m);
+ }
+
+ const mapped_type* findProduct(const_monomial a, const_monomial b) const {
+ return const_cast<MonomialMap&>(*this).findProduct(a, b);
+ }
+
+ size_t size() const {return map().size();}
+ void clear() {map().clear();}
+ void insert(const value_type& val) {
+ map().insert(val);
+ }
+
+ inline const PolyRing& ring() const {return mMap.ring();}
+
+private:
+ typename MapType::Map& map() {return mMap.map();}
+ const typename MapType::Map& map() const {return mMap.map();}
+
+ MapType mMap;
+};
+
+#endif
diff --git a/src/mathicgb/QuadMatrix.cpp b/src/mathicgb/QuadMatrix.cpp
index b3d9ff5..f846f8f 100755
--- a/src/mathicgb/QuadMatrix.cpp
+++ b/src/mathicgb/QuadMatrix.cpp
@@ -62,7 +62,7 @@ std::string QuadMatrix::toString() const {
return out.str();
}
-std::ostream& operator<<(std::ostream& out, const QuadMatrix qm) {
+std::ostream& operator<<(std::ostream& out, const QuadMatrix& qm) {
qm.print(out);
return out;
}
diff --git a/src/mathicgb/QuadMatrix.hpp b/src/mathicgb/QuadMatrix.hpp
index 345c234..a8fb565 100755
--- a/src/mathicgb/QuadMatrix.hpp
+++ b/src/mathicgb/QuadMatrix.hpp
@@ -17,6 +17,8 @@ class ostream;
QuadMatrixBuilder or F4MatrixBuilder. */
class QuadMatrix {
public:
+ QuadMatrix() {}
+
SparseMatrix topLeft;
SparseMatrix topRight;
SparseMatrix bottomLeft;
@@ -35,8 +37,12 @@ public:
#ifdef MATHICGB_DEBUG
bool debugAssertValid() const;
#endif
+
+private:
+ QuadMatrix(const QuadMatrix&); // not available
+ void operator=(const QuadMatrix&); // not available
};
-std::ostream& operator<<(std::ostream& out, const QuadMatrix qm);
+std::ostream& operator<<(std::ostream& out, const QuadMatrix& qm);
#endif
diff --git a/src/mathicgb/RawVector.hpp b/src/mathicgb/RawVector.hpp
new file mode 100755
index 0000000..aef9bc7
--- /dev/null
+++ b/src/mathicgb/RawVector.hpp
@@ -0,0 +1,286 @@
+#ifndef MATHICGB_RAW_VECTOR_GUARD
+#define MATHICGB_RAW_VECTOR_GUARD
+
+#include <iterator>
+#include <cstddef>
+#include <memory>
+#include <utility>
+#include <algorithm>
+#include <stdexcept>
+
+/// RawVector mimics std::vector except that it does not do memory management.
+/// So you can directly access the pointers to memory and you can replace them
+/// with any other pointers that you need to.
+///
+/// Warning: RawVector is called raw because it does not allocate any memory,
+/// so it is an error to insert elements that there is not enough space for -
+/// you need to make sure ahead of time that there is enough space.
+///
+/// This class makes it possible to have a lot of the convenience of the
+/// std::vector interface even in places where std::vector cannot be used
+/// because of a need for manual memory management.
+template<class T>
+class RawVector {
+public:
+ typedef T& reference;
+ typedef const T& const_reference;
+ typedef T* iterator;
+ typedef const T* const_iterator;
+ typedef size_t size_type;
+ typedef ptrdiff_t difference_type;
+ typedef T value_type;
+ typedef T* pointer;
+ typedef const T* const_pointer;
+ typedef std::reverse_iterator<iterator> reverse_iterator;
+ typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
+
+ /// Initializes to the null state.
+ RawVector(): mBegin(0), mEnd(0), mCapacityEnd(0) {}
+
+ /// Copies the pointers from v. It is a shallow copy.
+ RawVector(const RawVector& v):
+ mBegin(v.mBegin()),
+ mEnd(v.mEnd()),
+ mCapacityEnd(v.mCapacityEnd) {}
+
+ /// Copies the pointers from v. It is a shallow copy. Sets v to a null state.
+ RawVector(RawVector&& v):
+ mBegin(v.mBegin),
+ mEnd(v.mEnd),
+ mCapacityEnd(v.mCapacityEnd)
+ {
+ v.releaseMemory();
+ }
+
+ RawVector(pointer begin, pointer end):
+ mBegin(begin),
+ mEnd(end),
+ mCapacityEnd(end) {}
+
+ RawVector(pointer begin, pointer end, pointer capacityEnd):
+ mBegin(begin),
+ mEnd(end),
+ mCapacityEnd(capacityEnd) {}
+
+ /// As clear() --- does not free memory.
+ ~RawVector() {
+ clear();
+ }
+
+ /// There is no operator== to avoid confusion about whether it compares
+ /// pointers or contents. Instead, there is this.
+ bool contentsEqual(const RawVector<T>& v) const {
+ return std::equal(begin(), end(), v.begin());
+ }
+
+ /// Copies the pointers from v. It is a shallow copy.
+ RawVector<T>& operator=(const RawVector& v) {
+ mBegin = v.mBegin;
+ mEnd = v.mEnd;
+ mCapacityEnd = v.mCapacityEnd;
+ return *this;
+ }
+
+ /// Copies the pointers from v. It is a shallow copy. Sets v to the null state.
+ RawVector<T>& operator=(RawVector&& v) {
+ mBegin = v.mBegin;
+ mEnd = v.mEnd;
+ mCapacityEnd = v.mCapacityEnd;
+ v.releaseMemory();
+ return *this;
+ }
+
+ iterator begin() {return mBegin;}
+ const_iterator begin() const {return mBegin;}
+ iterator end() {return mEnd;}
+ const_iterator end() const {return mEnd;}
+
+ reverse_iterator rbegin() {return reverse_iterator(end());}
+ const_reverse_iterator rbegin() const {return const_reverse_iterator(end());}
+ reverse_iterator rend() {return reverse_iterator(begin());}
+ const_reverse_iterator rend() const {return const_reverse_iterator(begin());}
+
+ size_type size() const {return mEnd - mBegin;}
+ size_type max_size() const {return std::numeric_limits<size_type>::max();}
+
+ /// There must be enough capacity for the new size.
+ size_type resize(const size_type newSize) {
+ MATHICGB_ASSERT(newSize <= capacity());
+ while (newSize < size()) {
+ new (end) T();
+ ++end;
+ }
+ while (newSize > size()) {
+ --end;
+ end->~T();
+ }
+ MATHICGB_ASSERT(newSize == size());
+ }
+
+ size_type capacity() const {return mCapacityEnd - mBegin;}
+
+ bool empty() const {return mBegin == mEnd;}
+
+ reference operator[](size_type index) {
+ MATHICGB_ASSERT(index < size());
+ return mBegin[index];
+ }
+
+ const_reference operator[](size_type index) const {
+ MATHICGB_ASSERT(index < size());
+ return mBegin[index];
+ }
+
+ reference at(size_type index) {
+ if (index >= size())
+ throw std::out_of_range("Invalid index in RawVector::at.");
+ return mBegin[index];
+ }
+
+ const_reference at(size_type index) const {
+ if (index >= size())
+ throw std::out_of_range("Invalid index in RawVector::at.");
+ return mBegin[index];
+ }
+
+ reference front() {
+ MATHICGB_ASSERT(!empty());
+ return *mBegin;
+ }
+
+ const_reference front() const {
+ MATHICGB_ASSERT(!empty());
+ return *mBegin;
+ }
+
+ reference back() {
+ MATHICGB_ASSERT(!empty());
+ return *(mEnd - 1);
+ }
+
+ const_reference back() const {
+ MATHICGB_ASSERT(!empty());
+ return *(mEnd - 1);
+ }
+
+ /// There must be enough capacity for the new size.
+ template<class Iter>
+ void rawAssign(Iter begin, Iter end) {
+ const auto count = std::distance(begin, end);
+ MATHICGB_ASSERT(count <= capacity());
+ if (count > size())
+ resize(count);
+ std::copy(begin, end, mBegin);
+ }
+
+ /// There must be enough capacity for the new size.
+ void rawAssign(size_type count, const T& t) {
+ MATHICGB_ASSERT(count <= capacity());
+ if (count > size())
+ resize(count);
+ std::fill_n(begin(), count, t);
+ }
+
+ /// There must be enough capacity for the new size. This is why there is no
+ /// push_back and only a rawPushBack --- to remind you of the changed
+ /// contract.
+ void rawPushBack(const T& t) {
+ MATHICGB_ASSERT(size() < capacity());
+ new (mEnd) T(t);
+ ++mEnd;
+ }
+
+ void pop_back() {
+ MATHICGB_ASSERT(!empty());
+ mEnd->~T();
+ --mEnd;
+ }
+
+ void swap(RawVector<T>& v) {
+ std::swap(mBegin, v.mBegin);
+ std::swap(mEnd, v.mEnd);
+ std::swap(mCapacityEnd, v.mCapacityEnd);
+ }
+
+ void clear() {
+ while (!empty())
+ pop_back();
+ }
+
+ // **** Extra functionality not available on std::vector
+
+ /// Returns true if there is no more capacity left. That is, if
+ /// capacity() == size().
+ bool atCapacity() const {
+ return mEnd == mCapacityEnd;
+ }
+
+ /// Sets this object to its null state without destructing memory. Returns
+ /// a pointer to the previous beginning of the buffer.
+ pointer releaseMemory() {
+ const auto oldBegin = mBegin;
+ mBegin = 0;
+ mEnd = 0;
+ mCapacityEnd = 0;
+ return oldBegin;
+ }
+
+ /// Takes over the new memory that is passed in, copies the old values to the
+ /// new memory and destructs the old values. Returns a pointer to the previous
+ /// beginning of the buffer.
+ pointer setMemoryAndCopy(
+ pointer begin,
+ pointer capacityEnd
+ ) {
+ MATHICGB_ASSERT(size() <=
+ static_cast<size_type>(std::distance(begin, capacityEnd)));
+ const auto oldBegin = mBegin;
+ const auto end = std::copy(mBegin, mEnd, begin);
+ *this = RawVector<T>(begin, end, capacityEnd);
+ return oldBegin;
+ }
+
+ /// Takes over the new memory that is passed in without destructing memory.
+ /// Returns a pointer to the previous beginning of the buffer.
+ pointer releaseAndSetMemory(
+ pointer begin,
+ pointer end,
+ pointer capacityEnd
+ ) {
+ const auto oldBegin = mBegin;
+ mBegin = begin;
+ mEnd = end;
+ mCapacityEnd = capacityEnd;
+ return oldBegin;
+ }
+
+ /// Destructs memory and then takes over the new memory that is passed in.
+ /// Does not deallocate the backing memory. Returns a pointer to the previous
+ /// beginning of the buffer.
+ pointer clearAndSetMemory(
+ pointer begin,
+ pointer end,
+ pointer capacityEnd
+ ) {
+ clear();
+ const auto oldBegin = mBegin;
+ mBegin = begin;
+ mEnd = end;
+ mCapacityEnd = capacityEnd;
+ return oldBegin;
+ }
+
+private:
+ T* mBegin;
+ T* mEnd;
+ T* mCapacityEnd;
+};
+
+namespace std {
+ template<class T>
+ void swap(RawVector<T>& a, RawVector<T>& b) {
+ a.swap(b);
+ }
+}
+
+#endif
diff --git a/src/mathicgb/SparseMatrix.cpp b/src/mathicgb/SparseMatrix.cpp
index 4ea5e1b..83c7f74 100755
--- a/src/mathicgb/SparseMatrix.cpp
+++ b/src/mathicgb/SparseMatrix.cpp
@@ -55,7 +55,7 @@ void SparseMatrix::sortRowsByIncreasingPivots() {
void SparseMatrix::applyColumnMap(std::vector<ColIndex> colMap) {
MATHICGB_ASSERT(colMap.size() >= colCount());
- auto end = mColIndices.end();
+ const auto end = mColIndices.end();
for (auto it = mColIndices.begin(); it != end; ++it) {
MATHICGB_ASSERT(*it < colCount());
*it = colMap[*it];
@@ -194,7 +194,7 @@ bool SparseMatrix::appendRowWithModulusIfNonZero(std::vector<uint64> const& v, S
bool SparseMatrix::operator==(const SparseMatrix& mat) const {
return mColCount == mat.mColCount &&
- mColIndices == mat.mColIndices &&
- mEntries == mat.mEntries &&
+ mColIndices.contentsEqual(mat.mColIndices) &&
+ mEntries.contentsEqual(mat.mEntries) &&
mRowOffsets == mat.mRowOffsets;
}
diff --git a/src/mathicgb/SparseMatrix.hpp b/src/mathicgb/SparseMatrix.hpp
index 759aedc..b470825 100755
--- a/src/mathicgb/SparseMatrix.hpp
+++ b/src/mathicgb/SparseMatrix.hpp
@@ -1,6 +1,7 @@
#ifndef MATHICGB_SPARSE_MATRIX_GUARD
#define MATHICGB_SPARSE_MATRIX_GUARD
+#include "RawVector.hpp"
#include "PolyRing.hpp"
#include <mathic.h>
#include <vector>
@@ -44,10 +45,39 @@ class SparseMatrix {
typedef uint32 ColIndex;
typedef uint16 Scalar;
+ SparseMatrix(SparseMatrix&& matrix):
+ mColIndices(std::move(matrix.mColIndices)),
+ mEntries(std::move(matrix.mEntries)),
+ mRowOffsets(std::move(matrix.mRowOffsets)),
+ mColCount(matrix.mColCount)
+ {
+ }
+
+ SparseMatrix& operator=(SparseMatrix&& matrix) {
+ this->~SparseMatrix();
+ new (this) SparseMatrix(std::move(matrix));
+ return *this;
+ }
+
+ ~SparseMatrix() {
+ delete[] mEntries.releaseMemory();
+ }
+
/** Preallocate space for at least count entries. */
void reserveEntries(size_t count) {
- mEntries.reserve(count);
- mColIndices.reserve(count);
+ if (count < mEntries.capacity())
+ return;
+
+ {
+ const auto begin = new Scalar[count];
+ const auto capacityEnd = begin + count;
+ delete[] mEntries.setMemoryAndCopy(begin, capacityEnd);
+ }
+ {
+ const auto begin = new ColIndex[count];
+ const auto capacityEnd = begin + count;
+ delete[] mColIndices.setMemoryAndCopy(begin, capacityEnd);
+ }
}
/** Preallocate space for at least count rows. */
@@ -77,10 +107,8 @@ class SparseMatrix {
trimThisMany, even if the scalar of that entry is set to zero. */
void trimLeadingZeroColumns(ColIndex trimThisMany) {
MATHICGB_ASSERT(trimThisMany <= colCount());
- typedef std::vector<ColIndex>::iterator Iter;
- Iter end = mColIndices.end();
- Iter it = mColIndices.begin();
- for (; it != end; ++it) {
+ const auto end = mColIndices.end();
+ for (auto it = mColIndices.begin(); it != end; ++it) {
MATHICGB_ASSERT(*it >= trimThisMany);
*it -= trimThisMany;
}
@@ -88,7 +116,9 @@ class SparseMatrix {
}
/** Construct a matrix with no rows and colCount columns. */
- SparseMatrix(ColIndex colCount = 0): mColCount(colCount) {
+ SparseMatrix(ColIndex colCount = 0):
+ mColCount(colCount)
+ {
mRowOffsets.push_back(0);
}
@@ -99,9 +129,7 @@ class SparseMatrix {
}
/** Returns the number of entries in the whole matrix. */
- size_t entryCount() const {
- return mEntries.size();
- }
+ size_t entryCount() const {return mEntries.size();}
/** Prints the matrix in a human readable format to out. */
void print(std::ostream& out) const;
@@ -111,7 +139,7 @@ class SparseMatrix {
/** Adds a new row that contains all terms that have been appended
since the last time a row was added or the matrix was created. */
void rowDone() {
- MATHICGB_ASSERT(mColIndices.size() == mEntries.size());
+ MATHICGB_ASSERT(mColIndices.size() == entryCount());
mRowOffsets.push_back(mColIndices.size());
}
@@ -119,13 +147,19 @@ class SparseMatrix {
until rowDone is called. Do not call other methods that add rows
after calling this method until rowDone has been called. */
void appendEntry(ColIndex colIndex, Scalar scalar) {
- MATHICGB_ASSERT(mColIndices.size() == mEntries.size());
+ MATHICGB_ASSERT(mColIndices.size() == entryCount());
MATHICGB_ASSERT(colIndex < colCount());
- mColIndices.push_back(colIndex);
- mEntries.push_back(scalar);
+ MATHICGB_ASSERT(mEntries.atCapacity() == mColIndices.atCapacity());
+ if (mEntries.atCapacity())
+ growEntryCapacity();
+ MATHICGB_ASSERT(!mEntries.atCapacity());
+ MATHICGB_ASSERT(!mColIndices.atCapacity());
- MATHICGB_ASSERT(mColIndices.size() == mEntries.size());
+ mColIndices.rawPushBack(colIndex);
+ mEntries.rawPushBack(scalar);
+
+ MATHICGB_ASSERT(mColIndices.size() == entryCount());
}
void appendRowAndNormalize(const SparseMatrix& matrix, RowIndex row, Scalar modulus);
@@ -212,18 +246,34 @@ class SparseMatrix {
private:
friend class RowIterator;
- ColIndex indexAtOffset(size_t offset) const {
- MATHICGB_ASSERT(offset < mColIndices.size());
- return mColIndices[offset];
- }
+ SparseMatrix(const SparseMatrix&); // not available
+ void operator=(const SparseMatrix&); // not available
- Scalar scalarAtOffset(size_t offset) const {
- MATHICGB_ASSERT(offset < mEntries.size());
- return mEntries[offset];
+ ColIndex indexAtOffset(size_t offset) const {return mColIndices[offset];}
+ Scalar scalarAtOffset(size_t offset) const {return mEntries[offset];}
+
+ void growEntryCapacity() {
+ MATHICGB_ASSERT(mColIndices.size() == mEntries.size());
+ MATHICGB_ASSERT(mColIndices.capacity() == mEntries.capacity());
+
+ const size_t initialCapacity = 1 << 16;
+ const size_t growthFactor = 2;
+ const size_t newCapacity =
+ mEntries.empty() ? initialCapacity : mEntries.capacity() * growthFactor;
+ reserveEntries(newCapacity);
+
+ MATHICGB_ASSERT(mColIndices.size() == mEntries.size());
+ MATHICGB_ASSERT(mColIndices.capacity() == newCapacity);
+ MATHICGB_ASSERT(mEntries.capacity() == newCapacity);
}
- std::vector<ColIndex> mColIndices;
- std::vector<Scalar> mEntries;
+ /// We need a RawVector here to tie the checks for the need to reallocate
+ /// together between mColIndices and mEntries. We only need to check
+ /// the capacity once, which, believe it or not, is a significant performance
+ /// win. Not least because it decreases the amount of code and therefore
+ /// causes different compiler inlining decisions.
+ RawVector<Scalar> mEntries;
+ RawVector<ColIndex> mColIndices;
std::vector<size_t> mRowOffsets;
ColIndex mColCount;
};
--
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