Expression Template Study

main.cpp

#include <sys/time.h>

#include <iostream>
#include <vector>

#include <boost/format.hpp>

using std::size_t;

#define DYNAMIC_ALLOCATE
//#define DEBUG
#define FORCE_LAZY

template<typename T, size_t length, bool enable_et = false>
class vec
{

private:
  typedef vec<T, length, false> _vec;

#ifdef DYNAMIC_ALLOCATE
  std::vector<T> data;
#else
  T data[length];
#endif

public:
  vec()
#ifdef DYNAMIC_ALLOCATE
    : data(length)
#endif
  {
#ifdef DEBUG
    std::cerr << boost::format("ctor:\t%1%\n") % this;
#endif
  };

  ~vec()
  {
#ifdef DEBUG
    std::cerr << boost::format("dtor:\t%1%\n") % this;
#endif
  }


  T& operator[](const size_t index) { return data[index]; };
  T operator[](const size_t index) const { return data[index]; };

  _vec operator+(const _vec& r) const
  {
    _vec ret;
    for(size_t i=0; i<length; ++i) {
      ret.data[i] = data[i] + r.data[i];
    }
    return ret;
  }

  _vec operator-(const _vec& r) const
  {
    _vec ret;
    for(size_t i=0; i<length; ++i) {
      ret.data[i] = data[i] - r.data[i];
    }
    return ret;
  }

  _vec operator*(const T r) const
  {
    _vec ret;
    for(size_t i=0; i<length; ++i) {
      ret.data[i] = data[i] * r;
    }
    return ret;
  }

  _vec operator/(const T r) const
  {
    _vec ret;
    for(size_t i=0; i<length; ++i) {
      ret.data[i] = data[i] / r;
    }
    return ret;
  }

  _vec& operator+=(const _vec r)
  {
    for(size_t i=0; i<length; ++i) {
      data[i] += r.data[i];
    }
    return *this;
  }

  _vec& operator-=(const _vec r)
  {
    for(size_t i=0; i<length; ++i) {
      data[i] -= r.data[i];
    }
    return *this;
  }

  _vec& operator*=(const T r)
  {
    for(size_t i=0; i<length; ++i) {
      data[i] *= r;
    }
    return *this;
  }

  _vec& operator/=(const T r)
  {
    for(size_t i=0; i<length; ++i) {
      data[i] /= r;
    }
    return *this;
  }

  _vec& operator=(const T r)
  {
    for(size_t i=0; i<length; ++i) {
      data[i] = r;
    }
    return *this;
  }

  _vec& operator=(const _vec r)
  {
    for(size_t i=0; i<length; ++i) {
      data[i] = r.data[i];
    }
    return *this;
  }



};


template<typename T, size_t length>
class vec<T, length, true>
{

private:
  typedef vec<T, length, true> _vec;

#ifdef DYNAMIC_ALLOCATE
  std::vector<T> data;
#else
  T data[length];
#endif



public:
  vec()
#ifdef DYNAMIC_ALLOCATE
    : data(length)
#endif
  {
#ifdef DEBUG
    std::cerr << boost::format("ctor:\t%1%\n") % this;
#endif
  };

  ~vec()
  {
#ifdef DEBUG
    std::cerr << boost::format("dtor:\t%1%\n") % this;
#endif
  }

  T& operator[](const size_t index) { return data[index]; };
  T operator[](const size_t index) const { return data[index]; };


  class plus
  {
  public:
    static T op(const T l, const T r) { return l + r; }
  };

  class minus
  {
  public:
    static T op(const T l, const T r) { return l - r; }
  };

  class cross
  {
  public:
    static T op(const T l, const T r) { return l * r; }
  };

  class slash
  {
  public:
    static T op(const T l, const T r) { return l / r; }
  };


  class mono
  {

  private:
    const T val;

  public:
    mono(const T val) : val(val) {};
    T operator[](const size_t index) const { return val; }

  };

  template<typename op, typename left, typename right>
  class binary
  {

  private:
    typedef binary<op, left, right> _binary;

    const left& l;
    const right& r;

  public:
    binary(const left& l, const right& r) : l(l), r(r) {};
    
    T operator[](const size_t index) const
    {
      return op::op(l[index], r[index]);
    };
    
    template<typename right_>
    binary<plus, _binary, right_> operator+(const right_& r_) const
    {
      return binary<plus, _binary, right_>(*this, r_);
    }

    template<typename right_>
    binary<minus, _binary, right_> operator-(const right_& r_) const
    {
      return binary<minus, _binary, right_>(*this, r_);
    }

    binary<cross, _binary, mono> operator*(const T r_) const
    {
      return binary<cross, _binary, mono>(*this, mono(r_));
    }

    binary<slash, _binary, mono> operator/(const T r_) const
    {
      return binary<slash, _binary, mono>(*this, mono(r_));
    }

  };



  template<typename right>
  binary<plus, _vec, right> operator+(const right& r) const
  {
    return binary<plus, _vec, right>(*this, r);
  }

  template<typename right>
  binary<minus, _vec, right> operator-(const right& r) const
  {
    return binary<minus, _vec, right>(*this, r);
  }

  binary<cross, _vec, mono> operator*(const T r) const
  {
    return binary<cross, _vec, mono>(*this, mono(r));
  }

  binary<cross, _vec, mono> operator/(const T r) const
  {
    return binary<slash, _vec, mono>(*this, mono(r));
  }


  template<typename op, typename left_, typename right_>
  _vec& operator+=(const binary<op, left_, right_>& r)
  {
    for(size_t i=0; i<length; ++i) {
      this->data[i] += r[i];
    }
    return *this;
  }

  template<typename op, typename left_, typename right_>
  _vec& operator-=(const binary<op, left_, right_>& r)
  {
    for(size_t i=0; i<length; ++i) {
      this->data[i] -= r[i];
    }
    return *this;
  }

  _vec& operator*=(const T r)
  {
    for(size_t i=0; i<length; ++i) {
      this->data[i] *= r;
    }
    return *this;
  }

  _vec& operator/=(const T r)
  {
    for(size_t i=0; i<length; ++i) {
      this->data[i] /= r;
    }
    return *this;
  }

  template<typename op, typename left_, typename right_>
  _vec& operator=(const binary<op, left_, right_>& r)
  {
    for(size_t i=0; i<length; ++i) {
      this->data[i] = r[i];
    }
    return *this;
  }


};

int main()
{
  using std::cerr;
  using std::endl;

  typedef double float_t;

  const float_t dt = 0.1;
  const double mups_pred = 80;
  const size_t loops = (size_t)(1e6 * mups_pred);

  const size_t dimension = 3;
  vec<float_t, dimension, false> q, 

  for(size_t i=0; i<dimension; ++i) {
    q[i] = 0;
    p[i] = 0;
  }
  q[0] = 0.5403023058681398;
  q[1] = 0.8414709848078965;

  p[1] = 0.9092974268256817;
  p[2] = -0.4161468365471424;

  timeval start_time, end_time;
  for(;;){
    gettimeofday(&start_time, NULL);

    for(size_t l=0; l<loops; ++l) {
#ifdef FORCE_LAZY
      for(size_t d=0; d<dimension; ++d)
	p[d] = p[d] - q[d] * dt;
      for(size_t d=0; d<dimension; ++d)
	q[d] = q[d] + p[d] * dt;
#else
      p = p - q * dt;
      q = q + p * dt;
#endif
    }

    gettimeofday(&end_time, NULL);

    const double elapsed_time = ( ( end_time.tv_sec - start_time.tv_sec ) + (end_time.tv_usec - start_time.tv_usec) * (1.0/1000000) );
    const double mups = 1e-6 * loops / elapsed_time;
    cerr << boost::format("mups:\t%1%\n") % mups;
    cerr << boost::format("q:\t%1%\t%2%\t%3%\n") % q[0] % q[1] % q[2];
    cerr << boost::format("p:\t%1%\t%2%\t%3%\n") % p[0] % p[1] % p[2];
    float_t energy=0;
    for(size_t i=0; i<dimension; ++i) {
      const float_t p_i = p[i];
      const float_t q_i = q[i];
      energy += p_i * p_i + q_i * q_i;
    }
    energy /= 2;
    cerr << boost::format("energy:\t%1%\n") % energy;
  }
  return 0;
};