@@ -43,14 +43,14 @@ public:

     template<typename T>

     friend Archive& operator<<(Archive &ar, T val)

     {

-        ar.mStream.write(reinterpret_cast<char*>(&val), sizeof(T));

+        ar.mStream.write(reinterpret_cast<char*>(&val), sizeof(T)); // NOLINT

         return ar;

     }

     template<typename T>

     friend Archive& operator>>(Archive &ar, T &val)

     {

-        ar.mStream.read(reinterpret_cast<char*>(&val), sizeof(T));

+        ar.mStream.read(reinterpret_cast<char*>(&val), sizeof(T)); // NOLINT

         return ar;

     }

 };

@@ -196,5 +196,5 @@ void AtomicDomain::flushCache()

 // O(logN)

 void AtomicDomain::updateMass(uint64_t pos, float newMass)

 {

-    mAtoms[mAtomPositions.at(pos)].mass = newMass;

+    mAtoms[mAtomPositions.at(pos)].mass = newMass; // TODO at is C++11

 }

@@ -4,10 +4,11 @@

 #include "Archive.h"

 #include <boost/unordered_set.hpp>

-#include <stdint.h>

+

 #include <cstddef>

-#include <vector>

 #include <map>

+#include <stdint.h>

+#include <vector>

 class AtomicDomain;

@@ -1,4 +1,9 @@

 #include "GapsRunner.h"

+#include "math/SIMD.h"

+

+#ifdef __GAPS_OPENMP__

+    #include <omp.h>

+#endif

 GapsRunner::GapsRunner(const Rcpp::NumericMatrix &D, const Rcpp::NumericMatrix &S,

 unsigned nFactor, unsigned nEquil, unsigned nCool, unsigned nSample,

@@ -55,7 +55,7 @@ public:

     uint32_t mSeed;

     bpt::ptime mLastCheckpoint;

-    long mCheckpointInterval;

+    int64_t mCheckpointInterval;

     std::string mCheckpointFile;

     unsigned mNumUpdatesA;

@@ -45,7 +45,7 @@ void AmplitudeGibbsSampler::updateAPMatrix(unsigned row, unsigned col, float del

     unsigned size = mAPMatrix.nCol();

     gaps::simd::packedFloat pOther, pAP;

-    gaps::simd::Index i = 0;

+    gaps::simd::Index i(0);

     gaps::simd::packedFloat pDelta(delta);

     for (; i <= size - i.increment(); ++i)

     {

@@ -142,7 +142,7 @@ void PatternGibbsSampler::updateAPMatrix(unsigned row, unsigned col, float delta

     unsigned size = mAPMatrix.nRow();

     gaps::simd::packedFloat pOther, pAP;

-    gaps::simd::Index i = 0;

+    gaps::simd::Index i(0);

     gaps::simd::packedFloat pDelta(delta);

     for (; i <= size - i.increment(); ++i)

     {

@@ -1,15 +1,16 @@

 #ifndef __COGAPS_GIBBS_SAMPLER_H__

 #define __COGAPS_GIBBS_SAMPLER_H__

-#include "GapsAssert.h"

 #include "Archive.h"

+#include "AtomicDomain.h"

+#include "GapsAssert.h"

+#include "ProposalQueue.h"

 #include "data_structures/Matrix.h"

-#include "math/Random.h"

 #include "math/Algorithms.h"

-#include "ProposalQueue.h"

-#include "AtomicDomain.h"

+#include "math/Random.h"

 #include <Rcpp.h>

+

 #include <algorithm>

 // forward declarations needed for friend classes

@@ -56,8 +57,8 @@ protected:

     void death(uint64_t pos, float mass, unsigned row, unsigned col);

     void move(uint64_t src, float mass, uint64_t dest, unsigned r1, unsigned c1,

         unsigned r2, unsigned c2);

-    void exchange(uint64_t p1, float mass1, uint64_t p2, float mass2,

-        unsigned r1, unsigned c1, unsigned r2, unsigned c2);

+    void exchange(uint64_t p1, float m1, uint64_t p2, float m2, unsigned r1,

+        unsigned c1, unsigned r2, unsigned c2);

     float gibbsMass(unsigned row, unsigned col, float mass);

     float gibbsMass(unsigned r1, unsigned c1, float m1, unsigned r2,

@@ -65,7 +66,7 @@ protected:

     void addMass(uint64_t pos, float mass, unsigned row, unsigned col);

     void removeMass(uint64_t pos, float mass, unsigned row, unsigned col);

-    float updateAtomMass(uint64_t pos, float mass, float delta);

+    bool updateAtomMass(uint64_t pos, float mass, float delta);

     void acceptExchange(uint64_t p1, float m1, float d1, uint64_t p2, float m2,

         float d2, unsigned r1, unsigned c1, unsigned r2, unsigned c2);

@@ -150,8 +151,8 @@ GibbsSampler<T, MatA, MatB>::GibbsSampler(const Rcpp::NumericMatrix &D,

 const Rcpp::NumericMatrix &S, unsigned nrow, unsigned ncol, float alpha)

     :

 mMatrix(nrow, ncol), mOtherMatrix(NULL), mDMatrix(D), mSMatrix(S),

-mAPMatrix(D.nrow(), D.ncol()), mQueue(nrow * ncol, alpha),

-mAnnealingTemp(0.f), mNumRows(nrow), mNumCols(ncol),

+mAPMatrix(D.nrow(), D.ncol()), mQueue(nrow * ncol, alpha), mLambda(0.f),

+mMaxGibbsMass(0.f), mAnnealingTemp(0.f), mNumRows(nrow), mNumCols(ncol),

 mAvgQueue(0.f), mNumQueues(0.f)

 {

     mBinSize = std::numeric_limits<uint64_t>::max()

@@ -343,16 +344,13 @@ float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

             acceptExchange(p1, m1, delta, p2, m2, -delta, r1, c1, r2, c2);

             return;

         }

-        else // normal case

+        float deltaLL = impl()->computeDeltaLL(r1, c1, delta, r2, c2, -delta);

+        float priorLL = (pOld == 0.f) ? 1.f : pOld / pNew;

+        float u = std::log(gaps::random::uniform() * priorLL);

+        if (u < deltaLL * mAnnealingTemp)

         {

-            float deltaLL = impl()->computeDeltaLL(r1, c1, delta, r2, c2, -delta);

-            float priorLL = (pOld == 0.f) ? 1.f : pOld / pNew;

-            float u = std::log(gaps::random::uniform() * priorLL);

-            if (u < deltaLL * mAnnealingTemp)

-            {

-                acceptExchange(p1, m1, delta, p2, m2, -delta, r1, c1, r2, c2);

-                return;

-            }

+            acceptExchange(p1, m1, delta, p2, m2, -delta, r1, c1, r2, c2);

+            return;

         }

     }

     mQueue.rejectDeath();

@@ -361,7 +359,7 @@ float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

 // helper function for acceptExchange, used to conditionally update the mass

 // at a single position, deleting the atom if the new mass is too small

 template <class T, class MatA, class MatB>

-float GibbsSampler<T, MatA, MatB>::updateAtomMass(uint64_t pos, float mass,

+bool GibbsSampler<T, MatA, MatB>::updateAtomMass(uint64_t pos, float mass,

 float delta)

 {

     if (mass + delta < gaps::algo::epsilon)

@@ -370,11 +368,8 @@ float delta)

         mQueue.acceptDeath();

         return false;

     }

-    else

-    {

-        mDomain.updateMass(pos, mass + delta);

-        return true;

-    }

+    mDomain.updateMass(pos, mass + delta);

+    return true;

 }

 // helper function for exchange step, updates the atomic domain, matrix, and

@@ -1,4 +1,7 @@

-PKG_CXXFLAGS =  -DNSIMD -DGAPS_NDEBUG -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0 -Werror

+PKG_CFLAGS = $(SHLIB_OPENMP_CFLAGS)

+PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(SHLIB_OPENMP_CFLAGS)

+PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -DNSIMD -DGAPS_NDEBUG -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0 -Werror

+

 OBJECTS =   AtomicDomain.o \

             Cogaps.o \

             GapsRunner.o \

@@ -6,8 +6,9 @@

 #include "data_structures/EfficientSets.h"

 #include <boost/unordered_set.hpp>

-#include <stdint.h>

+

 #include <cstddef>

+#include <stdint.h>

 struct AtomicProposal

 {

@@ -65,8 +66,8 @@ private:

 public:

     ProposalQueue(unsigned nBins, float alpha)

-        : mMinAtoms(0), mMaxAtoms(0), mNumBins(nBins), mAlpha(alpha),

-        mUseCachedRng(false), mU1(0.f), mU2(0.f)

+        : mMinAtoms(0), mMaxAtoms(0), mNumBins(nBins), mDimensionSize(0),

+        mDomainSize(0), mAlpha(alpha), mUseCachedRng(false), mU1(0.f), mU2(0.f)

     {}

     // set parameters

@@ -33,8 +33,6 @@ private:

 public:

-    IntDenseOrderedSet() {}

-

     void insert(uint64_t p) {mVec.push_back(p);}

     void clear() {mVec.clear();}

@@ -52,12 +52,6 @@ RowMatrix::RowMatrix(const Rcpp::NumericMatrix &rmat)

     }

 }

-RowMatrix& RowMatrix::operator=(const RowMatrix &mat)

-{

-    copyMatrix(*this, mat);

-    return *this;

-}

-

 RowMatrix& RowMatrix::operator=(const ColMatrix &mat)

 {

     copyMatrix(*this, mat);

@@ -131,12 +125,6 @@ ColMatrix ColMatrix::operator/(float val) const

     return mat;

 }

-ColMatrix& ColMatrix::operator=(const ColMatrix &mat)

-{

-    copyMatrix(*this, mat);

-    return *this;

-}

-

 ColMatrix& ColMatrix::operator=(const RowMatrix &mat)

 {

     copyMatrix(*this, mat);

@@ -22,7 +22,7 @@ private:

 public:

     RowMatrix(unsigned nrow, unsigned ncol);

-    RowMatrix(const Rcpp::NumericMatrix &rmat);

+    explicit RowMatrix(const Rcpp::NumericMatrix &rmat);

     template <class Parser>

     RowMatrix(Parser &p, unsigned nrow, unsigned ncol);

@@ -41,7 +41,6 @@ public:

     RowMatrix operator/(float val) const;

-    RowMatrix& operator=(const RowMatrix &mat);

     RowMatrix& operator=(const ColMatrix &mat);

     Rcpp::NumericMatrix rMatrix() const;

@@ -60,7 +59,7 @@ private:

 public:

     ColMatrix(unsigned nrow, unsigned ncol);

-    ColMatrix(const Rcpp::NumericMatrix &rmat);

+    explicit ColMatrix(const Rcpp::NumericMatrix &rmat);

     template <class Parser>

     ColMatrix(Parser &p, unsigned nrow, unsigned ncol);

@@ -79,7 +78,6 @@ public:

     ColMatrix operator/(float val) const;

-    ColMatrix& operator=(const ColMatrix &mat);

     ColMatrix& operator=(const RowMatrix &mat);

     Rcpp::NumericMatrix rMatrix() const;

@@ -20,7 +20,7 @@ private:

 public:

-    CsvParser(const std::string &path);

+    explicit CsvParser(const std::string &path);

     unsigned nRow() const { return mNumRows; }

     unsigned nCol() const { return mNumCols; }

@@ -17,7 +17,7 @@ private:

 public:

-    MtxParser(const std::string &path);

+    explicit MtxParser(const std::string &path);

     unsigned nRow() const { return mNumRows; }

     unsigned nCol() const { return mNumCols; }

@@ -20,7 +20,7 @@ private:

 public:

-    TsvParser(const std::string &path);

+    explicit TsvParser(const std::string &path);

     unsigned nRow() const { return mNumRows; }

     unsigned nCol() const { return mNumCols; }

@@ -111,8 +111,8 @@ AlphaParameters gaps::algo::alphaParameters(unsigned size, const float *D,

 const float *S, const float *AP, const float *mat)

 {

     gaps::simd::packedFloat ratio, pMat, pD, pAP, pS;

-    gaps::simd::packedFloat partialS = 0.f, partialSU = 0.f;

-    gaps::simd::Index i = 0;

+    gaps::simd::packedFloat partialS(0.f), partialSU(0.f);

+    gaps::simd::Index i(0);

     for (; i <= size - i.increment(); ++i)

     {   

         pMat.load(mat + i);

@@ -138,8 +138,8 @@ AlphaParameters gaps::algo::alphaParameters(unsigned size, const float *D,

 const float *S, const float *AP, const float *mat1, const float *mat2)

 {

     gaps::simd::packedFloat ratio, pMat1, pMat2, pD, pAP, pS;

-    gaps::simd::packedFloat partialS = 0.f, partialSU = 0.f;

-    gaps::simd::Index i = 0;

+    gaps::simd::packedFloat partialS(0.f), partialSU(0.f);

+    gaps::simd::Index i(0);

     for (; i <= size - i.increment(); ++i)

     {   

         pMat1.load(mat1 + i);

@@ -165,9 +165,9 @@ const float *S, const float *AP, const float *mat1, const float *mat2)

 float gaps::algo::deltaLL(unsigned size, const float *D, const float *S,

 const float *AP, const float *mat, float delta)

 {

-    const gaps::simd::packedFloat pDelta = delta, two = 2.f;

-    gaps::simd::packedFloat d, pMat, pD, pAP, pS, partialSum = 0.f;

-    gaps::simd::Index i = 0;

+    const gaps::simd::packedFloat pDelta(delta), two(2.f);

+    gaps::simd::packedFloat d, pMat, pD, pAP, pS, partialSum(0.f);

+    gaps::simd::Index i(0);

     for (; i <= size - i.increment(); ++i)

     {   

         pMat.load(mat + i);

@@ -190,9 +190,9 @@ float gaps::algo::deltaLL(unsigned size, const float *D, const float *S,

 const float *AP, const float *mat1, float delta1, const float *mat2,

 float delta2)

 {

-    const gaps::simd::packedFloat pDelta1 = delta1, pDelta2 = delta2, two = 2.f;

-    gaps::simd::packedFloat d, pMat1, pMat2, pD, pAP, pS, partialSum = 0.f;

-    gaps::simd::Index i = 0;

+    const gaps::simd::packedFloat pDelta1(delta1), pDelta2(delta2), two(2.f);

+    gaps::simd::packedFloat d, pMat1, pMat2, pD, pAP, pS, partialSum(0.f);

+    gaps::simd::Index i(0);

     for (; i <= size - i.increment(); ++i)

     {   

         pMat1.load(mat1 + i);

@@ -116,7 +116,7 @@ float gaps::algo::nonZeroMean(const GenericMatrix &mat)

             }

         }

     }

-    return sum / (float)nNonZeros;

+    return sum / static_cast<float>(nNonZeros);

 }

 template<class GenericMatrix>

@@ -129,8 +129,8 @@ const GenericMatrix &meanMat, unsigned nUpdates)

     {

         for (unsigned c = 0; c < retMat.nCol(); ++c)

         {

-            meanTerm = meanMat(r,c) * meanMat(r,c) / (float)nUpdates;

-            retMat(r,c) = std::sqrt((stdMat(r,c) - meanTerm) / ((float)nUpdates - 1.f));

+            meanTerm = meanMat(r,c) * meanMat(r,c) / static_cast<float>(nUpdates);

+            retMat(r,c) = std::sqrt((stdMat(r,c) - meanTerm) / (static_cast<float>(nUpdates) - 1.f));

         }

     }

     return retMat;

@@ -8,16 +8,17 @@

 // spawn a random generator for each thread

 // - no way to acheive consistent results across different nThreads

-#include <boost/random/uniform_01.hpp>

-#include <boost/random/uniform_real_distribution.hpp>

-#include <boost/random/normal_distribution.hpp>

-#include <boost/random/poisson_distribution.hpp>

-#include <boost/random/exponential_distribution.hpp>

-#include <boost/random/mersenne_twister.hpp>

-#include <boost/math/distributions/normal.hpp>

 #include <boost/math/distributions/exponential.hpp>

 #include <boost/math/distributions/gamma.hpp>

+#include <boost/math/distributions/normal.hpp>

+

+#include <boost/random/exponential_distribution.hpp>

+#include <boost/random/mersenne_twister.hpp>

+#include <boost/random/normal_distribution.hpp>

+#include <boost/random/poisson_distribution.hpp>

+#include <boost/random/uniform_01.hpp>

+#include <boost/random/uniform_real_distribution.hpp>

 #include <stdint.h>

@@ -35,59 +36,60 @@

 //typedef boost::random::mt19937 RNGType;

 typedef boost::random::mt11213b RNGType; // should be faster

-static std::vector<RNGType> rng;

-//static RNGType rng;

+static std::vector<RNGType>& rng()

+{

+    static std::vector<RNGType> allRngs(omp_get_max_threads());

+    return allRngs;

+}

 void gaps::random::save(Archive &ar)

 {

-    ar << rng;

+    //ar << rng;

 }

 void gaps::random::load(Archive &ar)

 {

-    ar >> rng;

+    //ar >> rng;

 }

 void gaps::random::setSeed(uint32_t seed)

 {

     unsigned n = omp_get_max_threads();

-    rng.push_back(RNGType());

-    rng[0].seed(seed);

+    rng().at(0).seed(seed);

     boost::random::uniform_int_distribution<uint32_t> seedDist(0,

         std::numeric_limits<uint32_t>::max());

     for (unsigned i = 1; i < n; ++i)

     {

-        rng.push_back(RNGType());

-        uint32_t newSeed = seedDist(rng[i-1]);

-        rng[i].seed(newSeed);

+        uint32_t newSeed = seedDist(rng().at(i-1));

+        rng().at(i).seed(newSeed);

     }

 }

 float gaps::random::normal(float mean, float var)

 {

     boost::random::normal_distribution<float> dist(mean, var);

-    return dist(rng[omp_get_thread_num()]);

+    return dist(rng().at(omp_get_thread_num()));

 }

 int gaps::random::poisson(float lambda)

 {

     boost::random::poisson_distribution<> dist(lambda);

-    return dist(rng[omp_get_thread_num()]);

+    return dist(rng().at(omp_get_thread_num()));

 }

 float gaps::random::exponential(float lambda)

 {

     boost::random::exponential_distribution<> dist(lambda);

-    return dist(rng[omp_get_thread_num()]);

+    return dist(rng().at(omp_get_thread_num()));

 }

 float gaps::random::uniform()

 {

     float ret = 0.f;

-    boost::random::uniform_01<RNGType&> u01_dist(rng[omp_get_thread_num()]);

+    boost::random::uniform_01<RNGType&> u01_dist(rng().at(omp_get_thread_num()));

     return u01_dist();

 }

@@ -97,18 +99,15 @@ float gaps::random::uniform(float a, float b)

     {

         return a;

     }

-    else

-    {

-        boost::random::uniform_real_distribution<> dist(a,b);

-        return dist(rng[omp_get_thread_num()]);

-    }

+    boost::random::uniform_real_distribution<> dist(a,b);

+    return dist(rng().at(omp_get_thread_num()));

 }

 uint64_t gaps::random::uniform64()

 {

     boost::random::uniform_int_distribution<uint64_t> dist(0,

         std::numeric_limits<uint64_t>::max());

-    return dist(rng[omp_get_thread_num()]);

+    return dist(rng().at(omp_get_thread_num()));

 }

 uint64_t gaps::random::uniform64(uint64_t a, uint64_t b)

@@ -117,11 +116,8 @@ uint64_t gaps::random::uniform64(uint64_t a, uint64_t b)

     {

         return a;

     }

-    else

-    {

-        boost::random::uniform_int_distribution<uint64_t> dist(a,b);

-        return dist(rng[omp_get_thread_num()]);

-    }

+    boost::random::uniform_int_distribution<uint64_t> dist(a,b);

+    return dist(rng().at(omp_get_thread_num()));

 }

 float gaps::random::d_gamma(float d, float shape, float scale)

@@ -142,11 +138,8 @@ float gaps::random::q_gamma(float q, float shape, float scale)

     {

         return Q_GAMMA_MIN_VALUE;

     }

-    else

-    {

-        boost::math::gamma_distribution<> gam(shape, scale);

-        return quantile(gam, q);

-    }

+    boost::math::gamma_distribution<> gam(shape, scale);

+    return quantile(gam, q);

 }

 float gaps::random::d_norm(float d, float mean, float sd)

@@ -3,40 +3,38 @@

 #include "../Archive.h"

+#include <fstream>

 #include <stdint.h>

 #include <vector>

-#include <fstream>

 namespace gaps

 {

-

-namespace random

-{

-    void setSeed(uint32_t seed);

-

-    int poisson(float lambda);

-    uint64_t uniform64();

-    uint64_t uniform64(uint64_t a, uint64_t b);

-

-    float uniform();

-    float uniform(float a, float b);

-    float normal(float mean, float var);

-    float exponential(float lambda);

-

-    float d_gamma(float d, float shape, float scale);

-    float p_gamma(float p, float shape, float scale);

-    float q_gamma(float q, float shape, float scale);

-    float d_norm(float d, float mean, float sd);

-    float q_norm(float q, float mean, float sd);

-    float p_norm(float p, float mean, float sd);

-

-    float inverseNormSample(float a, float b, float mean, float sd);

-    float inverseGammaSample(float a, float b, float mean, float sd);

-

-    void save(Archive &ar);

-    void load(Archive &ar);

-}

-

-}

+    namespace random

+    {

+        void setSeed(uint32_t seed);

+

+        float uniform();

+        float uniform(float a, float b);

+        uint64_t uniform64();

+        uint64_t uniform64(uint64_t a, uint64_t b);

+

+        float exponential(float lambda);

+        float normal(float mean, float var);

+        int poisson(float lambda);

+

+        float inverseNormSample(float a, float b, float mean, float sd);

+        float inverseGammaSample(float a, float b, float mean, float sd);

+

+        float d_gamma(float d, float shape, float scale);

+        float p_gamma(float p, float shape, float scale);

+        float q_gamma(float q, float shape, float scale);

+        float d_norm(float d, float mean, float sd);

+        float q_norm(float q, float mean, float sd);

+        float p_norm(float p, float mean, float sd);

+

+        void save(Archive &ar);

+        void load(Archive &ar);

+    } // namespace random

+} // namespace gaps

 #endif

\ No newline at end of file

@@ -68,7 +68,7 @@ namespace simd

     const unsigned index_increment = 1;

     #define SET_SCALAR(x) x

     #define LOAD_PACKED(x) *(x)

-    #define STORE_PACKED(p,x) *p = x

+    #define STORE_PACKED(p,x) *(p) = (x)

     #define ADD_PACKED(a,b) ((a)+(b))

     #define SUB_PACKED(a,b) ((a)-(b))

     #define MUL_PACKED(a,b) ((a)*(b))

@@ -78,10 +78,13 @@ namespace simd

 class Index

 {

 private:

+

     unsigned index;

+

 public:

-    Index(unsigned i) : index(i) {}

-    void operator=(unsigned val) { index = val; }

+

+    explicit Index(unsigned i) : index(i) {}

+    Index& operator=(unsigned val) { index = val; return *this; }

     bool operator<(unsigned comp) { return index < comp; }

     bool operator<=(unsigned comp) { return index <= comp; }

     void operator++() { index += index_increment; }

@@ -102,16 +105,16 @@ private:

 public:

-    packedFloat() {}

-    packedFloat(float val) : mData(SET_SCALAR(val)) {}

+    packedFloat() : mData() {}

+    explicit packedFloat(float val) : mData(SET_SCALAR(val)) {}

 #if defined( __GAPS_SSE__ ) || defined( __GAPS_AVX__ )

-    packedFloat(gaps_packed_t val) : mData(val) {}

+    explicit packedFloat(gaps_packed_t val) : mData(val) {}

 #endif

-    packedFloat operator+(packedFloat b) const { return ADD_PACKED(mData, b.mData); }

-    packedFloat operator-(packedFloat b) const { return SUB_PACKED(mData, b.mData); }

-    packedFloat operator*(packedFloat b) const { return MUL_PACKED(mData, b.mData); }

-    packedFloat operator/(packedFloat b) const { return DIV_PACKED(mData, b.mData); }

+    packedFloat operator+(packedFloat b) const { return packedFloat(ADD_PACKED(mData, b.mData)); }

+    packedFloat operator-(packedFloat b) const { return packedFloat(SUB_PACKED(mData, b.mData)); }

+    packedFloat operator*(packedFloat b) const { return packedFloat(MUL_PACKED(mData, b.mData)); }

+    packedFloat operator/(packedFloat b) const { return packedFloat(DIV_PACKED(mData, b.mData)); }

     void operator+=(packedFloat val) { mData = ADD_PACKED(mData, val.mData); }

     void load(const float *ptr) { mData = LOAD_PACKED(ptr); }