@@ -20,14 +20,14 @@ TEST_CASE("Test SparseIterator.h - One Dimensional")

         v.insert(7, 8.f);

         v.insert(9, 10.f);

-        SparseIterator it(v);

-        REQUIRE(it.getValue() == 1.f);

+        TemplatedSparseIterator<1> it(v);

+        REQUIRE(get<1>(it) == 1.f);

         it.next();

-        REQUIRE(it.getValue() == 5.f);

+        REQUIRE(get<1>(it) == 5.f);

         it.next();

-        REQUIRE(it.getValue() == 8.f);

+        REQUIRE(get<1>(it) == 8.f);

         it.next();

-        REQUIRE(it.getValue() == 10.f);

+        REQUIRE(get<1>(it) == 10.f);

         it.next();

         REQUIRE(it.atEnd());

     }

@@ -51,10 +51,10 @@ TEST_CASE("Test SparseIterator.h - One Dimensional")

         for (unsigned j = 0; j < ref.nCol(); ++j)

         {

             float colSum = 0.f;

-            SparseIterator it(mat.getCol(j));

+            TemplatedSparseIterator<1> it(mat.getCol(j));

             while (!it.atEnd())

             {

-                colSum += it.getValue();

+                colSum += get<1>(it);

                 it.next();

             }

             REQUIRE(colSum == gaps::sum(ref.getCol(j)));

@@ -79,12 +79,12 @@ TEST_CASE("Test SparseIterator.h - Two Dimensional")

         hv.add(6, 5.f);

         hv.add(7, 6.f);

-        SparseIteratorTwo it(sv, hv);

-        REQUIRE(it.getValue_1() == 5.f);

-        REQUIRE(it.getValue_2() == 3.f);

+        TemplatedSparseIterator<2> it(sv, hv);

+        REQUIRE(get<1>(it) == 5.f);

+        REQUIRE(get<2>(it) == 3.f);

         it.next();

-        REQUIRE(it.getValue_1() == 8.f);

-        REQUIRE(it.getValue_2() == 6.f);

+        REQUIRE(get<1>(it) == 8.f);

+        REQUIRE(get<2>(it) == 6.f);

         it.next();

         REQUIRE(it.atEnd());

     }

@@ -107,9 +107,9 @@ TEST_CASE("Test SparseIterator.h - Two Dimensional")

         hv.add(8, 9.f);

         hv.add(75, 76.f);

-        SparseIteratorTwo it(sv, hv);

-        REQUIRE(it.getValue_1() == 75.f);

-        REQUIRE(it.getValue_2() == 76.f);

+        TemplatedSparseIterator<2> it(sv, hv);

+        REQUIRE(get<1>(it) == 75.f);

+        REQUIRE(get<2>(it) == 76.f);

         it.next();

         REQUIRE(it.atEnd());

     }

@@ -170,7 +170,7 @@ TEST_CASE("Test SparseIterator.h - Two Dimensional")

         // calculate dot product

         float sdot = 0.f, ddot = 0.f;

-        SparseIteratorTwo it(sv, hv);

+        TemplatedSparseIterator<2> it(sv, hv);

         unsigned i = 0;

         while (!it.atEnd())

         {

@@ -182,12 +182,12 @@ TEST_CASE("Test SparseIterator.h - Two Dimensional")

             if (i < dv1.size())

             {

                 ddot += dv1[i] * dv2[i];

-                REQUIRE(dv1[i] == it.getValue_1());

-                REQUIRE(dv2[i] == it.getValue_2());

+                REQUIRE(dv1[i] == get<1>(it));

+                REQUIRE(dv2[i] == get<2>(it));

                 ++i;

             }

-            sdot += it.getValue_1() * it.getValue_2();

+            sdot += get<1>(it) * get<2>(it);

             it.next();

         }

@@ -218,10 +218,10 @@ TEST_CASE("Test SparseIterator.h - Two Dimensional")

             for (unsigned j2 = j1; j2 < ref.nCol(); ++j2)

             {

                 float dot = 0.f;

-                SparseIteratorTwo it(sMat.getCol(j1), hMat.getCol(j2));

+                TemplatedSparseIterator<2> it(sMat.getCol(j1), hMat.getCol(j2));

                 while (!it.atEnd())

                 {

-                    dot += it.getValue_1() * it.getValue_2();

+                    dot += get<1>(it) * get<2>(it);

                     it.next();

                 }

                 REQUIRE(dot == gaps::dot(ref.getCol(j1), ref.getCol(j2)));

@@ -265,14 +265,14 @@ TEST_CASE("Test SparseIterator.h - Three Dimensional")

         hv2.add(8, 7.f);

         hv2.add(9, 8.f);

-        SparseIteratorThree it(sv, hv1, hv2);

-        REQUIRE(it.getValue_1() == 5.f); // 4

-        REQUIRE(it.getValue_2() == 3.f);

-        REQUIRE(it.getValue_3() == 6.f);

+        TemplatedSparseIterator<3> it(sv, hv1, hv2);

+        REQUIRE(get<1>(it) == 5.f); // 4

+        REQUIRE(get<2>(it) == 3.f);

+        REQUIRE(get<3>(it) == 6.f);

         it.next();

-        REQUIRE(it.getValue_1() == 10.f); // 9

-        REQUIRE(it.getValue_2() == 7.f);

-        REQUIRE(it.getValue_3() == 8.f);

+        REQUIRE(get<1>(it) == 10.f); // 9

+        REQUIRE(get<2>(it) == 7.f);

+        REQUIRE(get<3>(it) == 8.f);

         it.next();

         REQUIRE(it.atEnd());

     }

@@ -301,11 +301,11 @@ TEST_CASE("Test SparseIterator.h - Three Dimensional")

                 for (unsigned j3 = j2; j3 < ref.nCol(); ++j3)

                 {

                     float prod = 0.f;

-                    SparseIteratorThree it(sMat.getCol(j1), hMat.getCol(j2),

+                    TemplatedSparseIterator<3> it(sMat.getCol(j1), hMat.getCol(j2),

                         hMat.getCol(j3));

                     while (!it.atEnd())

                     {

-                        prod += it.getValue_1() * it.getValue_2() * it.getValue_3();

+                        prod += get<1>(it) * get<2>(it) * get<3>(it);

                         it.next();

                     }

                     REQUIRE(prod == tripleProduct(ref.getCol(j1),

@@ -21,6 +21,9 @@ public:

     friend class SparseIteratorTwo;

     friend class SparseIteratorThree;

+    

+    template <unsigned N>

+    friend class TemplatedSparseIterator;

     explicit HybridVector(unsigned size);

     explicit HybridVector(const std::vector<float> &v);

@@ -18,6 +18,183 @@ static uint64_t clearLowerBits(uint64_t u, unsigned pos)

     return u & ~((1ull << (pos + 1ull)) - 1ull);

 }

+template <class Iter>

+void gotoNextCommon(Iter &it)

+{

+    // get the common indices in this chunk

+    it.calculateCommonFlags();

+

+    // if nothing common in this chunk, find a chunk that has common indices

+    while (!it.mCommonFlags)

+    {

+        // first count how many sparse indices we are skipping

+        it.mSparseIndex += __builtin_popcountll(it.mSparseFlags);

+        

+        // advance to next chunk

+        if (++it.mBigIndex == it.mTotalIndices)

+        {   

+            it.mAtEnd = true;

+            return;

+        }

+

+        // update the flags

+        it.getFlags();

+        it.calculateCommonFlags();

+    }

+

+    // must have at least one common value, this is our index

+    it.mSmallIndex = __builtin_ffsll(it.mCommonFlags) - 1;

+

+    // find the number of skipped entries in the sparse vector

+    it.mSparseIndex += 1 + countLowerBits(it.mSparseFlags, it.mSmallIndex);

+

+    // clear out all skipped indices and the current index from the bitflags

+    it.mSparseFlags = clearLowerBits(it.mSparseFlags, it.mSmallIndex);

+}

+

+template<>

+float get<1>(const TemplatedSparseIterator<1> &it)

+{

+    return it.mSparse.mData[it.mSparseIndex];

+}

+

+template<>

+float get<1>(const TemplatedSparseIterator<2> &it)

+{

+    return it.mSparse.mData[it.mSparseIndex];

+}

+

+template<>

+float get<1>(const TemplatedSparseIterator<3> &it)

+{

+    return it.mSparse.mData[it.mSparseIndex];

+}

+

+template<>

+float get<2>(const TemplatedSparseIterator<2> &it)

+{

+    return it.mHybrid_1[64 * it.mBigIndex + it.mSmallIndex];

+}

+

+template<>

+float get<2>(const TemplatedSparseIterator<3> &it)

+{

+    return it.mHybrid_1[64 * it.mBigIndex + it.mSmallIndex];

+}

+

+template<>

+float get<3>(const TemplatedSparseIterator<3> &it)

+{

+    return it.mHybrid_2[64 * it.mBigIndex + it.mSmallIndex];

+}

+

+

+TemplatedSparseIterator<1>::TemplatedSparseIterator(const SparseVector &v)

+:

+mSparse(v),

+mSparseIndex(0)

+{}

+

+bool TemplatedSparseIterator<1>::atEnd() const

+{

+    return mSparseIndex == mSparse.mData.size();

+}

+

+void TemplatedSparseIterator<1>::next()

+{

+    ++mSparseIndex;

+}

+

+TemplatedSparseIterator<2>::TemplatedSparseIterator(const SparseVector &v, const HybridVector &h)

+    :

+mSparse(v),

+mHybrid_1(h),

+mSparseFlags(v.mIndexBitFlags[0]),

+mHybridFlags_1(h.mIndexBitFlags[0]),

+mCommonFlags(v.mIndexBitFlags[0] & h.mIndexBitFlags[0]),

+mTotalIndices(v.mIndexBitFlags.size()),

+mBigIndex(0),

+mSmallIndex(0),

+mSparseIndex(0),

+mAtEnd(false)

+{

+    GAPS_ASSERT(v.size() == h.size());

+

+    next();

+    mSparseIndex -= 1; // next puts us at position 1, this resets to 0

+}

+

+bool TemplatedSparseIterator<2>::atEnd() const

+{

+    return mAtEnd;

+}

+

+void TemplatedSparseIterator<2>::next()

+{

+    gotoNextCommon(*this);

+}

+

+void TemplatedSparseIterator<2>::calculateCommonFlags()

+{

+    mCommonFlags = mSparseFlags & mHybridFlags_1;

+}

+

+void TemplatedSparseIterator<2>::getFlags()

+{

+    mSparseFlags = mSparse.mIndexBitFlags[mBigIndex];

+    mHybridFlags_1 = mHybrid_1.mIndexBitFlags[mBigIndex];

+}

+

+TemplatedSparseIterator<3>::TemplatedSparseIterator(const SparseVector &v,

+const HybridVector &h1, const HybridVector &h2)

+    :

+mSparse(v),

+mHybrid_1(h1),

+mHybrid_2(h2),

+mSparseFlags(v.mIndexBitFlags[0]),

+mHybridFlags_1(h1.mIndexBitFlags[0]),

+mHybridFlags_2(h2.mIndexBitFlags[0]),

+mCommonFlags(v.mIndexBitFlags[0] & h1.mIndexBitFlags[0] & h2.mIndexBitFlags[0]),

+mTotalIndices(v.mIndexBitFlags.size()),

+mBigIndex(0),

+mSmallIndex(0),

+mSparseIndex(0),

+mAtEnd(false)

+{

+    GAPS_ASSERT(v.size() == h1.size());

+    GAPS_ASSERT(h1.size() == h2.size());

+

+    next();

+    mSparseIndex -= 1;

+}

+

+bool TemplatedSparseIterator<3>::atEnd() const

+{

+    return mAtEnd;

+}

+

+void TemplatedSparseIterator<3>::next()

+{

+    gotoNextCommon(*this);

+}

+

+void TemplatedSparseIterator<3>::calculateCommonFlags()

+{

+    mCommonFlags = mSparseFlags & mHybridFlags_1 & mHybridFlags_2;

+}

+

+void TemplatedSparseIterator<3>::getFlags()

+{

+    mSparseFlags = mSparse.mIndexBitFlags[mBigIndex];

+    mHybridFlags_1 = mHybrid_1.mIndexBitFlags[mBigIndex];

+    mHybridFlags_2 = mHybrid_2.mIndexBitFlags[mBigIndex];

+}

+

+////////////////////////////////////////////////////////////////////////////////

+////////////////////////////////////////////////////////////////////////////////

+////////////////////////////////////////////////////////////////////////////////

+////////////////////////////////////////////////////////////////////////////////

+

 SparseIterator::SparseIterator(const SparseVector &v)

     :

 mSparse(v),

@@ -4,7 +4,111 @@

 #include "HybridVector.h"

 #include "SparseVector.h"

-// TODO make these nicer with templates - make sure no performance lost

+template <unsigned N, class Iter>

+float get(const Iter &it);

+

+// only allow this class to constructed with N=1,2,3

+template <unsigned N>

+class TemplatedSparseIterator

+{

+private:

+    TemplatedSparseIterator() {}

+};

+

+template<>

+class TemplatedSparseIterator<1>

+{

+public:

+

+    TemplatedSparseIterator(const SparseVector &v);

+

+    bool atEnd() const;

+    void next();

+

+private:

+

+    friend float get<1>(const TemplatedSparseIterator<1> &it);

+

+    const SparseVector &mSparse;

+    unsigned mSparseIndex;

+};

+

+template<>

+class TemplatedSparseIterator<2>

+{

+public:

+

+    TemplatedSparseIterator(const SparseVector &v, const HybridVector &h);

+

+    bool atEnd() const;

+    void next();

+    void calculateCommonFlags();

+    void getFlags();

+

+private:

+

+    template <class Iter>

+    friend void gotoNextCommon(Iter &it);

+

+    friend float get<1>(const TemplatedSparseIterator<2> &it);

+    friend float get<2>(const TemplatedSparseIterator<2> &it);

+

+    const SparseVector &mSparse;  

+    const HybridVector &mHybrid_1;

+

+    uint64_t mSparseFlags;

+    uint64_t mHybridFlags_1;

+    uint64_t mCommonFlags;

+

+    unsigned mTotalIndices;

+    unsigned mBigIndex;

+    unsigned mSmallIndex;

+    unsigned mSparseIndex;

+    bool mAtEnd;

+};

+

+template<>

+class TemplatedSparseIterator<3>

+{

+public:

+

+    TemplatedSparseIterator(const SparseVector &v, const HybridVector &h1,

+    const HybridVector &h2);

+

+    bool atEnd() const;

+    void next();

+    void calculateCommonFlags();

+    void getFlags();

+

+private:

+

+    template <class Iter>

+    friend void gotoNextCommon(Iter &it);

+

+    friend float get<1>(const TemplatedSparseIterator<3> &it);

+    friend float get<2>(const TemplatedSparseIterator<3> &it);

+    friend float get<3>(const TemplatedSparseIterator<3> &it);

+

+    const SparseVector &mSparse;  

+    const HybridVector &mHybrid_1;

+    const HybridVector &mHybrid_2;

+

+    uint64_t mSparseFlags;

+    uint64_t mHybridFlags_1;

+    uint64_t mHybridFlags_2;

+    uint64_t mCommonFlags;

+

+    unsigned mTotalIndices;

+    unsigned mBigIndex;

+    unsigned mSmallIndex;

+    unsigned mSparseIndex;

+    bool mAtEnd;

+};

+

+////////////////////////////////////////////////////////////////////////////////

+////////////////////////////////////////////////////////////////////////////////

+////////////////////////////////////////////////////////////////////////////////

+////////////////////////////////////////////////////////////////////////////////

 class SparseIterator

 {