common.h (11376B)
1 /* File: common.h */ 2 /* Date: Mon Dec 21 18:00:32 2015 */ 3 /* Author: Fabian Wermelinger */ 4 /* Tag: common stuff */ 5 /* Copyright 2015 ETH Zurich. All Rights Reserved. */ 6 #ifndef COMMON_H_4JKNXGVM 7 #define COMMON_H_4JKNXGVM 8 9 #include <cassert> 10 #include <ostream> 11 #include <algorithm> 12 #include <cstdlib> 13 #include <cstring> 14 #include <cmath> 15 #include <omp.h> 16 17 #define LOOPSCHED 18 19 #ifndef _ALIGN_ 20 #define _ALIGN_ 16 21 #endif /* _ALIGN_ */ 22 23 #ifdef _FLOAT_PRECISION_ 24 using Real = float; 25 #else 26 using Real = double; 27 #endif 28 29 template <typename T> 30 inline T mysqrt(T a) { abort(); return a; } 31 32 template <> 33 inline float mysqrt(float a) { return sqrtf(a); } 34 template <> 35 inline double mysqrt(double a) { return sqrt(a); } 36 37 38 // state vector type 39 template <typename T, int _SIZE> 40 class State; 41 42 template <typename T, int _SIZE> 43 inline std::ostream& operator<<(std::ostream& out, const State<T,_SIZE>& s); 44 45 template <typename T, int _SIZE> 46 struct State 47 { 48 static constexpr int SIZE = _SIZE; 49 T s[_SIZE]; 50 51 State() : s{0} {} 52 State(const State& rhs) 53 { 54 for (int i = 0; i < _SIZE; ++i) 55 s[i] = rhs.s[i]; 56 } 57 58 inline T& operator[](const size_t i) { assert(i<_SIZE); return s[i]; } 59 inline T operator[](const size_t i) const { assert(i<_SIZE); return s[i]; } 60 61 // vector 62 inline State& operator=(const State& rhs) 63 { 64 for (int i = 0; i < _SIZE; ++i) 65 s[i] = rhs.s[i]; 66 return *this; 67 } 68 inline State& operator+=(const State& rhs) 69 { 70 for (int i = 0; i < _SIZE; ++i) 71 s[i] += rhs.s[i]; 72 return *this; 73 } 74 inline State& operator-=(const State& rhs) 75 { 76 for (int i = 0; i < _SIZE; ++i) 77 s[i] -= rhs.s[i]; 78 return *this; 79 } 80 inline State& operator*=(const State& rhs) 81 { 82 for (int i = 0; i < _SIZE; ++i) 83 s[i] *= rhs.s[i]; 84 return *this; 85 } 86 inline State& operator/=(const State& rhs) 87 { 88 for (int i = 0; i < _SIZE; ++i) 89 { 90 assert(rhs.s[i] != 0); 91 s[i] /= rhs.s[i]; 92 } 93 return *this; 94 } 95 inline State operator+(const State& rhs) const 96 { 97 State thisCopy(*this); 98 return thisCopy += rhs; 99 } 100 inline State operator-(const State& rhs) const 101 { 102 State thisCopy(*this); 103 return thisCopy -= rhs; 104 } 105 inline State operator*(const State& rhs) const 106 { 107 State thisCopy(*this); 108 return thisCopy *= rhs; 109 } 110 inline State operator/(const State& rhs) const 111 { 112 State thisCopy(*this); 113 return thisCopy /= rhs; 114 } 115 116 // scalar 117 inline State& operator=(const Real c) 118 { 119 for (int i = 0; i < _SIZE; ++i) 120 s[i] = c; 121 return *this; 122 } 123 inline State& operator+=(const Real c) 124 { 125 for (int i = 0; i < _SIZE; ++i) 126 s[i] += c; 127 return *this; 128 } 129 inline State& operator-=(const Real c) 130 { 131 for (int i = 0; i < _SIZE; ++i) 132 s[i] -= c; 133 return *this; 134 } 135 inline State& operator*=(const Real c) 136 { 137 for (int i = 0; i < _SIZE; ++i) 138 s[i] *= c; 139 return *this; 140 } 141 inline State& operator/=(const Real c) 142 { 143 assert(c != static_cast<Real>(0.0)); 144 for (int i = 0; i < _SIZE; ++i) 145 s[i] /= c; 146 return *this; 147 } 148 inline State operator+(const Real c) const 149 { 150 State thisCopy(*this); 151 return thisCopy += c; 152 } 153 inline State operator-(const Real c) const 154 { 155 State thisCopy(*this); 156 return thisCopy -= c; 157 } 158 inline State operator*(const Real c) const 159 { 160 State thisCopy(*this); 161 return thisCopy *= c; 162 } 163 inline State operator/(const Real c) const 164 { 165 State thisCopy(*this); 166 return thisCopy /= c; 167 } 168 169 // output operator 170 friend std::ostream& operator<< <>(std::ostream& out, const State& s); 171 }; 172 173 template <typename T, int _SIZE> 174 inline State<T,_SIZE> operator+(const Real c, State<T,_SIZE> const& SV) 175 { 176 State<T,_SIZE> SVcopy(SV); 177 return SVcopy += c; 178 } 179 180 template <typename T, int _SIZE> 181 inline State<T,_SIZE> operator-(const Real c, State<T,_SIZE> const& SV) 182 { 183 State<T,_SIZE> SVcopy(SV); 184 return SVcopy -= c; 185 } 186 187 template <typename T, int _SIZE> 188 inline State<T,_SIZE> operator*(const Real c, State<T,_SIZE> const& SV) 189 { 190 State<T,_SIZE> SVcopy(SV); 191 return SVcopy *= c; 192 } 193 194 template <typename T, int _SIZE> 195 inline std::ostream& operator<<(std::ostream& out, const State<T,_SIZE>& s) 196 { 197 assert(_SIZE > 0); 198 out << "(" << s[0]; 199 for (int i = 1; i < _SIZE; ++i) 200 out << ", " << s[i]; 201 out << ")"; 202 return out; 203 } 204 205 206 // generic 'lightweight' vector type with stencil support 207 template <typename T, int _SS=0, int _SE=0> 208 class LightVector; 209 210 template <typename T, int _SS, int _SE> 211 inline std::ostream& operator<<(std::ostream& out, const LightVector<T,_SS,_SE>& v); 212 213 template <typename T, int _SS, int _SE> 214 class LightVector 215 { 216 bool _bAllocated; 217 size_t _N; // number of elements w/o _SS and _SE 218 T* _data; 219 220 inline void _alloc(const size_t n) 221 { 222 if (!_bAllocated) 223 { 224 void* pmem; 225 posix_memalign(&pmem, _ALIGN_, n*sizeof(T)); 226 _bAllocated = true; 227 _data = static_cast<T*>(pmem); 228 } 229 } 230 231 inline void _dealloc() 232 { 233 if (_bAllocated) free(_data); 234 _N = 0; 235 _bAllocated = false; 236 _data = nullptr; 237 } 238 239 inline void _copy(const T* const src) 240 { 241 #pragma omp parallel for LOOPSCHED 242 for (size_t i = 0; i < _N+_SE-_SS; ++i) 243 _data[i] = src[i]; 244 } 245 246 inline void _clear(const Real v0=0.0) 247 { 248 #pragma omp parallel for LOOPSCHED 249 for (size_t i = 0; i < _N+_SE-_SS; ++i) 250 _data[i] = v0; 251 } 252 253 public: 254 LightVector() : _bAllocated(false), _N(0), _data(nullptr) {} 255 LightVector(const int n, const Real v0=0.0) : _bAllocated(false), _N(n), _data(nullptr) { _alloc(n+(_SE-_SS)); _clear(v0); } 256 LightVector(const LightVector& rhs) : _bAllocated(false), _N(rhs._N), _data(nullptr) { _alloc(rhs._N+(_SE-_SS)); _copy(rhs._data); } 257 virtual ~LightVector() { _dealloc(); } 258 259 static const int SS = _SS; 260 static const int SE = _SE; 261 using DataType = T; 262 263 // elementwise vector operations 264 LightVector& operator=(const LightVector& rhs) 265 { 266 if (this != &rhs) 267 { 268 assert(_N+_SE-_SS == rhs._N+_SE-_SS); 269 _copy(rhs._data); 270 } 271 return *this; 272 } 273 274 LightVector& operator+=(const LightVector& c) 275 { 276 #pragma omp parallel for LOOPSCHED 277 for (size_t i = 0; i < _N+_SE-_SS; ++i) 278 _data[i] += c._data[i]; 279 return *this; 280 } 281 282 LightVector& operator-=(const LightVector& c) 283 { 284 #pragma omp parallel for LOOPSCHED 285 for (size_t i = 0; i < _N+_SE-_SS; ++i) 286 _data[i] -= c._data[i]; 287 return *this; 288 } 289 290 LightVector& operator*=(const LightVector& c) 291 { 292 #pragma omp parallel for LOOPSCHED 293 for (size_t i = 0; i < _N+_SE-_SS; ++i) 294 _data[i] *= c._data[i]; 295 return *this; 296 } 297 298 LightVector& operator/=(const LightVector& c) 299 { 300 #pragma omp parallel for LOOPSCHED 301 for (size_t i = 0; i < _N+_SE-_SS; ++i) 302 _data[i] /= c._data[i]; 303 return *this; 304 } 305 306 LightVector operator+(const LightVector& c) const 307 { 308 LightVector thisCopy(*this); 309 return thisCopy += c; 310 } 311 312 LightVector operator-(const LightVector& c) const 313 { 314 LightVector thisCopy(*this); 315 return thisCopy -= c; 316 } 317 318 LightVector operator*(const LightVector& c) const 319 { 320 LightVector thisCopy(*this); 321 return thisCopy *= c; 322 } 323 324 LightVector operator/(const LightVector& c) const 325 { 326 LightVector thisCopy(*this); 327 return thisCopy /= c; 328 } 329 330 // other rhs types (NOTE: parameter are assed by value. The U-types are 331 // meant to be built-in or cheap to copy) 332 template <typename U> 333 LightVector& operator=(const U c) 334 { 335 #pragma omp parallel for LOOPSCHED 336 for (size_t i = 0; i < _N+_SE-_SS; ++i) 337 _data[i] = c; 338 return *this; 339 } 340 341 template <typename U> 342 LightVector& operator+=(const U c) 343 { 344 #pragma omp parallel for LOOPSCHED 345 for (size_t i = 0; i < _N+_SE-_SS; ++i) 346 _data[i] += c; 347 return *this; 348 } 349 350 template <typename U> 351 LightVector& operator-=(const U c) 352 { 353 #pragma omp parallel for LOOPSCHED 354 for (size_t i = 0; i < _N+_SE-_SS; ++i) 355 _data[i] -= c; 356 return *this; 357 } 358 359 template <typename U> 360 LightVector& operator*=(const U c) 361 { 362 #pragma omp parallel for LOOPSCHED 363 for (size_t i = 0; i < _N+_SE-_SS; ++i) 364 _data[i] *= c; 365 return *this; 366 } 367 368 template <typename U> 369 LightVector& operator/=(const U c) 370 { 371 #pragma omp parallel for LOOPSCHED 372 for (size_t i = 0; i < _N+_SE-_SS; ++i) 373 _data[i] /= c; 374 return *this; 375 } 376 377 template <typename U> 378 LightVector operator+(const U c) const 379 { 380 LightVector thisCopy(*this); 381 return thisCopy += c; 382 } 383 384 template <typename U> 385 LightVector operator-(const U c) const 386 { 387 LightVector thisCopy(*this); 388 return thisCopy -= c; 389 } 390 391 template <typename U> 392 LightVector operator*(const U c) const 393 { 394 LightVector thisCopy(*this); 395 return thisCopy *= c; 396 } 397 398 template <typename U> 399 LightVector operator/(const U c) const 400 { 401 LightVector thisCopy(*this); 402 return thisCopy /= c; 403 } 404 405 // output operator 406 friend std::ostream& operator<< <>(std::ostream& out, const LightVector& v); 407 408 inline void allocate(const size_t n) { _N=n; _alloc(n+(_SE-_SS)); _clear(); } 409 inline size_t size() const { return _N; } 410 inline size_t size_all() const { return _N+_SE-_SS; } 411 inline void clear() { _clear(); } 412 inline T& operator[](const size_t i) { assert(i<_N+_SE-_SS); return _data[i]; } 413 inline const T& operator[](const size_t i) const { assert(i<_N+_SE-_SS); return _data[i]; } 414 inline T& operator()(const int ix) { assert((ix-_SS >= 0)&&(ix-_SS < static_cast<int>(_N)+_SE-_SS)); return _data[ix-_SS]; } 415 inline const T& operator()(const int ix) const { assert((ix-_SS >= 0)&&(ix-_SS < static_cast<int>(_N)+_SE-_SS)); return _data[ix-_SS]; } 416 inline T* data() { return _data; } 417 inline T *data() const { return _data; } 418 }; 419 420 template <typename T, typename U, int _SS, int _SE> 421 LightVector<T,_SS,_SE> operator+(const U c, const LightVector<T,_SS,_SE>& v) 422 { 423 LightVector<T,_SS,_SE> vCopy(v); 424 return vCopy += c; 425 } 426 427 template <typename T, typename U, int _SS, int _SE> 428 LightVector<T,_SS,_SE> operator-(const U c, const LightVector<T,_SS,_SE>& v) 429 { 430 LightVector<T,_SS,_SE> vCopy(v); 431 return vCopy -= c; 432 } 433 434 template <typename T, typename U, int _SS, int _SE> 435 LightVector<T,_SS,_SE> operator*(const U c, const LightVector<T,_SS,_SE>& v) 436 { 437 LightVector<T,_SS,_SE> vCopy(v); 438 return vCopy *= c; 439 } 440 441 template <typename T, int _SS, int _SE> 442 inline std::ostream& operator<<(std::ostream& out, const LightVector<T,_SS,_SE>& v) 443 { 444 for (size_t j = 0; j < v.size(); ++j) 445 out << v[j] << std::endl; 446 return out; 447 } 448 449 #endif /* COMMON_H_4JKNXGVM */