1 /* 2 * This file is part of libcxxsupport. 3 * 4 * libcxxsupport is free software; you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation; either version 2 of the License, or 7 * (at your option) any later version. 8 * 9 * libcxxsupport is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with libcxxsupport; if not, write to the Free Software 16 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 17 */ 18 19 /* 20 * libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik 21 * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt 22 * (DLR). 23 */ 24 25 /*! \file arr.h 26 * Various high-performance array classes used by the Planck LevelS package. 27 * 28 * Copyright (C) 2002-2015 Max-Planck-Society 29 * \author Martin Reinecke 30 */ 31 32 #ifndef PLANCK_ARR_H 33 #define PLANCK_ARR_H 34 35 #include <algorithm> 36 #include <vector> 37 #include <cstdlib> 38 #include "alloc_utils.h" 39 #include "datatypes.h" 40 #include "math_utils.h" 41 42 /*! \defgroup arraygroup Array classes */ 43 /*! \{ */ 44 45 /*! View of a 1D array */ 46 template <typename T> class arr_ref 47 { 48 protected: 49 tsize s; 50 T *d; 51 52 public: 53 /*! Constructs an \a arr_ref of size \a s_, starting at \a d_. */ arr_ref(T * d_,tsize s_)54 arr_ref(T *d_, tsize s_) : s(s_),d(d_) {} 55 56 /*! Returns the current array size. */ size()57 tsize size() const { return s; } 58 59 /*! Writes \a val into every element of the array. */ fill(const T & val)60 void fill (const T &val) 61 { for (tsize m=0; m<s; ++m) d[m]=val; } 62 63 /*! Returns a reference to element \a n */ 64 template<typename T2> T &operator[] (T2 n) {return d[n];} 65 /*! Returns a constant reference to element \a n */ 66 template<typename T2> const T &operator[] (T2 n) const {return d[n];} 67 68 /*! Returns a pointer to the first array element, or NULL if the array 69 is zero-sized. */ begin()70 T *begin() { return d; } 71 /*! Returns a pointer to the one-past-last array element, or NULL if the 72 array is zero-sized. */ end()73 T *end() { return d+s; } 74 /*! Returns a constant pointer to the first array element, or NULL if the 75 array is zero-sized. */ begin()76 const T *begin() const { return d; } 77 /*! Returns a constant pointer to the one-past-last array element, or NULL 78 if the array is zero-sized. */ end()79 const T *end() const { return d+s; } 80 81 /*! Copies all array elements to \a ptr. */ copyToPtr(T * ptr)82 template<typename T2> void copyToPtr (T *ptr) const 83 { for (tsize m=0; m<s; ++m) ptr[m]=d[m]; } 84 85 /*! Sorts the elements in the array, in ascending order. */ sort()86 void sort() 87 { std::sort (d,d+s); } 88 89 /*! Sorts the elements in the array, such that \a comp(d[i],d[j])==true 90 for \a i<j. */ sort(Comp comp)91 template<typename Comp> void sort(Comp comp) 92 { std::sort (d,d+s,comp); } 93 94 /*! Helper function for linear interpolation (or extrapolation). 95 \a idx and \a val are computed such that 96 \a val=d[idx]+frac*(d[idx+1]-d[idx]). If \a val<d[0], \a frac will be 97 negative, if \a val>d[s-1], frac will be larger than 1. In all other 98 cases \a 0<=frac<=1. 99 100 The array must be ordered in ascending order; no two values may be 101 equal. */ interpol_helper(const T & val,tsize & idx,double & frac)102 void interpol_helper (const T &val, tsize &idx, double &frac) const 103 { ::interpol_helper (d, d+s, val, idx, frac); } 104 105 /*! Helper function for linear interpolation (or extrapolation). 106 \a idx and \a val are computed such that 107 \a val=d[idx]+frac*(d[idx+1]-d[idx]). If \a comp(val,d[0])==true, 108 \a frac will be negative, if \a comp(val,d[s-1])==false, frac will be 109 larger than 1. In all other cases \a 0<=frac<=1. 110 111 The array must be ordered such that \a comp(d[i],d[j])==true 112 for \a i<j; no two values may be equal. */ interpol_helper(const T & val,Comp comp,tsize & idx,double & frac)113 template<typename Comp> void interpol_helper (const T &val, Comp comp, 114 tsize &idx, double &frac) const 115 { ::interpol_helper (d, d+s, val, comp, idx, frac); } 116 117 /*! Returns the minimum and maximum entry in \a minv and \a maxv, 118 respectively. Throws an exception if the array is zero-sized. */ minmax(T & minv,T & maxv)119 void minmax (T &minv, T &maxv) const 120 { 121 planck_assert(s>0,"trying to find min and max of a zero-sized array"); 122 minv=maxv=d[0]; 123 for (tsize m=1; m<s; ++m) 124 { 125 if (d[m]<minv) minv=d[m]; 126 else if (d[m]>maxv) maxv=d[m]; 127 } 128 } 129 130 /*! Returns \a true, if \a val is found in the array, else \a false. */ contains(const T & val)131 bool contains (const T &val) const 132 { 133 for (tsize m=0; m<s; ++m) 134 if (d[m]==val) return true; 135 return false; 136 } 137 138 /*! Returns the index of the first occurrence of \a val in the array. 139 If it is not found, an exception is thrown. */ find(const T & val)140 tsize find (const T &val) const 141 { 142 for (tsize m=0; m<s; ++m) 143 if (d[m]==val) return m; 144 planck_fail ("entry not found in array"); 145 } 146 147 /*! Returns \a true if the array has the same size as \a other and all 148 elements of both arrays are equal, else \a false. */ contentsEqual(const arr_ref & other)149 bool contentsEqual(const arr_ref &other) const 150 { 151 if (s!=other.s) return false; 152 for (tsize i=0; i<s; ++i) 153 if (d[i]!=other.d[i]) return false; 154 return true; 155 } 156 }; 157 158 /*! An array whose size is known at compile time. Very useful for storing 159 small arrays on the stack, without need for \a new and \a delete(). */ 160 template <typename T, tsize sz> class fix_arr 161 { 162 private: 163 T d[sz]; 164 165 public: 166 /*! Returns the size of the array. */ size()167 tsize size() const { return sz; } 168 169 /*! Returns a reference to element \a n */ 170 template<typename T2> T &operator[] (T2 n) {return d[n];} 171 /*! Returns a constant reference to element \a n */ 172 template<typename T2> const T &operator[] (T2 n) const {return d[n];} 173 }; 174 175 176 /*! One-dimensional array type, with selectable storage management. */ 177 template <typename T, typename stm> class arrT: public arr_ref<T> 178 { 179 private: 180 bool own; 181 reset()182 void reset() 183 { this->d=0; this->s=0; own=true; } 184 185 public: 186 /*! Creates a zero-sized array. */ arrT()187 arrT() : arr_ref<T>(0,0), own(true) {} 188 /*! Creates an array with \a sz entries. */ arrT(tsize sz)189 explicit arrT(tsize sz) : arr_ref<T>(stm::alloc(sz),sz), own(true) {} 190 /*! Creates an array with \a sz entries, and initializes them with 191 \a inival. */ arrT(tsize sz,const T & inival)192 arrT(tsize sz, const T &inival) : arr_ref<T>(stm::alloc(sz),sz), own(true) 193 { this->fill(inival); } 194 /*! Creates an array with \a sz entries, which uses the memory pointed 195 to by \a ptr. 196 \note \a ptr will <i>not</i> be deallocated by the destructor. 197 \warning Only use this if you REALLY know what you are doing. 198 In particular, this is only safely usable if 199 <ul> 200 <li>\a T is a POD type</li> 201 <li>\a ptr survives during the lifetime of the array object</li> 202 <li>\a ptr is not subject to garbage collection</li> 203 </ul> 204 Other restrictions may apply. You have been warned. */ arrT(T * ptr,tsize sz)205 arrT (T *ptr, tsize sz): arr_ref<T>(ptr,sz), own(false) {} 206 /*! Creates an array which is a copy of \a orig. The data in \a orig 207 is duplicated. */ arrT(const arrT & orig)208 arrT (const arrT &orig): arr_ref<T>(stm::alloc(orig.s),orig.s), own(true) 209 { for (tsize m=0; m<this->s; ++m) this->d[m] = orig.d[m]; } 210 /*! Frees the memory allocated by the object. */ ~arrT()211 ~arrT() { if (own) stm::dealloc(this->d); } 212 213 /*! Allocates space for \a sz elements. The content of the array is 214 undefined on exit. \a sz can be 0. If \a sz is the 215 same as the current size, no reallocation is performed. */ alloc(tsize sz)216 void alloc (tsize sz) 217 { 218 if (sz==this->s) return; 219 if (own) stm::dealloc(this->d); 220 this->s = sz; 221 this->d = stm::alloc(sz); 222 own = true; 223 } 224 /*! Allocates space for \a sz elements. If \a sz is the 225 same as the current size, no reallocation is performed. 226 All elements are set to \a inival. */ allocAndFill(tsize sz,const T & inival)227 void allocAndFill (tsize sz, const T &inival) 228 { alloc(sz); this->fill(inival); } 229 /*! Deallocates the memory held by the array, and sets the array size 230 to 0. */ dealloc()231 void dealloc() {if (own) stm::dealloc(this->d); reset();} 232 /*! Resizes the array to hold \a sz elements. The existing content of the 233 array is copied over to the new array to the extent possible. 234 \a sz can be 0. If \a sz is the same as the current size, no 235 reallocation is performed. */ resize(tsize sz)236 void resize (tsize sz) 237 { 238 using namespace std; 239 if (sz==this->s) return; 240 T *tmp = stm::alloc(sz); 241 for (tsize m=0; m<min(sz,this->s); ++m) 242 tmp[m]=this->d[m]; 243 if (own) stm::dealloc(this->d); 244 this->s = sz; 245 this->d = tmp; 246 own = true; 247 } 248 249 /*! Changes the array to be a copy of \a orig. */ 250 arrT &operator= (const arrT &orig) 251 { 252 if (this==&orig) return *this; 253 alloc (orig.s); 254 for (tsize m=0; m<this->s; ++m) this->d[m] = orig.d[m]; 255 return *this; 256 } 257 258 /*! Changes the array to be a copy of the std::vector \a orig. */ copyFrom(const std::vector<T2> & orig)259 template<typename T2> void copyFrom (const std::vector<T2> &orig) 260 { 261 alloc (orig.size()); 262 for (tsize m=0; m<this->s; ++m) this->d[m] = orig[m]; 263 } 264 /*! Changes the std::vector \a vec to be a copy of the object. */ copyTo(std::vector<T2> & vec)265 template<typename T2> void copyTo (std::vector<T2> &vec) const 266 { 267 vec.clear(); vec.reserve(this->s); 268 for (tsize m=0; m<this->s; ++m) vec.push_back(this->d[m]); 269 } 270 271 /*! Reserves space for \a sz elements, then copies \a sz elements 272 from \a ptr into the array. */ copyFromPtr(const T2 * ptr,tsize sz)273 template<typename T2> void copyFromPtr (const T2 *ptr, tsize sz) 274 { 275 alloc(sz); 276 for (tsize m=0; m<this->s; ++m) this->d[m]=ptr[m]; 277 } 278 279 /*! Assigns the contents and size of \a other to the array. 280 \note On exit, \a other is zero-sized! */ transfer(arrT & other)281 void transfer (arrT &other) 282 { 283 if (own) stm::dealloc(this->d); 284 this->d=other.d; 285 this->s=other.s; 286 own=other.own; 287 other.reset(); 288 } 289 /*! Swaps contents and size with \a other. */ swap(arrT & other)290 void swap (arrT &other) 291 { 292 std::swap(this->d,other.d); 293 std::swap(this->s,other.s); 294 std::swap(own,other.own); 295 } 296 }; 297 298 /*! One-dimensional array type. */ 299 template <typename T> 300 class arr: public arrT<T,normalAlloc__<T> > 301 { 302 public: 303 /*! Creates a zero-sized array. */ arr()304 arr() : arrT<T,normalAlloc__<T> >() {} 305 /*! Creates an array with \a sz entries. */ arr(tsize sz)306 explicit arr(tsize sz) : arrT<T,normalAlloc__<T> >(sz) {} 307 /*! Creates an array with \a sz entries, and initializes them with 308 \a inival. */ arr(tsize sz,const T & inival)309 arr(tsize sz, const T &inival) : arrT<T,normalAlloc__<T> >(sz,inival) {} 310 /*! Creates an array with \a sz entries, which uses the memory pointed 311 to by \a ptr. 312 \note \a ptr will <i>not</i> be deallocated by the destructor. 313 \warning Only use this if you REALLY know what you are doing. 314 In particular, this is only safely usable if 315 <ul> 316 <li>\a T is a POD type</li> 317 <li>\a ptr survives during the lifetime of the array object</li> 318 <li>\a ptr is not subject to garbage collection</li> 319 </ul> 320 Other restrictions may apply. You have been warned. */ arr(T * ptr,tsize sz)321 arr (T *ptr, tsize sz): arrT<T,normalAlloc__<T> >(ptr,sz) {} 322 /*! Creates an array which is a copy of \a orig. The data in \a orig 323 is duplicated. */ arr(const arr & orig)324 arr (const arr &orig): arrT<T,normalAlloc__<T> >(orig) {} 325 }; 326 327 /*! One-dimensional array type, with selectable storage alignment. */ 328 template <typename T, int align> 329 class arr_align: public arrT<T,alignAlloc__<T,align> > 330 { 331 public: 332 /*! Creates a zero-sized array. */ arr_align()333 arr_align() : arrT<T,alignAlloc__<T,align> >() {} 334 /*! Creates an array with \a sz entries. */ arr_align(tsize sz)335 explicit arr_align(tsize sz) : arrT<T,alignAlloc__<T,align> >(sz) {} 336 /*! Creates an array with \a sz entries, and initializes them with 337 \a inival. */ arr_align(tsize sz,const T & inival)338 arr_align(tsize sz, const T &inival) 339 : arrT<T,alignAlloc__<T,align> >(sz,inival) {} 340 }; 341 342 343 /*! Two-dimensional array type, with selectable storage management. 344 The storage ordering is the same as in C. 345 An entry is located by address arithmetic, not by double dereferencing. 346 The indices start at zero. */ 347 template <typename T, typename storageManager> class arr2T 348 { 349 private: 350 tsize s1, s2; 351 arrT<T, storageManager> d; 352 353 public: 354 /*! Creates a zero-sized array. */ arr2T()355 arr2T() : s1(0), s2(0) {} 356 /*! Creates an array with the dimensions \a sz1 and \a sz2. */ arr2T(tsize sz1,tsize sz2)357 arr2T(tsize sz1, tsize sz2) 358 : s1(sz1), s2(sz2), d(s1*s2) {} 359 /*! Creates an array with the dimensions \a sz1 and \a sz2 360 and initializes them with \a inival. */ 361 /*! Creates an array with the dimensions \a sz1 and \a sz2 from existing 362 pointer. */ arr2T(T * p,tsize sz1,tsize sz2)363 arr2T(T* p, tsize sz1, tsize sz2) 364 : s1(sz1), s2(sz2), d(p, s1*s2) {} arr2T(tsize sz1,tsize sz2,const T & inival)365 arr2T(tsize sz1, tsize sz2, const T &inival) 366 : s1(sz1), s2(sz2), d (s1*s2) 367 { fill(inival); } 368 /*! Creates the array as a copy of \a orig. */ arr2T(const arr2T & orig)369 arr2T(const arr2T &orig) 370 : s1(orig.s1), s2(orig.s2), d(orig.d) {} 371 /*! Frees the memory associated with the array. */ ~arr2T()372 ~arr2T() {} 373 374 /*! Returns the first array dimension. */ size1()375 tsize size1() const { return s1; } 376 /*! Returns the second array dimension. */ size2()377 tsize size2() const { return s2; } 378 /*! Returns the total array size, i.e. the product of both dimensions. */ size()379 tsize size () const { return s1*s2; } 380 381 /*! Allocates space for an array with \a sz1*sz2 elements. 382 The content of the array is undefined on exit. 383 \a sz1 or \a sz2 can be 0. If \a sz1*sz2 is the same as the 384 currently allocated space, no reallocation is performed. */ alloc(tsize sz1,tsize sz2)385 void alloc (tsize sz1, tsize sz2) 386 { 387 if (sz1*sz2 != d.size()) 388 d.alloc(sz1*sz2); 389 s1=sz1; s2=sz2; 390 } 391 /*! Allocates space for an array with \a sz1*sz2 elements. 392 All elements are set to \a inival. 393 \a sz1 or \a sz2 can be 0. If \a sz1*sz2 is the same as the 394 currently allocated space, no reallocation is performed. */ allocAndFill(tsize sz1,tsize sz2,const T & inival)395 void allocAndFill (tsize sz1, tsize sz2, const T &inival) 396 { alloc(sz1,sz2); fill(inival); } 397 /*! Allocates space for an array with \a sz1*sz2 elements. 398 The content of the array is undefined on exit. 399 \a sz1 or \a sz2 can be 0. If \a sz1*sz2 is smaller than the 400 currently allocated space, no reallocation is performed. */ fast_alloc(tsize sz1,tsize sz2)401 void fast_alloc (tsize sz1, tsize sz2) 402 { 403 if (sz1*sz2<=d.size()) 404 { s1=sz1; s2=sz2; } 405 else 406 alloc(sz1,sz2); 407 } 408 /*! Deallocates the space and makes the array zero-sized. */ dealloc()409 void dealloc () {d.dealloc(); s1=0; s2=0;} 410 411 /*! Sets all array elements to \a val. */ fill(const T & val)412 void fill (const T &val) 413 { for (tsize m=0; m<s1*s2; ++m) d[m]=val; } 414 415 /*! Multiplies all array elements by \a val. */ scale(const T & val)416 void scale (const T &val) 417 { for (tsize m=0; m<s1*s2; ++m) d[m]*=val; } 418 419 /*! Changes the array to be a copy of \a orig. */ 420 arr2T &operator= (const arr2T &orig) 421 { 422 if (this==&orig) return *this; 423 alloc (orig.s1, orig.s2); 424 d = orig.d; 425 return *this; 426 } 427 428 /*! Returns a pointer to the beginning of slice \a n. */ 429 template<typename T2> T *operator[] (T2 n) {return &d[n*s2];} 430 /*! Returns a constant pointer to the beginning of slice \a n. */ 431 template<typename T2> const T *operator[] (T2 n) const {return &d[n*s2];} 432 433 /*! Returns a reference to the element with the indices \a n1 and \a n2. */ operator()434 template<typename T2, typename T3> T &operator() (T2 n1, T3 n2) 435 {return d[n1*s2 + n2];} 436 /*! Returns a constant reference to the element with the indices 437 \a n1 and \a n2. */ operator()438 template<typename T2, typename T3> const T &operator() (T2 n1, T3 n2) const 439 {return d[n1*s2 + n2];} 440 441 /*! Returns the minimum and maximum entry in \a minv and \a maxv, 442 respectively. Throws an exception if the array is zero-sized. */ minmax(T & minv,T & maxv)443 void minmax (T &minv, T &maxv) const 444 { 445 planck_assert(s1*s2>0, 446 "trying to find min and max of a zero-sized array"); 447 minv=maxv=d[0]; 448 for (tsize m=1; m<s1*s2; ++m) 449 { 450 if (d[m]<minv) minv=d[m]; 451 if (d[m]>maxv) maxv=d[m]; 452 } 453 } 454 455 /*! Swaps contents and sizes with \a other. */ swap(arr2T & other)456 void swap (arr2T &other) 457 { 458 d.swap(other.d); 459 std::swap(s1,other.s1); 460 std::swap(s2,other.s2); 461 } 462 463 /*! Returns \c true if the array and \a other have the same dimensions, 464 else \c false. */ conformable(const arr2T<T2,T3> & other)465 template<typename T2, typename T3> bool conformable 466 (const arr2T<T2,T3> &other) const 467 { return (other.size1()==s1) && (other.size2()==s2); } 468 }; 469 470 /*! Two-dimensional array type. The storage ordering is the same as in C. 471 An entry is located by address arithmetic, not by double dereferencing. 472 The indices start at zero. */ 473 template <typename T> 474 class arr2: public arr2T<T,normalAlloc__<T> > 475 { 476 public: 477 /*! Creates a zero-sized array. */ arr2()478 arr2() : arr2T<T,normalAlloc__<T> > () {} 479 /*! Creates an array with the dimensions \a sz1 and \a sz2. */ arr2(tsize sz1,tsize sz2)480 arr2(tsize sz1, tsize sz2) : arr2T<T,normalAlloc__<T> > (sz1,sz2) {} 481 /*! Creates an array with the dimensions \a sz1 and \a sz2 from existing 482 pointer. */ arr2(T * p,tsize sz1,tsize sz2)483 arr2(T* p, tsize sz1, tsize sz2) : arr2T<T,normalAlloc__<T> > (p,sz1,sz2) {} 484 /*! Creates an array with the dimensions \a sz1 and \a sz2 485 and initializes them with \a inival. */ arr2(tsize sz1,tsize sz2,const T & inival)486 arr2(tsize sz1, tsize sz2, const T &inival) 487 : arr2T<T,normalAlloc__<T> > (sz1,sz2,inival) {} 488 }; 489 490 /*! Two-dimensional array type, with selectable storage alignment. 491 The storage ordering is the same as in C. 492 An entry is located by address arithmetic, not by double dereferencing. 493 The indices start at zero. */ 494 template <typename T, int align> 495 class arr2_align: public arr2T<T,alignAlloc__<T,align> > 496 { 497 public: 498 /*! Creates a zero-sized array. */ arr2_align()499 arr2_align() : arr2T<T,alignAlloc__<T,align> > () {} 500 /*! Creates an array with the dimensions \a sz1 and \a sz2. */ arr2_align(tsize sz1,tsize sz2)501 arr2_align(tsize sz1, tsize sz2) 502 : arr2T<T,alignAlloc__<T,align> > (sz1,sz2) {} 503 /*! Creates an array with the dimensions \a sz1 and \a sz2 504 and initializes them with \a inival. */ arr2_align(tsize sz1,tsize sz2,const T & inival)505 arr2_align(tsize sz1, tsize sz2, const T &inival) 506 : arr2T<T,alignAlloc__<T,align> > (sz1,sz2,inival) {} 507 }; 508 509 /*! Two-dimensional array type. An entry is located by double dereferencing, 510 i.e. via an array of pointers. The indices start at zero. */ 511 template <typename T> class arr2b 512 { 513 private: 514 tsize s1, s2; 515 arr<T> d; 516 arr<T *> d1; 517 fill_d1()518 void fill_d1() 519 { for (tsize m=0; m<s1; ++m) d1[m] = &d[m*s2]; } 520 521 public: 522 /*! Creates a zero-sized array. */ arr2b()523 arr2b() : s1(0), s2(0), d(0), d1(0) {} 524 /*! Creates an array with the dimensions \a sz1 and \a sz2. */ arr2b(tsize sz1,tsize sz2)525 arr2b(tsize sz1, tsize sz2) 526 : s1(sz1), s2(sz2), d(s1*s2), d1(s1) 527 { fill_d1(); } 528 /*! Creates the array as a copy of \a orig. */ arr2b(const arr2b & orig)529 arr2b(const arr2b &orig) 530 : s1(orig.s1), s2(orig.s2), d(orig.d), d1(s1) 531 { fill_d1(); } 532 /*! Frees the memory associated with the array. */ ~arr2b()533 ~arr2b() {} 534 535 /*! Returns the first array dimension. */ size1()536 tsize size1() const { return s1; } 537 /*! Returns the second array dimension. */ size2()538 tsize size2() const { return s2; } 539 /*! Returns the total array size, i.e. the product of both dimensions. */ size()540 tsize size () const { return s1*s2; } 541 542 /*! Allocates space for an array with \a sz1*sz2 elements. 543 The content of the array is undefined on exit. */ alloc(tsize sz1,tsize sz2)544 void alloc (tsize sz1, tsize sz2) 545 { 546 if ((s1==sz1) && (s2==sz2)) return; 547 s1=sz1; s2=sz2; 548 d.alloc(s1*s2); 549 d1.alloc(s1); 550 fill_d1(); 551 } 552 /*! Deallocates the space and makes the array zero-sized. */ dealloc()553 void dealloc () {d.dealloc(); d1.dealloc(); s1=0; s2=0;} 554 555 /*! Sets all array elements to \a val. */ fill(const T & val)556 void fill (const T &val) 557 { d.fill(val); } 558 559 /*! Changes the array to be a copy of \a orig. */ 560 arr2b &operator= (const arr2b &orig) 561 { 562 if (this==&orig) return *this; 563 alloc (orig.s1, orig.s2); 564 for (tsize m=0; m<s1*s2; ++m) d[m] = orig.d[m]; 565 return *this; 566 } 567 568 /*! Returns a pointer to the beginning of slice \a n. */ 569 template<typename T2> T *operator[] (T2 n) {return d1[n];} 570 /*! Returns a constant pointer to the beginning of slice \a n. */ 571 template<typename T2> const T *operator[] (T2 n) const {return d1[n];} 572 573 /*! Returns a pointer to the beginning of the pointer array. */ p0()574 T **p0() {return &d1[0];} 575 }; 576 577 578 /*! Three-dimensional array type. The storage ordering is the same as in C. 579 An entry is located by address arithmetic, not by multiple dereferencing. 580 The indices start at zero. */ 581 template <typename T> class arr3 582 { 583 private: 584 tsize s1, s2, s3, s2s3; 585 arr<T> d; 586 587 public: 588 /*! Creates a zero-sized array. */ arr3()589 arr3() : s1(0), s2(0), s3(0), s2s3(0), d(0) {} 590 /*! Creates an array with the dimensions \a sz1, \a sz2 and \a sz3. */ arr3(tsize sz1,tsize sz2,tsize sz3)591 arr3(tsize sz1, tsize sz2, tsize sz3) 592 : s1(sz1), s2(sz2), s3(sz3), s2s3(s2*s3), d(s1*s2*s3) {} 593 /*! Creates the array as a copy of \a orig. */ arr3(const arr3 & orig)594 arr3(const arr3 &orig) 595 : s1(orig.s1), s2(orig.s2), s3(orig.s3), s2s3(orig.s2s3), d(orig.d) {} 596 /*! Frees the memory associated with the array. */ ~arr3()597 ~arr3() {} 598 599 /*! Returns the first array dimension. */ size1()600 tsize size1() const { return s1; } 601 /*! Returns the second array dimension. */ size2()602 tsize size2() const { return s2; } 603 /*! Returns the third array dimension. */ size3()604 tsize size3() const { return s3; } 605 /*! Returns the total array size, i.e. the product of all dimensions. */ size()606 tsize size () const { return s1*s2*s3; } 607 608 /*! Allocates space for an array with \a sz1*sz2*sz3 elements. 609 The content of the array is undefined on exit. */ alloc(tsize sz1,tsize sz2,tsize sz3)610 void alloc (tsize sz1, tsize sz2, tsize sz3) 611 { 612 d.alloc(sz1*sz2*sz3); 613 s1=sz1; s2=sz2; s3=sz3; s2s3=s2*s3; 614 } 615 /*! Deallocates the space and makes the array zero-sized. */ dealloc()616 void dealloc () {d.dealloc(); s1=0; s2=0; s3=0; s2s3=0;} 617 618 /*! Sets all array elements to \a val. */ fill(const T & val)619 void fill (const T &val) 620 { d.fill(val); } 621 622 /*! Changes the array to be a copy of \a orig. */ 623 arr3 &operator= (const arr3 &orig) 624 { 625 if (this==&orig) return *this; 626 alloc (orig.s1, orig.s2, orig.s3); 627 d = orig.d; 628 return *this; 629 } 630 631 /*! Returns a reference to the element with the indices 632 \a n1, \a n2 and \a n3. */ operator()633 template<typename T2, typename T3, typename T4> T &operator() 634 (T2 n1, T3 n2, T4 n3) 635 {return d[n1*s2s3 + n2*s3 + n3];} 636 /*! Returns a constant reference to the element with the indices 637 \a n1, \a n2 and \a n3. */ operator()638 template<typename T2, typename T3, typename T4> const T &operator() 639 (T2 n1, T3 n2, T4 n3) const 640 {return d[n1*s2s3 + n2*s3 + n3];} 641 642 /*! Swaps contents and sizes with \a other. */ swap(arr3 & other)643 void swap (arr3 &other) 644 { 645 d.swap(other.d); 646 std::swap(s1,other.s1); 647 std::swap(s2,other.s2); 648 std::swap(s3,other.s3); 649 std::swap(s2s3,other.s2s3); 650 } 651 652 /*! Returns \c true if the array and \a other have the same dimensions, 653 else \c false. */ conformable(const arr3<T2> & other)654 template<typename T2> bool conformable (const arr3<T2> &other) const 655 { return (other.size1()==s1)&&(other.size2()==s2)&&(other.size3()==s3); } 656 }; 657 658 /*! \} */ 659 660 #endif 661