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