[mathicgb] 90/393: Got a 2% speed-up by processing two hash table monomial product queries at the same time.

Doug Torrance dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:38 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 cb11a5b2f8059a8458ef9540be5c3a8790c1a3b5
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date:   Thu Nov 1 16:10:33 2012 +0100

    Got a 2% speed-up by processing two hash table monomial product queries at the same time.
---
 src/mathicgb/F4MatrixBuilder.cpp   | 56 ++++++++++++++++++++++++++++++++++++--
 src/mathicgb/F4MatrixBuilder.hpp   |  8 ++++++
 src/mathicgb/MonomialMap.hpp       | 46 ++++++++++++++++++++++++++++++-
 src/mathicgb/Poly.hpp              | 17 ++++++++++--
 src/mathicgb/PolyRing.hpp          | 37 +++++++++++++++++++++++--
 src/mathicgb/QuadMatrixBuilder.hpp | 24 +++++++++++++---
 6 files changed, 176 insertions(+), 12 deletions(-)

diff --git a/src/mathicgb/F4MatrixBuilder.cpp b/src/mathicgb/F4MatrixBuilder.cpp
index 9df7d74..99315f3 100755
--- a/src/mathicgb/F4MatrixBuilder.cpp
+++ b/src/mathicgb/F4MatrixBuilder.cpp
@@ -5,7 +5,7 @@
 #include <omp.h>
 #endif
 
-MATHIC_INLINE QuadMatrixBuilder::LeftRightColIndex
+MATHICGB_INLINE QuadMatrixBuilder::LeftRightColIndex
   F4MatrixBuilder::findOrCreateColumn
 (
   const const_monomial monoA,
@@ -20,6 +20,35 @@ MATHIC_INLINE QuadMatrixBuilder::LeftRightColIndex
   return createColumn(builder, monoA, monoB);
 }
 
+MATHICGB_INLINE const std::pair<
+  QuadMatrixBuilder::LeftRightColIndex,
+  QuadMatrixBuilder::LeftRightColIndex
+> F4MatrixBuilder::findOrCreateTwoColumns
+(
+  const const_monomial monoA1,
+  const const_monomial monoA2,
+  const const_monomial monoB,
+  QuadMatrixBuilder& builder
+) {
+  MATHICGB_ASSERT(!monoA1.isNull());
+  MATHICGB_ASSERT(!monoA2.isNull());
+  MATHICGB_ASSERT(!monoB.isNull());
+  MATHICGB_ASSERT(!ring().monomialEQ(monoA1, monoA2));
+  auto colPair = builder.findTwoColumnsProduct(monoA1, monoA2, monoB);
+  if (!colPair.first.valid()) {
+    colPair.first = createColumn(builder, monoA1, monoB);
+    // Syncing builder to mBuilder could have created a col for monoA2*monoB.
+    if (&mBuilder != &builder && !colPair.second.valid())
+      colPair.second = builder.findColumnProduct(monoA2, monoB);
+  }
+  if (!colPair.second.valid())
+    colPair.second = createColumn(builder, monoA2, monoB);
+  MATHICGB_ASSERT(colPair == std::make_pair(
+    findOrCreateColumn(monoA1, monoB, builder),
+    findOrCreateColumn(monoA2, monoB, builder)));
+  return colPair;
+}
+
 F4MatrixBuilder::F4MatrixBuilder(
   const PolyBasis& basis,
   const int threadCount,
@@ -327,12 +356,33 @@ void F4MatrixBuilder::appendRowTop(
   MATHICGB_ASSERT(&poly != 0);
   MATHICGB_ASSERT(&builder != 0);
 
-  const Poly::const_iterator end = poly.end();
-  for (Poly::const_iterator it = poly.begin(); it != end; ++it) {
+  auto it = poly.begin();
+  const auto end = poly.end();
+  if ((std::distance(it, end) % 2) == 1) {
     const auto col = findOrCreateColumn(it.getMonomial(), multiple, builder);
 	MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
     MATHICGB_ASSERT(it.getCoefficient());
     builder.appendEntryTop(col, static_cast<Scalar>(it.getCoefficient()));
+    ++it;
+  }
+  MATHICGB_ASSERT((std::distance(it, end) % 2) == 0);
+  while (it != end) {
+	MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
+    MATHICGB_ASSERT(it.getCoefficient() != 0);
+    const auto scalar1 = static_cast<Scalar>(it.getCoefficient());
+    const const_monomial mono1 = it.getMonomial();
+    ++it;
+
+	MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
+    MATHICGB_ASSERT(it.getCoefficient() != 0);
+    const auto scalar2 = static_cast<Scalar>(it.getCoefficient());
+    const const_monomial mono2 = it.getMonomial();
+    ++it;
+
+    const auto colPair =
+      findOrCreateTwoColumns(mono1, mono2, multiple, builder);
+    builder.appendEntryTop(colPair.first, scalar1);
+    builder.appendEntryTop(colPair.second, scalar2);
   }
   builder.rowDoneTopLeftAndRight();
 }
diff --git a/src/mathicgb/F4MatrixBuilder.hpp b/src/mathicgb/F4MatrixBuilder.hpp
index 2d98368..51f797f 100755
--- a/src/mathicgb/F4MatrixBuilder.hpp
+++ b/src/mathicgb/F4MatrixBuilder.hpp
@@ -100,6 +100,14 @@ private:
   MATHIC_INLINE LeftRightColIndex findOrCreateColumn
     (const_monomial monoA, const_monomial monoB, QuadMatrixBuilder& builder);
 
+  MATHIC_INLINE const std::pair<LeftRightColIndex, LeftRightColIndex>
+  findOrCreateTwoColumns(
+    const const_monomial monoA1,
+    const const_monomial monoA2,
+    const const_monomial monoB,
+    QuadMatrixBuilder& builder
+  );
+
 
   const int mThreadCount;
   monomial mTmp;
diff --git a/src/mathicgb/MonomialMap.hpp b/src/mathicgb/MonomialMap.hpp
index f073e0f..f3c399e 100755
--- a/src/mathicgb/MonomialMap.hpp
+++ b/src/mathicgb/MonomialMap.hpp
@@ -80,7 +80,10 @@ namespace MonomialMapInternal {
       return 0;
     }
 
-    mapped_type* findProduct(const_monomial a, const_monomial b) {
+    MATHICGB_INLINE mapped_type* findProduct(
+      const const_monomial a,
+      const const_monomial b
+    ) {
       const HashValue abHash = mRing.monomialHashOfProduct(a, b);
       Node* node = entry(hashToIndex(abHash));
       for (; node != 0; node = node->next) {
@@ -94,6 +97,25 @@ namespace MonomialMapInternal {
       return 0;
     }
 
+    /// As findProduct but looks for a1*b and a2*b at one time.
+    MATHICGB_INLINE std::pair<mapped_type*, mapped_type*> findTwoProducts(
+      const const_monomial a1,
+      const const_monomial a2,
+      const const_monomial b
+    ) {
+      const HashValue a1bHash = mRing.monomialHashOfProduct(a1, b);
+      const HashValue a2bHash = mRing.monomialHashOfProduct(a2, b);
+      Node* const node1 = entry(hashToIndex(a1bHash));
+      Node* const node2 = entry(hashToIndex(a2bHash));
+
+      if (node1 != 0 && node2 != 0 && mRing.monomialIsTwoProductsOfHintTrue
+        (a1, a2, b, node1->mono, node2->mono)
+      )
+        return std::make_pair(&node1->value, &node2->value);
+      else
+        return std::make_pair(findProduct(a1, b), findProduct(a2, b));
+    }
+
     void insert(const value_type& value) {
       Node* node = static_cast<Node*>(mNodeAlloc.alloc());
       const size_t index = hashToIndex(mRing.monomialHashValue(value.first));
@@ -130,6 +152,11 @@ namespace MonomialMapInternal {
       return mTable[index];
     }
 
+    const Node* entry(size_t index) const {
+      MATHICGB_ASSERT(index < mTable.size());
+      return mTable[index];
+    }
+
     void growTable() {
       // Determine parameters for larger hash table
       const size_t initialTableSize = 1 << 16; // must be a power of two!!!
@@ -340,6 +367,14 @@ public:
     return mMap.findProduct(a, b);
   }
 
+  MATHICGB_INLINE std::pair<mapped_type*, mapped_type*> findTwoProducts(
+    const const_monomial a1,
+    const const_monomial a2,
+    const const_monomial b
+  ) {
+    return mMap.findTwoProducts(a1, a2, b);
+  }
+
   const mapped_type* find(const_monomial m) const {
     return const_cast<MonomialMap&>(*this).find(m);
   }
@@ -348,6 +383,15 @@ public:
     return const_cast<MonomialMap&>(*this).findProduct(a, b);
   }
 
+  MATHICGB_INLINE const std::pair<const mapped_type*, const mapped_type*>
+  findTwoProducts(
+    const const_monomial a1,
+    const const_monomial a2,
+    const const_monomial b
+  ) const {
+    return const_cast<MonomialMap&>(*this).findTwoProducts(a1, a2, b);
+  }
+
   size_t size() const {return map().size();}
   void clear() {map().clear();}
   void insert(const value_type& val) {
diff --git a/src/mathicgb/Poly.hpp b/src/mathicgb/Poly.hpp
index d3be160..61a2c5b 100755
--- a/src/mathicgb/Poly.hpp
+++ b/src/mathicgb/Poly.hpp
@@ -8,6 +8,7 @@
 #include <ostream>
 #include <utility>
 #include <cstdio>
+#include <iterator>
 
 class Poly {
 public:
@@ -29,6 +30,12 @@ public:
     iterator(Poly& f) : monsize(f.getRing()->maxMonomialSize()), ic(f.coeffs.begin()), im(f.monoms.begin()) {}
     iterator(Poly& f,int) : ic(f.coeffs.end()), im() {}
   public:
+    typedef std::random_access_iterator_tag iterator_category;
+    typedef std::pair<coefficient, const const_monomial> value_type;
+    typedef ptrdiff_t difference_type;
+    typedef value_type* pointer; // todo: is this OK?
+    typedef std::pair<coefficient&, const const_monomial> reference;
+
     iterator() {}
     iterator operator++() { ++ic; im += monsize; return *this; }
     coefficient &getCoefficient() const { return *ic; }
@@ -36,7 +43,7 @@ public:
     size_t operator-(const iterator &b) const { return ic - b.ic; }
     friend bool operator==(const iterator &a, const iterator &b);
     friend bool operator!=(const iterator &a, const iterator &b);
-    std::pair<coefficient&, monomial> operator*() const {
+    reference operator*() const {
       return std::pair<coefficient&, monomial>(getCoefficient(), getMonomial());
     }
   };
@@ -51,6 +58,12 @@ public:
     const_iterator(const Poly& f) : monsize(f.getRing()->maxMonomialSize()), ic(f.coeffs.begin()), im(f.monoms.begin()) {}
     const_iterator(const Poly& f,int) : ic(f.coeffs.end()), im() {}
   public:
+    typedef std::random_access_iterator_tag iterator_category;
+    typedef std::pair<const coefficient, const const_monomial> value_type;
+    typedef ptrdiff_t difference_type;
+    typedef value_type* pointer; // todo: is this OK?
+    typedef std::pair<const coefficient&, const const_monomial> reference;
+
     const_iterator() {}
     const_iterator operator++() { ++ic; im += monsize; return *this; }
     coefficient getCoefficient() const { return *ic; }
@@ -58,7 +71,7 @@ public:
     size_t operator-(const const_iterator &b) const { return ic - b.ic; }
     friend bool operator==(const const_iterator &a, const const_iterator &b);
     friend bool operator!=(const const_iterator &a, const const_iterator &b);
-    std::pair<coefficient, const_monomial> operator*() const {
+    const value_type operator*() const {
       return std::pair<coefficient, const_monomial>
         (getCoefficient(), getMonomial());
     }
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index 71668e7..0f357c4 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -416,11 +416,20 @@ public:
   bool monomialIsProductOf
     (ConstMonomial a, ConstMonomial b, ConstMonomial ab) const;
 
-  /**As monomialIsProductOf but optimized for the case that the result
-    is true. */
+  /// As monomialIsProductOf but optimized for the case that the result
+  /// is true.
   bool monomialIsProductOfHintTrue
     (ConstMonomial a, ConstMonomial b, ConstMonomial ab) const;
 
+  /// As monomialIsProductOfHintTwo(), but checks two products are equal.
+  /// The return value is true if a1*b = a1b and a2*b = a2b.
+  MATHICGB_INLINE bool monomialIsTwoProductsOfHintTrue(
+    ConstMonomial a1,
+    ConstMonomial a2,
+    ConstMonomial b,
+    ConstMonomial a1b,
+    ConstMonomial a2b) const;
+
   /// Returns the hash of the product of a and b.
   HashValue monomialHashOfProduct(ConstMonomial a, ConstMonomial b) const {
     return static_cast<HashValue>(a[mHashIndex]) +
@@ -672,6 +681,30 @@ inline bool PolyRing::monomialIsProductOfHintTrue(
   return orOfXor == 0; 
 }
 
+MATHICGB_INLINE bool PolyRing::monomialIsTwoProductsOfHintTrue(
+  const ConstMonomial a1,
+  const ConstMonomial a2,
+  const ConstMonomial b,
+  const ConstMonomial a1b,
+  const ConstMonomial a2b
+) const {
+  uint64 orOfXor = 0;
+  for (size_t i = mNumVars / 2; i != static_cast<size_t>(-1); --i) {
+    uint64 A1, A2, B, A1B, A2B;
+    std::memcpy(&A1, &a1[i*2], 8);
+    std::memcpy(&A2, &a2[i*2], 8);
+    std::memcpy(&B, &b[i*2], 8);
+    std::memcpy(&A1B, &a1b[i*2], 8);
+    std::memcpy(&A2B, &a2b[i*2], 8);
+    orOfXor |= (A1B ^ (A1 + B)) | (A2B ^ (A2 + B));
+  }
+  MATHICGB_ASSERT((orOfXor == 0) ==
+    monomialIsProductOf(a1, b, a1b) &&
+    monomialIsProductOf(a2, b, a2b));
+
+  return orOfXor == 0;
+}
+
 inline bool PolyRing::monomialIsProductOf(
   ConstMonomial a, 
   ConstMonomial b, 
diff --git a/src/mathicgb/QuadMatrixBuilder.hpp b/src/mathicgb/QuadMatrixBuilder.hpp
index cb0ca74..7b050d8 100755
--- a/src/mathicgb/QuadMatrixBuilder.hpp
+++ b/src/mathicgb/QuadMatrixBuilder.hpp
@@ -194,16 +194,32 @@ class QuadMatrixBuilder {
       return LeftRightColIndex();
   }
 
-  /// As findColumn, but looks for the monomial that is the product of a and b.
-  LeftRightColIndex findColumnProduct(const_monomial a, const_monomial b) const
-  {
-    auto it = mMonomialToCol.findProduct(a, b);
+  /// As findColumn(), but looks for a*b. This is faster than computing a*b
+  /// and then looking that up.
+  LeftRightColIndex findColumnProduct(
+    const const_monomial a,
+    const const_monomial b
+  ) const {
+    const auto it = mMonomialToCol.findProduct(a, b);
     if (it != 0)
       return *it;
     else
       return LeftRightColIndex();
   }
 
+  /// As findColumnProduct(), but looks for a1*b and a2*b at the same time.
+  MATHICGB_INLINE std::pair<LeftRightColIndex, LeftRightColIndex>
+  findTwoColumnsProduct(
+    const const_monomial a1,
+    const const_monomial a2,
+    const const_monomial b
+  ) const {
+    const auto it = mMonomialToCol.findTwoProducts(a1, a2, b);
+    return std::make_pair(
+      it.first != 0 ? *it.first : LeftRightColIndex(),
+      it.second != 0 ? *it.second : LeftRightColIndex());
+  }
+
   const_monomial monomialOfLeftCol(ColIndex col) const {
     MATHICGB_ASSERT(col < mMonomialsLeft.size());
     return mMonomialsLeft[col];

-- 
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