Skip to content

Instantly share code, notes, and snippets.

@aont
Created Aug 15, 2011
Embed
What would you like to do?
Expression Template Study
#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;
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment