ode-toolbox

ODE integration tools
git clone https://git.0xfab.ch/ode-toolbox.git
Log | Files | Refs | README | LICENSE

commit 1c50b530ffe3e5dc4e75aa95e0b3e466408f3f17
parent 51af79e158759bbea9ab09fa5853093a91937f5f
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date:   Sun, 11 Jun 2017 12:03:31 -0700

got rid of StateVector type by using templated operators in LightVector

Diffstat:
MTimeStepper/orderVerification/orderVerification.cpp | 2+-
Mcommon.h | 314+++++++++++++++++++++++++++++++++++++++++++++++++------------------------------
2 files changed, 195 insertions(+), 121 deletions(-)

diff --git a/TimeStepper/orderVerification/orderVerification.cpp b/TimeStepper/orderVerification/orderVerification.cpp @@ -64,7 +64,7 @@ public: // types we use typedef State<Real,1> State1; -using vec_t = StateVector<State1>; +using vec_t = LightVector<State1>; using baseKernel_t = KernelBase<vec_t>*; using baseStepper_t = StepperBase<vec_t>*; diff --git a/common.h b/common.h @@ -31,6 +31,13 @@ template <> inline double mysqrt(double a) { return sqrt(a); } +// state vector type +template <typename T, int _SIZE> +class State; + +template <typename T, int _SIZE> +inline std::ostream& operator<<(std::ostream& out, const State<T,_SIZE>& s); + template <typename T, int _SIZE> struct State { @@ -47,18 +54,13 @@ struct State inline T& operator[](const size_t i) { assert(i<_SIZE); return s[i]; } inline T operator[](const size_t i) const { assert(i<_SIZE); return s[i]; } + // vector inline State& operator=(const State& rhs) { for (int i = 0; i < _SIZE; ++i) s[i] = rhs.s[i]; return *this; } - inline State& operator=(const Real c) - { - for (int i = 0; i < _SIZE; ++i) - s[i] = c; - return *this; - } inline State& operator+=(const State& rhs) { for (int i = 0; i < _SIZE; ++i) @@ -106,28 +108,105 @@ struct State State thisCopy(*this); return thisCopy /= rhs; } + + // scalar + inline State& operator=(const Real c) + { + for (int i = 0; i < _SIZE; ++i) + s[i] = c; + return *this; + } + inline State& operator+=(const Real c) + { + for (int i = 0; i < _SIZE; ++i) + s[i] += c; + return *this; + } + inline State& operator-=(const Real c) + { + for (int i = 0; i < _SIZE; ++i) + s[i] -= c; + return *this; + } inline State& operator*=(const Real c) { for (int i = 0; i < _SIZE; ++i) s[i] *= c; return *this; } + inline State& operator/=(const Real c) + { + assert(c != static_cast<Real>(0.0)); + for (int i = 0; i < _SIZE; ++i) + s[i] /= c; + return *this; + } + inline State operator+(const Real c) const + { + State thisCopy(*this); + return thisCopy += c; + } + inline State operator-(const Real c) const + { + State thisCopy(*this); + return thisCopy -= c; + } inline State operator*(const Real c) const { State thisCopy(*this); return thisCopy *= c; } + inline State operator/(const Real c) const + { + State thisCopy(*this); + return thisCopy /= c; + } + + // output operator + friend std::ostream& operator<< <>(std::ostream& out, const State& s); }; template <typename T, int _SIZE> -inline State<T,_SIZE> operator*(Real const c, State<T,_SIZE> const& SV) +inline State<T,_SIZE> operator+(const Real c, State<T,_SIZE> const& SV) +{ + State<T,_SIZE> SVcopy(SV); + return SVcopy += c; +} + +template <typename T, int _SIZE> +inline State<T,_SIZE> operator-(const Real c, State<T,_SIZE> const& SV) +{ + State<T,_SIZE> SVcopy(SV); + return SVcopy -= c; +} + +template <typename T, int _SIZE> +inline State<T,_SIZE> operator*(const Real c, State<T,_SIZE> const& SV) { State<T,_SIZE> SVcopy(SV); return SVcopy *= c; } +template <typename T, int _SIZE> +inline std::ostream& operator<<(std::ostream& out, const State<T,_SIZE>& s) +{ + assert(_SIZE > 0); + out << "(" << s[0]; + for (int i = 1; i < _SIZE; ++i) + out << ", " << s[i]; + out << ")"; + return out; +} + +// generic 'lightweight' vector type with stencil support template <typename T, int _SS=0, int _SE=0> +class LightVector; + +template <typename T, int _SS, int _SE> +inline std::ostream& operator<<(std::ostream& out, const LightVector<T,_SS,_SE>& v); + +template <typename T, int _SS, int _SE> class LightVector { bool _bAllocated; @@ -160,16 +239,16 @@ class LightVector _data[i] = src[i]; } - inline void _clear() + inline void _clear(const Real v0=0.0) { #pragma omp parallel for LOOPSCHED for (size_t i = 0; i < _N+_SE-_SS; ++i) - _data[i] = 0.0; + _data[i] = v0; } public: LightVector() : _bAllocated(false), _N(0), _data(nullptr) {} - LightVector(const int n) : _bAllocated(false), _N(n), _data(nullptr) { _alloc(n+(_SE-_SS)); _clear(); } + LightVector(const int n, const Real v0=0.0) : _bAllocated(false), _N(n), _data(nullptr) { _alloc(n+(_SE-_SS)); _clear(v0); } LightVector(const LightVector& rhs) : _bAllocated(false), _N(rhs._N), _data(nullptr) { _alloc(rhs._N+(_SE-_SS)); _copy(rhs._data); } virtual ~LightVector() { _dealloc(); } @@ -177,6 +256,7 @@ public: static const int SE = _SE; using DataType = T; + // elementwise vector operations LightVector& operator=(const LightVector& rhs) { if (this != &rhs) @@ -187,185 +267,179 @@ public: return *this; } - LightVector& operator=(const T& c) + LightVector& operator+=(const LightVector& c) { #pragma omp parallel for LOOPSCHED for (size_t i = 0; i < _N+_SE-_SS; ++i) - _data[i] = c; + _data[i] += c._data[i]; return *this; } - LightVector& operator*=(const T& c) + LightVector& operator-=(const LightVector& c) { #pragma omp parallel for LOOPSCHED for (size_t i = 0; i < _N+_SE-_SS; ++i) - _data[i] *= c; + _data[i] -= c._data[i]; return *this; } - LightVector& operator+=(const LightVector& c) + LightVector& operator*=(const LightVector& c) { #pragma omp parallel for LOOPSCHED for (size_t i = 0; i < _N+_SE-_SS; ++i) - _data[i] += c._data[i]; + _data[i] *= c._data[i]; return *this; } - inline void allocate(const size_t n) { _N=n; _alloc(n+(_SE-_SS)); _clear(); } - inline size_t size() const { return _N; } - inline size_t size_all() const { return _N+_SE-_SS; } - inline void clear() { _clear(); } - inline T& operator[](const size_t i) { assert(i<_N+_SE-_SS); return _data[i]; } - inline const T& operator[](const size_t i) const { assert(i<_N+_SE-_SS); return _data[i]; } - inline T& operator()(const int ix) { assert((ix-_SS >= 0)&&(ix-_SS < static_cast<int>(_N)+_SE-_SS)); return _data[ix-_SS]; } - inline const T& operator()(const int ix) const { assert((ix-_SS >= 0)&&(ix-_SS < static_cast<int>(_N)+_SE-_SS)); return _data[ix-_SS]; } - inline T* data() { return _data; } - inline const T * const data() const { return _data; } -}; - -template <typename T> -class StateVector; - -template <typename T> -inline std::ostream& operator<<(std::ostream& out, StateVector<T> const& SV); - -template <typename T> -class StateVector -{ -public: - StateVector(size_t const N=1) : m_state(N), m_N(N) {} - StateVector(StateVector const& rhs) : m_state(rhs.m_N), m_N(rhs.m_N) { _copy(rhs); } - ~StateVector() {} - - using DataType = T; - - // accessors - inline T& operator[](const int i) + LightVector& operator/=(const LightVector& c) { - assert(0 <= i); - assert(i < m_N); - return m_state[i]; +#pragma omp parallel for LOOPSCHED + for (size_t i = 0; i < _N+_SE-_SS; ++i) + _data[i] /= c._data[i]; + return *this; } - inline T operator[](const int i) const + LightVector operator+(const LightVector& c) const { - assert(0 <= i); - assert(i < m_N); - return m_state[i]; + LightVector thisCopy(*this); + return thisCopy += c; } - inline Real* const data() { return &m_state[0].s[0]; } - inline Real const * const data() const { return &m_state[0].s[0]; } - - inline size_t size() const { return m_N; } - - // operators - StateVector& operator=(StateVector const& rhs) + LightVector operator-(const LightVector& c) const { - if (this != &rhs) _copy(rhs); - return *this; + LightVector thisCopy(*this); + return thisCopy -= c; } - StateVector& operator=(Real const rhs) + + LightVector operator*(const LightVector& c) const { - for (size_t i = 0; i < m_state.size(); ++i) - m_state[i] = rhs; - return *this; + LightVector thisCopy(*this); + return thisCopy *= c; } - StateVector& operator+=(StateVector const& rhs) + LightVector operator/(const LightVector& c) const { - for (size_t j=0; j < m_state.size(); ++j) - m_state[j] += rhs.m_state[j]; - return *this; + LightVector thisCopy(*this); + return thisCopy /= c; } - StateVector& operator-=(StateVector const& rhs) + // other rhs types (NOTE: parameter are assed by value. The U-types are + // meant to be built-in or cheap to copy) + template <typename U> + LightVector& operator=(const U c) { - for (size_t j=0; j < m_state.size(); ++j) - m_state[j] -= rhs.m_state[j]; +#pragma omp parallel for LOOPSCHED + for (size_t i = 0; i < _N+_SE-_SS; ++i) + _data[i] = c; return *this; } - StateVector& operator*=(StateVector const& rhs) + template <typename U> + LightVector& operator+=(const U c) { - for (size_t j=0; j < m_state.size(); ++j) - m_state[j] *= rhs.m_state[j]; +#pragma omp parallel for LOOPSCHED + for (size_t i = 0; i < _N+_SE-_SS; ++i) + _data[i] += c; return *this; } - StateVector& operator/=(StateVector const& rhs) + template <typename U> + LightVector& operator-=(const U c) { - for (size_t j=0; j < m_state.size(); ++j) - m_state[j] /= rhs.m_state[j]; +#pragma omp parallel for LOOPSCHED + for (size_t i = 0; i < _N+_SE-_SS; ++i) + _data[i] -= c; return *this; } - inline StateVector operator+(StateVector const& rhs) const + template <typename U> + LightVector& operator*=(const U c) { - StateVector thisCopy(*this); - return thisCopy += rhs; +#pragma omp parallel for LOOPSCHED + for (size_t i = 0; i < _N+_SE-_SS; ++i) + _data[i] *= c; + return *this; } - inline StateVector operator-(StateVector const& rhs) const + template <typename U> + LightVector& operator/=(const U c) { - StateVector thisCopy(*this); - return thisCopy -= rhs; +#pragma omp parallel for LOOPSCHED + for (size_t i = 0; i < _N+_SE-_SS; ++i) + _data[i] /= c; + return *this; } - inline StateVector operator*(StateVector const& rhs) const + template <typename U> + LightVector operator+(const U c) const { - StateVector thisCopy(*this); - return thisCopy *= rhs; + LightVector thisCopy(*this); + return thisCopy += c; } - inline StateVector operator/(StateVector const& rhs) const + template <typename U> + LightVector operator-(const U c) const { - StateVector thisCopy(*this); - return thisCopy /= rhs; + LightVector thisCopy(*this); + return thisCopy -= c; } - inline StateVector& operator*=(Real const c) + template <typename U> + LightVector operator*(const U c) const { - for (size_t i = 0; i < m_state.size(); ++i) - m_state[i] *= c; - return *this; + LightVector thisCopy(*this); + return thisCopy *= c; } - inline StateVector operator*(Real const c) const + template <typename U> + LightVector operator/(const U c) const { - StateVector thisCopy(*this); - return thisCopy *= c; + LightVector thisCopy(*this); + return thisCopy /= c; } - friend std::ostream& operator<< <>(std::ostream& out, StateVector const& SV); - -private: - LightVector<T> m_state; - size_t m_N; + // output operator + friend std::ostream& operator<< <>(std::ostream& out, const LightVector& v); - inline void _copy(StateVector const& rhs) - { - assert(m_N == rhs.m_N); - m_N = rhs.m_N; - for(size_t j=0; j < m_state.size(); ++j) - m_state[j] = rhs.m_state[j]; - } + inline void allocate(const size_t n) { _N=n; _alloc(n+(_SE-_SS)); _clear(); } + inline size_t size() const { return _N; } + inline size_t size_all() const { return _N+_SE-_SS; } + inline void clear() { _clear(); } + inline T& operator[](const size_t i) { assert(i<_N+_SE-_SS); return _data[i]; } + inline const T& operator[](const size_t i) const { assert(i<_N+_SE-_SS); return _data[i]; } + inline T& operator()(const int ix) { assert((ix-_SS >= 0)&&(ix-_SS < static_cast<int>(_N)+_SE-_SS)); return _data[ix-_SS]; } + inline const T& operator()(const int ix) const { assert((ix-_SS >= 0)&&(ix-_SS < static_cast<int>(_N)+_SE-_SS)); return _data[ix-_SS]; } + inline T* data() { return _data; } + inline const T * const data() const { return _data; } }; -template <typename T> -inline std::ostream& operator<<(std::ostream& out, StateVector<T> const& SV) +template <typename T, typename U, int _SS, int _SE> +LightVector<T,_SS,_SE> operator+(const U c, const LightVector<T,_SS,_SE>& v) { - for (size_t j = 0; j < SV.m_state.size(); ++j) - for (size_t i = 0; i < T::SIZE; ++i) - out << '\t' << SV.m_state[j].s[i]; - return out; + LightVector<T,_SS,_SE> vCopy(v); + return vCopy += c; } -template <typename T> -inline StateVector<T> operator*(Real const c, StateVector<T> const& SV) +template <typename T, typename U, int _SS, int _SE> +LightVector<T,_SS,_SE> operator-(const U c, const LightVector<T,_SS,_SE>& v) { - StateVector<T> SVcopy(SV); - return SVcopy *= c; + LightVector<T,_SS,_SE> vCopy(v); + return vCopy -= c; +} + +template <typename T, typename U, int _SS, int _SE> +LightVector<T,_SS,_SE> operator*(const U c, const LightVector<T,_SS,_SE>& v) +{ + LightVector<T,_SS,_SE> vCopy(v); + return vCopy *= c; +} + +template <typename T, int _SS, int _SE> +inline std::ostream& operator<<(std::ostream& out, const LightVector<T,_SS,_SE>& v) +{ + for (size_t j = 0; j < v.size(); ++j) + out << v[j] << std::endl; + return out; } #endif /* COMMON_H_4JKNXGVM */