Marray
marray.hxx
1 
2 
3 
4 
5 
6 #pragma once
7 #ifndef ANDRES_MARRAY_HXX
8 #define ANDRES_MARRAY_HXX
9 
10 #include <cassert>
11 #include <cstddef>
12 #include <stdexcept>
13 #include <limits>
14 #include <string>
15 #include <sstream>
16 #include <cstring> // memcpy
17 #include <iterator>
18 #include <vector>
19 #include <set>
20 #include <iostream>
21 #include <memory> // allocator
22 #include <numeric> // accumulate
23 #include <functional> // std::multiplies
24 #ifdef HAVE_CPP11_INITIALIZER_LISTS
25  #include <initializer_list>
26 #endif
27 
29 namespace andres {
30 
31 enum StringStyle {TableStyle, MatrixStyle};
32 enum CoordinateOrder {FirstMajorOrder, LastMajorOrder};
33 struct InitializationSkipping { };
34 
35 static const bool Const = true;
36 static const bool Mutable = false;
37 static const CoordinateOrder defaultOrder = LastMajorOrder;
38 static const InitializationSkipping SkipInitialization = InitializationSkipping();
39 
40 template<class E, class T>
41  class ViewExpression;
42 // \cond suppress_doxygen
43 template<class E, class T, class UnaryFunctor>
44  class UnaryViewExpression;
45 template<class E1, class T1, class E2, class T2, class BinaryFunctor>
46  class BinaryViewExpression;
47 template<class E, class T, class S, class BinaryFunctor>
48  class BinaryViewExpressionScalarFirst;
49 template<class E, class T, class S, class BinaryFunctor>
50  class BinaryViewExpressionScalarSecond;
51 // \endcond suppress_doxygen
52 template<class T, bool isConst = false, class A = std::allocator<std::size_t> >
53  class View;
54 #ifdef HAVE_CPP11_TEMPLATE_ALIASES
55  template<class T, class A> using ConstView = View<T, true, A>;
56 #endif
57 template<class T, bool isConst, class A = std::allocator<std::size_t> >
58  class Iterator;
59 template<class T, class A = std::allocator<std::size_t> > class Marray;
60 
61 // assertion testing
62 #ifdef NDEBUG
63  const bool MARRAY_NO_DEBUG = true;
64  const bool MARRAY_NO_ARG_TEST = true;
65 #else
66  const bool MARRAY_NO_DEBUG = false;
67  const bool MARRAY_NO_ARG_TEST = false;
68 #endif
69 
70 // \cond suppress_doxygen
71 namespace marray_detail {
72  // meta-programming
73  template <bool PREDICATE, class TRUECASE, class FALSECASE>
74  struct IfBool;
75  template <class TRUECASE, class FALSECASE>
76  struct IfBool<true, TRUECASE, FALSECASE>
77  { typedef TRUECASE type; };
78  template <class TRUECASE, class FALSECASE>
79  struct IfBool<false, TRUECASE, FALSECASE>
80  { typedef FALSECASE type; };
81 
82  template <class T1, class T2>
83  struct IsEqual
84  { static const bool type = false; };
85  template <class T>
86  struct IsEqual<T, T>
87  { static const bool type = true; };
88 
89  template<class T> struct TypeTraits
90  { static const unsigned char position = 255; };
91  template<> struct TypeTraits<char>
92  { static const unsigned char position = 0; };
93  template<> struct TypeTraits<unsigned char>
94  { static const unsigned char position = 1; };
95  template<> struct TypeTraits<short>
96  { static const unsigned char position = 2; };
97  template<> struct TypeTraits<unsigned short>
98  { static const unsigned char position = 3; };
99  template<> struct TypeTraits<int>
100  { static const unsigned char position = 4; };
101  template<> struct TypeTraits<unsigned int>
102  { static const unsigned char position = 5; };
103  template<> struct TypeTraits<long>
104  { static const unsigned char position = 6; };
105  template<> struct TypeTraits<unsigned long>
106  { static const unsigned char position = 7; };
107  template<> struct TypeTraits<float>
108  { static const unsigned char position = 8; };
109  template<> struct TypeTraits<double>
110  { static const unsigned char position = 9; };
111  template<> struct TypeTraits<long double>
112  { static const unsigned char position = 10; };
113  template<class A, class B> struct PromoteType
114  { typedef typename IfBool<TypeTraits<A>::position >= TypeTraits<B>::position, A, B>::type type; };
115 
116  // assertion testing
117  template<class A> inline void Assert(A assertion) {
118  if(!assertion) throw std::runtime_error("Assertion failed.");
119  }
120 
121  // geometry of views
122  template<class A = std::allocator<std::size_t> > class Geometry;
123  template<class ShapeIterator, class StridesIterator>
124  inline void stridesFromShape(ShapeIterator, ShapeIterator,
125  StridesIterator, const CoordinateOrder& = defaultOrder);
126 
127  // operations on entries of views
128  template<class Functor, class T, class A>
129  inline void operate(View<T, false, A>&, Functor);
130  template<class Functor, class T, class A>
131  inline void operate(View<T, false, A>&, const T&, Functor);
132  template<class Functor, class T1, class T2, bool isConst, class A>
133  inline void operate(View<T1, false, A>&, const View<T2, isConst, A>&, Functor);
134  template<class Functor, class T1, class A, class E, class T2>
135  inline void operate(View<T1, false, A>& v, const ViewExpression<E, T2>& expression, Functor f);
136 
137  // helper classes
138  template<unsigned short N, class Functor, class T, class A>
139  struct OperateHelperUnary;
140  template<unsigned short N, class Functor, class T1, class T2, class A>
141  struct OperateHelperBinaryScalar;
142  template<unsigned short N, class Functor, class T1, class T2,
143  bool isConst, class A1, class A2>
144  struct OperateHelperBinary;
145  template<bool isConstTo, class TFrom, class TTo, class AFrom, class ATo>
146  struct AssignmentOperatorHelper;
147  template<bool isIntegral>
148  struct AccessOperatorHelper;
149 
150  // unary in-place functors
151  template<class T>
152  struct Negative { void operator()(T& x) { x = -x; } };
153  template<class T>
154  struct PrefixIncrement { void operator()(T& x) { ++x; } };
155  template<class T>
156  struct PostfixIncrement { void operator()(T& x) { x++; } };
157  template<class T>
158  struct PrefixDecrement { void operator()(T& x) { --x; } };
159  template<class T>
160  struct PostfixDecrement { void operator()(T& x) { x--; } };
161 
162  // binary in-place functors
163  template<class T1, class T2>
164  struct Assign { void operator()(T1& x, const T2& y) { x = y; } };
165  template<class T1, class T2>
166  struct PlusEqual { void operator()(T1& x, const T2& y) { x += y; } };
167  template<class T1, class T2>
168  struct MinusEqual { void operator()(T1& x, const T2& y) { x -= y; } };
169  template<class T1, class T2>
170  struct TimesEqual { void operator()(T1& x, const T2& y) { x *= y; } };
171  template<class T1, class T2>
172  struct DividedByEqual { void operator()(T1& x, const T2& y) { x /= y; } };
173 
174  // unary functors
175  template<class T>
176  struct Negate { T operator()(const T& x) const { return -x; } };
177 
178  // binary functors
179  template<class T1, class T2, class U>
180  struct Plus { U operator()(const T1& x, const T2& y) const { return x + y; } };
181  template<class T1, class T2, class U>
182  struct Minus { U operator()(const T1& x, const T2& y) const { return x - y; } };
183  template<class T1, class T2, class U>
184  struct Times { U operator()(const T1& x, const T2& y) const { return x * y; } };
185  template<class T1, class T2, class U>
186  struct DividedBy { U operator()(const T1& x, const T2& y) const { return x / y; } };
187 }
188 // \endcond suppress_doxygen
189 
208 template<class T, bool isConst, class A>
209 class View
210 : public ViewExpression<View<T, isConst, A>, T>
211 {
212 public:
213  typedef T value_type;
214  typedef typename marray_detail::IfBool<isConst, const T*, T*>::type pointer;
215  typedef const T* const_pointer;
216  typedef typename marray_detail::IfBool<isConst, const T&, T&>::type reference;
217  typedef const T& const_reference;
218  typedef Iterator<T, isConst, A> iterator;
219  typedef Iterator<T, true, A> const_iterator;
220  typedef std::reverse_iterator<iterator> reverse_iterator;
221  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
222  typedef ViewExpression<View<T, isConst, A>, T> base;
223  typedef typename A::template rebind<value_type>::other allocator_type;
224 
225  // construction
226  View(const allocator_type& = allocator_type());
227  View(pointer, const allocator_type& = allocator_type());
228  View(const View<T, false, A>&);
229  template<class ShapeIterator>
230  View(ShapeIterator, ShapeIterator, pointer,
231  const CoordinateOrder& = defaultOrder,
232  const CoordinateOrder& = defaultOrder,
233  const allocator_type& = allocator_type());
234  template<class ShapeIterator, class StrideIterator>
235  View(ShapeIterator, ShapeIterator, StrideIterator,
236  pointer, const CoordinateOrder&,
237  const allocator_type& = allocator_type());
238  #ifdef HAVE_CPP11_INITIALIZER_LISTS
239  View(std::initializer_list<std::size_t>, pointer,
240  const CoordinateOrder& = defaultOrder,
241  const CoordinateOrder& = defaultOrder,
242  const allocator_type& = allocator_type());
243  View(std::initializer_list<std::size_t>, std::initializer_list<std::size_t>,
244  pointer, const CoordinateOrder&,
245  const allocator_type& = allocator_type());
246  #endif
247 
248  // assignment
249  View<T, isConst, A>& operator=(const T&);
250  View<T, isConst, A>& operator=(const View<T, true, A>&); // over-write default
251  View<T, isConst, A>& operator=(const View<T, false, A>&); // over-write default
252  template<class TLocal, bool isConstLocal, class ALocal>
253  View<T, isConst, A>& operator=(const View<TLocal, isConstLocal, ALocal>&);
254  template<class E, class Te>
255  View<T, isConst, A>&
256  operator=(const ViewExpression<E, Te>&);
257 
258  void assign(const allocator_type& = allocator_type());
259  template<class ShapeIterator>
260  void assign(ShapeIterator, ShapeIterator, pointer,
261  const CoordinateOrder& = defaultOrder,
262  const CoordinateOrder& = defaultOrder,
263  const allocator_type& = allocator_type());
264  template<class ShapeIterator, class StrideIterator>
265  void assign(ShapeIterator, ShapeIterator, StrideIterator,
266  pointer, const CoordinateOrder&,
267  const allocator_type& = allocator_type());
268  #ifdef HAVE_CPP11_INITIALIZER_LISTS
269  void assign(std::initializer_list<std::size_t>, pointer,
270  const CoordinateOrder& = defaultOrder,
271  const CoordinateOrder& = defaultOrder,
272  const allocator_type& = allocator_type());
273  void assign(std::initializer_list<std::size_t>,
274  std::initializer_list<std::size_t>, pointer,
275  const CoordinateOrder&,
276  const allocator_type& = allocator_type());
277  #endif
278 
279  // query
280  const std::size_t dimension() const;
281  const std::size_t size() const;
282  const std::size_t shape(const std::size_t) const;
283  const std::size_t* shapeBegin() const;
284  const std::size_t* shapeEnd() const;
285  const std::size_t strides(const std::size_t) const;
286  const std::size_t* stridesBegin() const;
287  const std::size_t* stridesEnd() const;
288  const CoordinateOrder& coordinateOrder() const;
289  const bool isSimple() const;
290  template<class TLocal, bool isConstLocal, class ALocal>
291  bool overlaps(const View<TLocal, isConstLocal, ALocal>&) const;
292 
293  // element access
294  template<class U> reference operator()(U);
295  template<class U> reference operator()(U) const;
296  #ifndef HAVE_CPP11_VARIADIC_TEMPLATES
297  reference operator()(const std::size_t, const std::size_t);
298  reference operator()(const std::size_t, const std::size_t) const;
299  reference operator()(const std::size_t, const std::size_t, const std::size_t);
300  reference operator()(const std::size_t, const std::size_t, const std::size_t) const;
301  reference operator()(const std::size_t, const std::size_t, const std::size_t,
302  const std::size_t);
303  reference operator()(const std::size_t, const std::size_t, const std::size_t,
304  const std::size_t) const;
305  reference operator()(const std::size_t, const std::size_t, const std::size_t,
306  const std::size_t, const std::size_t);
307  reference operator()(const std::size_t, const std::size_t, const std::size_t,
308  const std::size_t, const std::size_t) const;
309  reference operator()(const std::size_t, const std::size_t, const std::size_t,
310  const std::size_t, const std::size_t, const std::size_t, const std::size_t,
311  const std::size_t, const std::size_t, const std::size_t);
312  reference operator()(const std::size_t, const std::size_t, const std::size_t,
313  const std::size_t, const std::size_t, const std::size_t, const std::size_t,
314  const std::size_t, const std::size_t, const std::size_t) const;
315  #else
316  reference operator()(const std::size_t);
317  reference operator()(const std::size_t) const;
318  template<typename... Args>
319  reference operator()(const std::size_t, const Args...);
320  template<typename... Args>
321  reference operator()(const std::size_t, const Args...) const;
322  private:
323  std::size_t elementAccessHelper(const std::size_t, const std::size_t);
324  std::size_t elementAccessHelper(const std::size_t, const std::size_t) const;
325  template<typename... Args>
326  std::size_t elementAccessHelper(const std::size_t, const std::size_t,
327  const Args...);
328  template<typename... Args>
329  std::size_t elementAccessHelper(const std::size_t, const std::size_t,
330  const Args...) const;
331  public:
332  #endif
333 
334  // sub-views
335  template<class BaseIterator, class ShapeIterator>
336  void view(BaseIterator, ShapeIterator, View<T, isConst, A>&) const;
337  template<class BaseIterator, class ShapeIterator>
338  void view(BaseIterator, ShapeIterator, const CoordinateOrder&,
339  View<T, isConst, A>&) const;
340  template<class BaseIterator, class ShapeIterator>
341  View<T, isConst, A> view(BaseIterator, ShapeIterator) const;
342  template<class BaseIterator, class ShapeIterator>
343  View<T, isConst, A> view(BaseIterator, ShapeIterator,
344  const CoordinateOrder&) const;
345  template<class BaseIterator, class ShapeIterator>
346  void constView(BaseIterator, ShapeIterator, View<T, true, A>&) const;
347  template<class BaseIterator, class ShapeIterator>
348  void constView(BaseIterator, ShapeIterator, const CoordinateOrder&,
349  View<T, true, A>&) const;
350  template<class BaseIterator, class ShapeIterator>
351  View<T, true, A> constView(BaseIterator, ShapeIterator) const;
352  template<class BaseIterator, class ShapeIterator>
353  View<T, true, A> constView(BaseIterator, ShapeIterator,
354  const CoordinateOrder&) const;
355  #ifdef HAVE_CPP11_INITIALIZER_LISTS
356  void view(std::initializer_list<std::size_t>,
357  std::initializer_list<std::size_t>, View<T, isConst, A>&) const;
358  void view(std::initializer_list<std::size_t>,
359  std::initializer_list<std::size_t>, const CoordinateOrder&,
360  View<T, isConst, A>&) const;
361  void constView(std::initializer_list<std::size_t>,
362  std::initializer_list<std::size_t>, View<T, true, A>&) const;
363  void constView(std::initializer_list<std::size_t>,
364  std::initializer_list<std::size_t>, const CoordinateOrder&,
365  View<T, true, A>&) const;
366  #endif
367 
368  // iterator access
369  iterator begin();
370  iterator end();
371  const_iterator begin() const;
372  const_iterator end() const;
373  reverse_iterator rbegin();
374  reverse_iterator rend();
375  const_reverse_iterator rbegin() const;
376  const_reverse_iterator rend() const;
377 
378  // coordinate transformation
379  template<class ShapeIterator>
380  void reshape(ShapeIterator, ShapeIterator);
381  template<class CoordinateIterator>
382  void permute(CoordinateIterator);
383  void transpose(const std::size_t, const std::size_t);
384  void transpose();
385  void shift(const int);
386  void squeeze();
387 
388  template<class ShapeIterator>
389  View<T, isConst, A> reshapedView(ShapeIterator, ShapeIterator) const;
390  template<class CoordinateIterator>
391  View<T, isConst, A> permutedView(CoordinateIterator) const;
392  View<T, isConst, A> transposedView(const std::size_t, const std::size_t) const;
393  View<T, isConst, A> transposedView() const;
394  View<T, isConst, A> shiftedView(const int) const;
395  View<T, isConst, A> boundView(const std::size_t, const std::size_t = 0) const;
396  View<T, isConst, A> squeezedView() const;
397 
398  #ifdef HAVE_CPP11_INITIALIZER_LISTS
399  void reshape(std::initializer_list<std::size_t>);
400  void permute(std::initializer_list<std::size_t>);
401 
402  View<T, isConst, A> reshapedView(std::initializer_list<std::size_t>) const;
403  View<T, isConst, A> permutedView(std::initializer_list<std::size_t>) const;
404  #endif
405 
406  // conversion between coordinates, index and offset
407  template<class CoordinateIterator>
408  void coordinatesToIndex(CoordinateIterator, std::size_t&) const;
409  template<class CoordinateIterator>
410  void coordinatesToOffset(CoordinateIterator, std::size_t&) const;
411  template<class CoordinateIterator>
412  void indexToCoordinates(std::size_t, CoordinateIterator) const;
413  void indexToOffset(std::size_t, std::size_t&) const;
414  #ifdef HAVE_CPP11_INITIALIZER_LISTS
415  void coordinatesToIndex(std::initializer_list<std::size_t>,
416  std::size_t&) const;
417  void coordinatesToOffset(std::initializer_list<std::size_t>,
418  std::size_t&) const;
419  #endif
420 
421  // output as string
422  std::string asString(const StringStyle& = MatrixStyle) const;
423 
424 private:
425  typedef typename marray_detail::Geometry<A> geometry_type;
426 
427  View(pointer, const geometry_type&);
428 
429  void updateSimplicity();
430  void testInvariant() const;
431 
432  // unsafe direct memory access
433  const T& operator[](const std::size_t) const;
434  T& operator[](const std::size_t);
435 
436  // data and memory address
437  pointer data_;
438  geometry_type geometry_;
439 
440 template<class TLocal, bool isConstLocal, class ALocal>
441  friend class View;
442 template<class TLocal, class ALocal>
443  friend class Marray;
444 // \cond suppress_doxygen
445 template<bool isConstTo, class TFrom, class TTo, class AFrom, class ATo>
446  friend struct marray_detail::AssignmentOperatorHelper;
447 friend struct marray_detail::AccessOperatorHelper<true>;
448 friend struct marray_detail::AccessOperatorHelper<false>;
449 
450 template<class E, class U>
451  friend class ViewExpression;
452 template<class E, class U, class UnaryFunctor>
453  friend class UnaryViewExpression;
454 template<class E1, class T1, class E2, class T2, class BinaryFunctor>
455  friend class BinaryViewExpression;
456 template<class E, class U, class S, class BinaryFunctor>
457  friend class BinaryViewExpressionScalarFirst;
458 template<class E, class U, class S, class BinaryFunctor>
459  friend class BinaryViewExpressionScalarSecond;
460 
461 template<class Functor, class T1, class Alocal, class E, class T2>
462  friend void marray_detail::operate(View<T1, false, Alocal>& v, const ViewExpression<E, T2>& expression, Functor f);
463 // \endcond end suppress_doxygen
464 };
465 
471 template<class T, bool isConst, class A>
472 class Iterator
473 {
474 public:
475  // STL random access iterator typedefs
476  typedef typename std::random_access_iterator_tag iterator_category;
477  typedef T value_type;
478  typedef std::ptrdiff_t difference_type;
479  typedef typename marray_detail::IfBool<isConst, const T*, T*>::type pointer;
480  typedef typename marray_detail::IfBool<isConst, const T&, T&>::type reference;
481 
482  // non-standard typedefs
483  typedef typename marray_detail::IfBool<isConst, const View<T, true, A>*,
484  View<T, false, A>*>::type view_pointer;
485  typedef typename marray_detail::IfBool<isConst, const View<T, true, A>&,
486  View<T, false, A>&>::type view_reference;
487 
488  // construction
489  Iterator();
490  Iterator(const View<T, false, A>&, const std::size_t = 0);
491  Iterator(View<T, false, A>&, const std::size_t = 0);
492  Iterator(const View<T, true, A>&, const std::size_t = 0);
493  Iterator(const Iterator<T, false, A>&);
494  // conversion from mutable to const
495 
496  // STL random access iterator operations
497  reference operator*() const;
498  pointer operator->() const;
499  reference operator[](const std::size_t) const;
500  Iterator<T, isConst, A>& operator+=(const difference_type&);
501  Iterator<T, isConst, A>& operator-=(const difference_type&);
502  Iterator<T, isConst, A>& operator++(); // prefix
503 
504  Iterator<T, isConst, A>& operator--(); // prefix
505  Iterator<T, isConst, A> operator++(int); // postfix
506  Iterator<T, isConst, A> operator--(int); // postfix
507  Iterator<T, isConst, A> operator+(const difference_type&) const;
508  Iterator<T, isConst, A> operator-(const difference_type&) const;
509  template<bool isConstLocal>
510  difference_type operator-(const Iterator<T, isConstLocal, A>&) const;
511  template<bool isConstLocal>
512  bool operator==(const Iterator<T, isConstLocal, A>&) const;
513  template<bool isConstLocal>
514  bool operator!=(const Iterator<T, isConstLocal, A>&) const;
515  template<bool isConstLocal>
516  bool operator<(const Iterator<T, isConstLocal, A>&) const;
517  template<bool isConstLocal>
518  bool operator>(const Iterator<T, isConstLocal, A>&) const;
519  template<bool isConstLocal>
520  bool operator<=(const Iterator<T, isConstLocal, A>&) const;
521  template<bool isConstLocal>
522  bool operator>=(const Iterator<T, isConstLocal, A>&) const;
523 
524  // operations beyond the STL standard
525  bool hasMore() const;
526  std::size_t index() const;
527  template<class CoordinateIterator>
528  void coordinate(CoordinateIterator) const;
529 
530 private:
531  void testInvariant() const;
532 
533  // attributes
534  view_pointer view_;
535  pointer pointer_;
536  std::size_t index_;
537  std::vector<std::size_t> coordinates_;
538 
539 friend class Marray<T, A>;
540 friend class Iterator<T, !isConst, A>; // for comparison operators
541 };
542 
544 template<class T, class A>
545 class Marray
546 : public View<T, false, A>
547 {
548 public:
549  typedef View<T, false, A> base;
550  typedef typename base::value_type value_type;
551  typedef typename base::pointer pointer;
552  typedef typename base::const_pointer const_pointer;
553  typedef typename base::reference reference;
554  typedef typename base::const_reference const_reference;
555  typedef typename base::iterator iterator;
556  typedef typename base::reverse_iterator reverse_iterator;
557  typedef typename base::const_iterator const_iterator;
558  typedef typename base::const_reverse_iterator const_reverse_iterator;
559  typedef typename A::template rebind<value_type>::other allocator_type;
560 
561  // constructors and destructor
562  Marray(const allocator_type& = allocator_type());
563  Marray(const T&, const CoordinateOrder& = defaultOrder,
564  const allocator_type& = allocator_type());
565  template<class ShapeIterator>
566  Marray(ShapeIterator, ShapeIterator, const T& = T(),
567  const CoordinateOrder& = defaultOrder,
568  const allocator_type& = allocator_type());
569  template<class ShapeIterator>
570  Marray(const InitializationSkipping&, ShapeIterator, ShapeIterator,
571  const CoordinateOrder& = defaultOrder,
572  const allocator_type& = allocator_type());
573  #ifdef HAVE_CPP11_INITIALIZER_LISTS
574  Marray(std::initializer_list<std::size_t>, const T& = T(),
575  const CoordinateOrder& = defaultOrder,
576  const allocator_type& = allocator_type());
577  #endif
578  Marray(const Marray<T, A>&);
579  template<class E, class Te>
580  Marray(const ViewExpression<E, Te>&,
581  const allocator_type& = allocator_type());
582  template<class TLocal, bool isConstLocal, class ALocal>
583  Marray(const View<TLocal, isConstLocal, ALocal>&);
584  ~Marray();
585 
586  // assignment
587  Marray<T, A>& operator=(const T&);
588  Marray<T, A>& operator=(const Marray<T, A>&); // over-write default
589  template<class TLocal, bool isConstLocal, class ALocal>
590  Marray<T, A>& operator=(const View<TLocal, isConstLocal, ALocal>&);
591  template<class E, class Te>
592  Marray<T, A>& operator=(const ViewExpression<E, Te>&);
593  void assign(const allocator_type& = allocator_type());
594 
595  // resize
596  template<class ShapeIterator>
597  void resize(ShapeIterator, ShapeIterator, const T& = T());
598  template<class ShapeIterator>
599  void resize(const InitializationSkipping&, ShapeIterator, ShapeIterator);
600  #ifdef HAVE_CPP11_INITIALIZER_LISTS
601  void resize(std::initializer_list<std::size_t>, const T& = T());
602  void resize(const InitializationSkipping&, std::initializer_list<std::size_t>);
603  #endif
604 
605 private:
606  typedef typename base::geometry_type geometry_type;
607 
608  void testInvariant() const;
609  template<bool SKIP_INITIALIZATION, class ShapeIterator>
610  void resizeHelper(ShapeIterator, ShapeIterator, const T& = T());
611 
612  allocator_type dataAllocator_;
613 };
614 
615 // implementation of View
616 
617 #ifdef HAVE_CPP11_INITIALIZER_LISTS
618 
619 
620 
621 
622 
623 
624 template<class T, bool isConst, class A>
625 inline void
626 View<T, isConst, A>::coordinatesToIndex
627 (
628  std::initializer_list<std::size_t> coordinate,
629  std::size_t& out
630 ) const
631 {
632  coordinatesToIndex(coordinate.begin(), out);
633 }
634 #endif
635 
642 template<class T, bool isConst, class A>
643 template<class CoordinateIterator>
644 inline void
645 View<T, isConst, A>::coordinatesToIndex
646 (
647  CoordinateIterator it,
648  std::size_t& out
649 ) const
650 {
651  testInvariant();
652  out = 0;
653  for(std::size_t j=0; j<this->dimension(); ++j, ++it) {
654  marray_detail::Assert(MARRAY_NO_ARG_TEST || static_cast<std::size_t>(*it) < shape(j));
655  out += static_cast<std::size_t>(*it) * geometry_.shapeStrides(j);
656  }
657 }
658 
659 #ifdef HAVE_CPP11_INITIALIZER_LISTS
660 
661 
662 
663 
664 
665 
666 template<class T, bool isConst, class A>
667 inline void
668 View<T, isConst, A>::coordinatesToOffset
669 (
670  std::initializer_list<std::size_t> coordinate,
671  std::size_t& out
672 ) const
673 {
674  coordinatesToOffset(coordinate.begin(), out);
675 }
676 #endif
677 
684 template<class T, bool isConst, class A>
685 template<class CoordinateIterator>
686 inline void
687 View<T, isConst, A>::coordinatesToOffset
688 (
689  CoordinateIterator it,
690  std::size_t& out
691 ) const
692 {
693  testInvariant();
694  out = 0;
695  for(std::size_t j=0; j<this->dimension(); ++j, ++it) {
696  marray_detail::Assert(MARRAY_NO_ARG_TEST || static_cast<std::size_t>(*it) < shape(j));
697  out += static_cast<std::size_t>(*it) * strides(j);
698  }
699 }
700 
708 template<class T, bool isConst, class A>
709 template<class CoordinateIterator>
710 inline void
711 View<T, isConst, A>::indexToCoordinates
712 (
713  std::size_t index, // copy to work on
714  CoordinateIterator outit
715 ) const
716 {
717  testInvariant();
718  marray_detail::Assert(MARRAY_NO_DEBUG || this->dimension() > 0);
719  marray_detail::Assert(MARRAY_NO_ARG_TEST || index < this->size());
720  if(coordinateOrder() == FirstMajorOrder) {
721  for(std::size_t j=0; j<this->dimension(); ++j, ++outit) {
722  *outit = std::size_t(index / geometry_.shapeStrides(j));
723  index = index % geometry_.shapeStrides(j);
724  }
725  }
726  else { // last major order
727  std::size_t j = this->dimension()-1;
728  outit += j;
729  for(;;) {
730  *outit = std::size_t(index / geometry_.shapeStrides(j));
731  index = index % geometry_.shapeStrides(j);
732  if(j == 0) {
733  break;
734  }
735  else {
736  --outit;
737  --j;
738  }
739  }
740  }
741 }
742 
749 template<class T, bool isConst, class A>
750 inline void
751 View<T, isConst, A>::indexToOffset
752 (
753  std::size_t index,
754  std::size_t& out
755 ) const
756 {
757  testInvariant();
758  marray_detail::Assert(MARRAY_NO_ARG_TEST || index < this->size());
759  if(isSimple()) {
760  out = index;
761  }
762  else {
763  out = 0;
764  if(coordinateOrder() == FirstMajorOrder) {
765  for(std::size_t j=0; j<this->dimension(); ++j) {
766  out += geometry_.strides(j) * (index / geometry_.shapeStrides(j));
767  index = index % geometry_.shapeStrides(j);
768  }
769  }
770  else { // last major order
771  if(this->dimension() == 0) {
772  marray_detail::Assert(MARRAY_NO_ARG_TEST || index == 0);
773  return;
774  }
775  else {
776  std::size_t j = this->dimension()-1;
777  for(;;) {
778  out += geometry_.strides(j) * (index / geometry_.shapeStrides(j));
779  index = index % geometry_.shapeStrides(j);
780  if(j == 0) {
781  break;
782  }
783  else {
784  --j;
785  }
786  }
787  }
788  }
789  }
790 }
791 
799 template<class T, bool isConst, class A>
800 inline
801 View<T, isConst, A>::View
802 (
803  const allocator_type& allocator
804 )
805 : data_(0),
806  geometry_(geometry_type(allocator))
807 {
808  testInvariant();
809 }
810 
811 // private constructor
812 template<class T, bool isConst, class A>
813 inline
814 View<T, isConst, A>::View
815 (
816  pointer data,
817  const geometry_type& geometry
818 )
819 : data_(data),
820  geometry_(geometry)
821 {
822  testInvariant();
823 }
824 
830 template<class T, bool isConst, class A>
831 inline
832 View<T, isConst, A>::View
833 (
834  pointer data,
835  const allocator_type& allocator
836 )
837 : data_(data),
838  geometry_(geometry_type(0, defaultOrder, 1, true, allocator))
839 {
840  testInvariant();
841 }
842 
847 template<class T, bool isConst, class A>
848 inline
849 View<T, isConst, A>::View
850 (
851  const View<T, false, A>& in
852 )
853 : data_(in.data_),
854  geometry_(in.geometry_)
855 {
856  testInvariant();
857 }
858 
871 
872 template<class T, bool isConst, class A>
873 template<class ShapeIterator>
874 inline
875 View<T, isConst, A>::View
876 (
877  ShapeIterator begin,
878  ShapeIterator end,
879  pointer data,
880  const CoordinateOrder& externalCoordinateOrder,
881  const CoordinateOrder& internalCoordinateOrder,
882  const allocator_type& allocator
883 )
884 : data_(data),
885  geometry_(begin, end, externalCoordinateOrder,
886  internalCoordinateOrder, allocator)
887 {
888  testInvariant();
889 }
890 
903 template<class T, bool isConst, class A>
904 template<class ShapeIterator, class StrideIterator>
905 inline
906 View<T, isConst, A>::View
907 (
908  ShapeIterator begin,
909  ShapeIterator end,
910  StrideIterator it,
911  pointer data,
912  const CoordinateOrder& internalCoordinateOrder,
913  const allocator_type& allocator
914 )
915 : data_(data),
916  geometry_(begin, end, it, internalCoordinateOrder, allocator)
917 {
918  testInvariant();
919 }
920 
921 #ifdef HAVE_CPP11_INITIALIZER_LISTS
922 
923 
924 
925 
926 
927 
928 
929 
930 
931 
932 template<class T, bool isConst, class A>
933 inline
934 View<T, isConst, A>::View
935 (
936  std::initializer_list<std::size_t> shape,
937  pointer data,
938  const CoordinateOrder& externalCoordinateOrder,
939  const CoordinateOrder& internalCoordinateOrder,
940  const allocator_type& allocator
941 )
942 : data_(data),
943  geometry_(shape.begin(), shape.end(), externalCoordinateOrder,
944  internalCoordinateOrder, allocator)
945 {
946  testInvariant();
947 }
948 
957 template<class T, bool isConst, class A>
958 inline
959 View<T, isConst, A>::View
960 (
961  std::initializer_list<std::size_t> shape,
962  std::initializer_list<std::size_t> strides,
963  pointer data,
964  const CoordinateOrder& internalCoordinateOrder,
965  const allocator_type& allocator
966 )
967 : data_(data),
968  geometry_(shape.begin(), shape.end(), strides.begin(),
969  internalCoordinateOrder, allocator)
970 {
971  testInvariant();
972 }
973 #endif
974 
982 template<class T, bool isConst, class A>
983 inline void
984 View<T, isConst, A>::assign
985 (
986  const allocator_type& allocator
987 )
988 {
989  geometry_ = geometry_type(allocator);
990  data_ = 0;
991  testInvariant();
992 }
993 
1006 template<class T, bool isConst, class A>
1007 template<class ShapeIterator>
1008 inline void
1009 View<T, isConst, A>::assign
1010 (
1011  ShapeIterator begin,
1012  ShapeIterator end,
1013  pointer data,
1014  const CoordinateOrder& externalCoordinateOrder,
1015  const CoordinateOrder& internalCoordinateOrder,
1016  const allocator_type& allocator
1017 )
1018 {
1019  // the invariant is not tested as a pre-condition of this
1020  // function to allow for unsafe manipulations prior to its
1021  // call
1022  geometry_ = typename marray_detail::Geometry<A>(begin, end,
1023  externalCoordinateOrder, internalCoordinateOrder, allocator);
1024  data_ = data;
1025  testInvariant();
1026 }
1027 
1040 template<class T, bool isConst, class A>
1041 template<class ShapeIterator, class StrideIterator>
1042 inline void
1043 View<T, isConst, A>::assign
1044 (
1045  ShapeIterator begin,
1046  ShapeIterator end,
1047  StrideIterator it,
1048  pointer data,
1049  const CoordinateOrder& internalCoordinateOrder,
1050  const allocator_type& allocator
1051 )
1052 {
1053  // the invariant is not tested as a pre-condition of this
1054  // function to allow for unsafe manipulations prior to its
1055  // call
1056  geometry_ = typename marray_detail::Geometry<A>(begin, end,
1057  it, internalCoordinateOrder, allocator);
1058  data_ = data;
1059  testInvariant();
1060 }
1061 
1062 #ifdef HAVE_CPP11_INITIALIZER_LISTS
1063 
1064 
1065 
1066 
1067 
1068 
1069 
1070 
1071 
1072 
1073 template<class T, bool isConst, class A>
1074 inline void
1075 View<T, isConst, A>::assign
1076 (
1077  std::initializer_list<std::size_t> shape,
1078  pointer data,
1079  const CoordinateOrder& externalCoordinateOrder,
1080  const CoordinateOrder& internalCoordinateOrder,
1081  const allocator_type& allocator
1082 )
1083 {
1084  assign(shape.begin(), shape.end(), data, externalCoordinateOrder,
1085  internalCoordinateOrder, allocator);
1086 }
1087 
1098 template<class T, bool isConst, class A>
1099 inline void
1100 View<T, isConst, A>::assign
1101 (
1102  std::initializer_list<std::size_t> shape,
1103  std::initializer_list<std::size_t> strides,
1104  pointer data,
1105  const CoordinateOrder& internalCoordinateOrder,
1106  const allocator_type& allocator
1107 )
1108 {
1109  assign(shape.begin(), shape.end(), strides.begin(), data,
1110  internalCoordinateOrder, allocator);
1111 }
1112 #endif
1113 
1121 template<class T, bool isConst, class A>
1122 template<class U>
1123 inline typename View<T, isConst, A>::reference
1124 View<T, isConst, A>::operator()
1125 (
1126  U u
1127 )
1128 {
1129  return marray_detail::AccessOperatorHelper<std::numeric_limits<U>::is_integer>::execute(*this, u);
1130 }
1131 
1139 template<class T, bool isConst, class A>
1140 template<class U>
1141 inline typename View<T, isConst, A>::reference
1142 View<T, isConst, A>::operator()
1143 (
1144  U u
1145 ) const
1146 {
1147  return marray_detail::AccessOperatorHelper<std::numeric_limits<U>::is_integer>::execute(*this, u);
1148 }
1149 
1150 #ifndef HAVE_CPP11_VARIADIC_TEMPLATES
1151 
1160 template<class T, bool isConst, class A>
1161 inline typename View<T, isConst, A>::reference
1162 View<T, isConst, A>::operator()
1163 (
1164  const std::size_t c0,
1165  const std::size_t c1
1166 )
1167 {
1168  testInvariant();
1169  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 2));
1170  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)));
1171  return data_[c0*strides(0) + c1*strides(1)];
1172 }
1173 
1182 template<class T, bool isConst, class A>
1183 inline typename View<T, isConst, A>::reference
1184 View<T, isConst, A>::operator()
1185 (
1186  const std::size_t c0,
1187  const std::size_t c1
1188 ) const
1189 {
1190  testInvariant();
1191  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 2));
1192  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)));
1193  return data_[c0*strides(0) + c1*strides(1)];
1194 }
1195 
1205 template<class T, bool isConst, class A>
1206 inline typename View<T, isConst, A>::reference
1207 View<T, isConst, A>::operator()
1208 (
1209  const std::size_t c0,
1210  const std::size_t c1,
1211  const std::size_t c2
1212 )
1213 {
1214  testInvariant();
1215  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 3));
1216  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
1217  && c2 < shape(2)));
1218  return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)];
1219 }
1220 
1230 template<class T, bool isConst, class A>
1231 inline typename View<T, isConst, A>::reference
1232 View<T, isConst, A>::operator()
1233 (
1234  const std::size_t c0,
1235  const std::size_t c1,
1236  const std::size_t c2
1237 ) const
1238 {
1239  testInvariant();
1240  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 3));
1241  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
1242  && c2 < shape(2)));
1243  return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)];
1244 }
1245 
1256 template<class T, bool isConst, class A>
1257 inline typename View<T, isConst, A>::reference
1258 View<T, isConst, A>::operator()
1259 (
1260  const std::size_t c0,
1261  const std::size_t c1,
1262  const std::size_t c2,
1263  const std::size_t c3
1264 )
1265 {
1266  testInvariant();
1267  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 4));
1268  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
1269  && c2 < shape(2) && c3 < shape(3)));
1270  return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
1271  + c3*strides(3)];
1272 }
1273 
1284 template<class T, bool isConst, class A>
1285 inline typename View<T, isConst, A>::reference
1286 View<T, isConst, A>::operator()
1287 (
1288  const std::size_t c0,
1289  const std::size_t c1,
1290  const std::size_t c2,
1291  const std::size_t c3
1292 ) const
1293 {
1294  testInvariant();
1295  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 4));
1296  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
1297  && c2 < shape(2) && c3 < shape(3)));
1298  return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
1299  + c3*strides(3)];
1300 }
1301 
1313 template<class T, bool isConst, class A>
1314 inline typename View<T, isConst, A>::reference
1315 View<T, isConst, A>::operator()
1316 (
1317  const std::size_t c0,
1318  const std::size_t c1,
1319  const std::size_t c2,
1320  const std::size_t c3,
1321  const std::size_t c4
1322 )
1323 {
1324  testInvariant();
1325  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 5));
1326  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
1327  && c2 < shape(2) && c3 < shape(3) && c4 < shape(4)));
1328  return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
1329  + c3*strides(3) + c4*strides(4)];
1330 }
1331 
1343 template<class T, bool isConst, class A>
1344 inline typename View<T, isConst, A>::reference
1345 View<T, isConst, A>::operator()
1346 (
1347  const std::size_t c0,
1348  const std::size_t c1,
1349  const std::size_t c2,
1350  const std::size_t c3,
1351  const std::size_t c4
1352 ) const
1353 {
1354  testInvariant();
1355  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 5));
1356  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
1357  && c2 < shape(2) && c3 < shape(3) && c4 < shape(4)));
1358  return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
1359  + c3*strides(3) + c4*strides(4)];
1360 }
1361 
1365 
1379 template<class T, bool isConst, class A>
1380 inline typename View<T, isConst, A>::reference
1381 View<T, isConst, A>::operator()
1382 (
1383  const std::size_t c0,
1384  const std::size_t c1,
1385  const std::size_t c2,
1386  const std::size_t c3,
1387  const std::size_t c4,
1388  const std::size_t c5,
1389  const std::size_t c6,
1390  const std::size_t c7,
1391  const std::size_t c8,
1392  const std::size_t c9
1393 )
1394 {
1395  testInvariant();
1396  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 10));
1397  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
1398  && c2 < shape(2) && c3 < shape(3) && c4 < shape(4)
1399  && c5 < shape(5) && c6 < shape(6) && c7 < shape(7)
1400  && c8 < shape(8) && c9 < shape(9)));
1401  return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
1402  + c3*strides(3) + c4*strides(4) + c5*strides(5) + c6*strides(6)
1403  + c7*strides(7) + c8*strides(8) + c9*strides(9)];
1404 }
1405 
1422 template<class T, bool isConst, class A>
1423 inline typename View<T, isConst, A>::reference
1424 View<T, isConst, A>::operator()
1425 (
1426  const std::size_t c0,
1427  const std::size_t c1,
1428  const std::size_t c2,
1429  const std::size_t c3,
1430  const std::size_t c4,
1431  const std::size_t c5,
1432  const std::size_t c6,
1433  const std::size_t c7,
1434  const std::size_t c8,
1435  const std::size_t c9
1436 ) const
1437 {
1438  testInvariant();
1439  marray_detail::Assert(MARRAY_NO_DEBUG || (data_ != 0 && dimension() == 10));
1440  marray_detail::Assert(MARRAY_NO_ARG_TEST || (c0 < shape(0) && c1 < shape(1)
1441  && c2 < shape(2) && c3 < shape(3) && c4 < shape(4)
1442  && c5 < shape(5) && c6 < shape(6) && c7 < shape(7)
1443  && c8 < shape(8) && c9 < shape(9)));
1444  return data_[c0*strides(0) + c1*strides(1) + c2*strides(2)
1445  + c3*strides(3) + c4*strides(4) + c5*strides(5) + c6*strides(6)
1446  + c7*strides(7) + c8*strides(8) + c9*strides(9)];
1447 }
1448 
1449 #else
1450 
1451 template<class T, bool isConst, class A>
1452 inline std::size_t
1453 View<T, isConst, A>::elementAccessHelper
1454 (
1455  const std::size_t Dim,
1456  const std::size_t value
1457 )
1458 {
1459  marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1) ) );
1460  return strides(Dim-1) * value;
1461 }
1462 
1463 template<class T, bool isConst, class A>
1464 inline std::size_t
1465 View<T, isConst, A>::elementAccessHelper
1466 (
1467  const std::size_t Dim,
1468  const std::size_t value
1469 ) const
1470 {
1471  marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1) ) );
1472  return strides(Dim-1) * value;
1473 }
1474 
1475 template<class T, bool isConst, class A>
1476 template<typename... Args>
1477 inline std::size_t
1478 View<T, isConst, A>::elementAccessHelper
1479 (
1480  const std::size_t Dim,
1481  const std::size_t value,
1482  const Args... args
1483 )
1484 {
1485  marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1-sizeof...(args)) ) );
1486  return value * strides(Dim-1-sizeof...(args)) + elementAccessHelper(Dim, args...);
1487 }
1488 
1489 template<class T, bool isConst, class A>
1490 template<typename... Args>
1491 inline std::size_t
1492 View<T, isConst, A>::elementAccessHelper
1493 (
1494  const std::size_t Dim,
1495  const std::size_t value,
1496  const Args... args
1497 ) const
1498 {
1499  marray_detail::Assert(MARRAY_NO_ARG_TEST || (value < shape(Dim-1-sizeof...(args)) ) );
1500  return value * strides(Dim-1-sizeof...(args)) + elementAccessHelper(Dim, args...);
1501 }
1502 
1503 template<class T, bool isConst, class A>
1504 inline typename View<T, isConst, A>::reference
1505 View<T, isConst, A>::operator()
1506 (
1507  const std::size_t value
1508 )
1509 {
1510  testInvariant();
1511  if(dimension() == 0) {
1512  marray_detail::Assert(MARRAY_NO_ARG_TEST || value == 0);
1513  return data_[0];
1514  }
1515  else {
1516  std::size_t offset;
1517  indexToOffset(value, offset);
1518  return data_[offset];
1519  }
1520 }
1521 
1522 template<class T, bool isConst, class A>
1523 inline typename View<T, isConst, A>::reference
1524 View<T, isConst, A>::operator()
1525 (
1526  const std::size_t value
1527 ) const
1528 {
1529  testInvariant();
1530  if(dimension() == 0) {
1531  marray_detail::Assert(MARRAY_NO_ARG_TEST || value == 0);
1532  return data_[0];
1533  }
1534  else {
1535  std::size_t offset;
1536  indexToOffset(value, offset);
1537  return data_[offset];
1538  }
1539 }
1540 
1541 template<class T, bool isConst, class A>
1542 template<typename... Args>
1543 inline typename View<T, isConst, A>::reference
1544 View<T, isConst, A>::operator()
1545 (
1546  const std::size_t value,
1547  const Args... args
1548 )
1549 {
1550  testInvariant();
1551  marray_detail::Assert( MARRAY_NO_DEBUG || ( data_ != 0 && sizeof...(args)+1 == dimension() ) );
1552  return data_[strides(0)*value + elementAccessHelper(sizeof...(args)+1, args...) ];
1553 }
1554 
1555 template<class T, bool isConst, class A>
1556 template<typename... Args>
1557 inline typename View<T, isConst, A>::reference
1558 View<T, isConst, A>::operator()
1559 (
1560  const std::size_t value,
1561  const Args... args
1562 ) const
1563 {
1564  testInvariant();
1565  marray_detail::Assert( MARRAY_NO_DEBUG || ( data_ != 0 && sizeof...(args)+1 == dimension() ) );
1566  return data_[ strides(0) * static_cast<std::size_t>(value)
1567  + static_cast<std::size_t>(elementAccessHelper(sizeof...(args)+1, args...)) ];
1568 }
1569 
1570 #endif // #ifndef HAVE_CPP11_VARIADIC_TEMPLATES
1571 
1576 template<class T, bool isConst, class A>
1577 inline const std::size_t
1578 View<T, isConst, A>::size() const
1579 {
1580  return geometry_.size();
1581 }
1582 
1589 template<class T, bool isConst, class A>
1590 inline const std::size_t
1591 View<T, isConst, A>::dimension() const
1592 {
1593  marray_detail::Assert(MARRAY_NO_DEBUG || this->data_ != 0);
1594  return geometry_.dimension();
1595 }
1596 
1602 template<class T, bool isConst, class A>
1603 inline const std::size_t
1604 View<T, isConst, A>::shape
1605 (
1606  const std::size_t dimension
1607 ) const
1608 {
1609  testInvariant();
1610  marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
1611  marray_detail::Assert(MARRAY_NO_ARG_TEST || dimension < this->dimension());
1612  return geometry_.shape(dimension);
1613 }
1614 
1620 template<class T, bool isConst, class A>
1621 inline const std::size_t*
1622 View<T, isConst, A>::shapeBegin() const
1623 {
1624  testInvariant();
1625  marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
1626  return geometry_.shapeBegin();
1627 }
1628 
1634 template<class T, bool isConst, class A>
1635 inline const std::size_t*
1636 View<T, isConst, A>::shapeEnd() const
1637 {
1638  testInvariant();
1639  marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
1640  return geometry_.shapeEnd();
1641 }
1642 
1648 template<class T, bool isConst, class A>
1649 inline const std::size_t
1650 View<T, isConst, A>::strides
1651 (
1652  const std::size_t dimension
1653 ) const
1654 {
1655  testInvariant();
1656  marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
1657  marray_detail::Assert(MARRAY_NO_ARG_TEST || dimension < this->dimension());
1658  return geometry_.strides(dimension);
1659 }
1660 
1666 template<class T, bool isConst, class A>
1667 inline const std::size_t*
1668 View<T, isConst, A>::stridesBegin() const
1669 {
1670  testInvariant();
1671  marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
1672  return geometry_.stridesBegin();
1673 }
1674 
1680 template<class T, bool isConst, class A>
1681 inline const std::size_t*
1682 View<T, isConst, A>::stridesEnd() const
1683 {
1684  testInvariant();
1685  marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
1686  return geometry_.stridesEnd();
1687 }
1688 
1693 template<class T, bool isConst, class A>
1694 inline const CoordinateOrder&
1695 View<T, isConst, A>::coordinateOrder() const
1696 {
1697  testInvariant();
1698  return geometry_.coordinateOrder();
1699 }
1700 
1705 template<class T, bool isConst, class A>
1706 inline const bool
1707 View<T, isConst, A>::isSimple() const
1708 {
1709  testInvariant();
1710  return geometry_.isSimple();
1711 }
1712 
1752 template<class T, bool isConst, class A>
1753 inline View<T, isConst, A>&
1754 View<T, isConst, A>::operator=
1755 (
1756  const View<T, true, A>& in
1757 )
1758 {
1759  testInvariant();
1760  marray_detail::AssignmentOperatorHelper<isConst, T, T, A, A>::execute(in, *this);
1761  testInvariant();
1762  return *this;
1763 }
1764 
1767 template<class T, bool isConst, class A>
1768 inline View<T, isConst, A>&
1769 View<T, isConst, A>::operator=
1770 (
1771  const View<T, false, A>& in
1772 )
1773 {
1774  testInvariant();
1775  marray_detail::AssignmentOperatorHelper<isConst, T, T, A, A>::execute(in, *this);
1776  testInvariant();
1777  return *this;
1778 }
1779 
1782 template<class T, bool isConst, class A>
1783 template<class TLocal, bool isConstLocal, class ALocal>
1784 inline View<T, isConst, A>&
1785 View<T, isConst, A>::operator=
1786 (
1787  const View<TLocal, isConstLocal, ALocal>& in
1788 )
1789 {
1790  testInvariant();
1791  marray_detail::AssignmentOperatorHelper<isConst, TLocal, T, ALocal, A>::execute(in, *this);
1792  testInvariant();
1793  return *this;
1794 }
1795 
1802 template<class T, bool isConst, class A>
1803 inline View<T, isConst, A>&
1804 View<T, isConst, A>::operator=
1805 (
1806  const T& value
1807 )
1808 {
1809  marray_detail::Assert(MARRAY_NO_DEBUG || data_ != 0);
1810  if(isSimple()) {
1811  for(std::size_t j=0; j<geometry_.size(); ++j) {
1812  data_[j] = value;
1813  }
1814  }
1815  else if(dimension() == 1)
1816  marray_detail::OperateHelperBinaryScalar<1, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1817  else if(dimension() == 2)
1818  marray_detail::OperateHelperBinaryScalar<2, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1819  else if(dimension() == 3)
1820  marray_detail::OperateHelperBinaryScalar<3, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1821  else if(dimension() == 4)
1822  marray_detail::OperateHelperBinaryScalar<4, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1823  else if(dimension() == 5)
1824  marray_detail::OperateHelperBinaryScalar<5, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1825  else if(dimension() == 6)
1826  marray_detail::OperateHelperBinaryScalar<6, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1827  else if(dimension() == 7)
1828  marray_detail::OperateHelperBinaryScalar<7, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1829  else if(dimension() == 8)
1830  marray_detail::OperateHelperBinaryScalar<8, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1831  else if(dimension() == 9)
1832  marray_detail::OperateHelperBinaryScalar<9, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1833  else if(dimension() == 10)
1834  marray_detail::OperateHelperBinaryScalar<10, marray_detail::Assign<T, T>, T, T, A>::operate(*this, value, marray_detail::Assign<T, T>(), data_);
1835  else {
1836  for(iterator it = begin(); it.hasMore(); ++it) {
1837  *it = value;
1838  }
1839  }
1840  return *this;
1841 }
1842 
1843 template<class T, bool isConst, class A>
1844 template<class E, class Te>
1845 inline View<T, isConst, A>&
1846 View<T, isConst, A>::operator=
1847 (
1848  const ViewExpression<E, Te>& expression
1849 )
1850 {
1851  marray_detail::operate(*this, expression, marray_detail::Assign<T, Te>());
1852  return *this;
1853 }
1854 
1863 template<class T, bool isConst, class A>
1864 template<class BaseIterator, class ShapeIterator>
1865 inline void
1866 View<T, isConst, A>::view
1867 (
1868  BaseIterator bit,
1869  ShapeIterator sit,
1870  View<T, isConst, A>& out
1871 ) const
1872 {
1873  view(bit, sit, coordinateOrder(), out);
1874 }
1875 
1886 template<class T, bool isConst, class A>
1887 template<class BaseIterator, class ShapeIterator>
1888 inline void
1889 View<T, isConst, A>::view
1890 (
1891  BaseIterator bit,
1892  ShapeIterator sit,
1893  const CoordinateOrder& internalCoordinateOrder,
1894  View<T, isConst, A>& out
1895 ) const
1896 {
1897  testInvariant();
1898  std::size_t offset = 0;
1899  coordinatesToOffset(bit, offset);
1900  out.assign(sit, sit+dimension(), geometry_.stridesBegin(),
1901  data_+offset, internalCoordinateOrder);
1902 }
1903 
1912 template<class T, bool isConst, class A>
1913 template<class BaseIterator, class ShapeIterator>
1914 inline View<T, isConst, A>
1915 View<T, isConst, A>::view
1916 (
1917  BaseIterator bit,
1918  ShapeIterator sit
1919 ) const
1920 {
1921  View<T, isConst, A> v;
1922  this->view(bit, sit, v);
1923  return v;
1924 }
1925 
1936 template<class T, bool isConst, class A>
1937 template<class BaseIterator, class ShapeIterator>
1938 inline View<T, isConst, A>
1939 View<T, isConst, A>::view
1940 (
1941  BaseIterator bit,
1942  ShapeIterator sit,
1943  const CoordinateOrder& internalCoordinateOrder
1944 ) const
1945 {
1946  View<T, isConst, A> v;
1947  this->view(bit, sit, internalCoordinateOrder, v);
1948  return v;
1949 }
1950 
1960 template<class T, bool isConst, class A>
1961 template<class BaseIterator, class ShapeIterator>
1962 inline void
1963 View<T, isConst, A>::constView
1964 (
1965  BaseIterator bit,
1966  ShapeIterator sit,
1967  View<T, true, A>& out
1968 ) const
1969 {
1970  constView(bit, sit, coordinateOrder(), out);
1971 }
1972 
1983 template<class T, bool isConst, class A>
1984 template<class BaseIterator, class ShapeIterator>
1985 inline void
1986 View<T, isConst, A>::constView
1987 (
1988  BaseIterator bit,
1989  ShapeIterator sit,
1990  const CoordinateOrder& internalCoordinateOrder,
1991  View<T, true, A>& out
1992 ) const
1993 {
1994  testInvariant();
1995  std::size_t offset = 0;
1996  coordinatesToOffset(bit, offset);
1997  out.assign(sit, sit+dimension(),
1998  geometry_.stridesBegin(),
1999  static_cast<const T*>(data_) + offset,
2000  internalCoordinateOrder);
2001 }
2002 
2012 template<class T, bool isConst, class A>
2013 template<class BaseIterator, class ShapeIterator>
2014 inline View<T, true, A>
2015 View<T, isConst, A>::constView
2016 (
2017  BaseIterator bit,
2018  ShapeIterator sit
2019 ) const
2020 {
2021  View<T, true, A> v;
2022  this->constView(bit, sit, v);
2023  return v;
2024 }
2025 
2036 template<class T, bool isConst, class A>
2037 template<class BaseIterator, class ShapeIterator>
2038 inline View<T, true, A>
2039 View<T, isConst, A>::constView
2040 (
2041  BaseIterator bit,
2042  ShapeIterator sit,
2043  const CoordinateOrder& internalCoordinateOrder
2044 ) const
2045 {
2046  View<T, true, A> v;
2047  this->constView(bit, sit, internalCoordinateOrder, v);
2048  return v;
2049 }
2050 
2051 #ifdef HAVE_CPP11_INITIALIZER_LISTS
2052 
2053 
2054 
2055 
2056 
2057 
2058 
2059 
2060 
2061 template<class T, bool isConst, class A>
2062 inline void
2063 View<T, isConst, A>::view
2064 (
2065  std::initializer_list<std::size_t> b,
2066  std::initializer_list<std::size_t> s,
2067  const CoordinateOrder& internalCoordinateOrder,
2068  View<T, isConst, A>& out
2069 ) const
2070 {
2071  view(b.begin(), s.begin(), internalCoordinateOrder, out);
2072 }
2073 
2082 template<class T, bool isConst, class A>
2083 inline void
2084 View<T, isConst, A>::view
2085 (
2086  std::initializer_list<std::size_t> b,
2087  std::initializer_list<std::size_t> s,
2088  View<T, isConst, A>& out
2089 ) const
2090 {
2091  view(b.begin(), s.begin(), coordinateOrder(), out);
2092 }
2093 
2103 template<class T, bool isConst, class A>
2104 inline void
2105 View<T, isConst, A>::constView
2106 (
2107  std::initializer_list<std::size_t> b,
2108  std::initializer_list<std::size_t> s,
2109  const CoordinateOrder& internalCoordinateOrder,
2110  View<T, true, A>& out
2111 ) const
2112 {
2113  constView(b.begin(), s.begin(), internalCoordinateOrder, out);
2114 }
2115 
2125 template<class T, bool isConst, class A>
2126 inline void
2127 View<T, isConst, A>::constView
2128 (
2129  std::initializer_list<std::size_t> b,
2130  std::initializer_list<std::size_t> s,
2131  View<T, true, A>& out
2132 ) const
2133 {
2134  constView(b.begin(), s.begin(), coordinateOrder(), out);
2135 }
2136 #endif
2137 
2151 template<class T, bool isConst, class A>
2152 template<class ShapeIterator>
2153 inline void
2154 View<T, isConst, A>::reshape
2155 (
2156  ShapeIterator begin,
2157  ShapeIterator end
2158 )
2159 {
2160  testInvariant();
2161  marray_detail::Assert(MARRAY_NO_DEBUG || isSimple());
2162  if(!MARRAY_NO_ARG_TEST) {
2163  std::size_t size = std::accumulate(begin, end, static_cast<std::size_t>(1),
2164  std::multiplies<std::size_t>());
2165  marray_detail::Assert(size == this->size());
2166  }
2167  assign(begin, end, data_, coordinateOrder(), coordinateOrder());
2168  testInvariant();
2169 }
2170 
2184 template<class T, bool isConst, class A>
2185 template<class ShapeIterator>
2186 inline View<T, isConst, A>
2187 View<T, isConst, A>::reshapedView
2188 (
2189  ShapeIterator begin,
2190  ShapeIterator end
2191 ) const
2192 {
2193  View<T, isConst, A> out = *this;
2194  out.reshape(begin, end);
2195  return out;
2196 }
2197 
2198 #ifdef HAVE_CPP11_INITIALIZER_LISTS
2199 
2200 
2201 
2202 
2203 
2204 
2205 
2206 
2207 
2208 
2209 
2210 template<class T, bool isConst, class A>
2211 inline void
2212 View<T, isConst, A>::reshape
2213 (
2214  std::initializer_list<std::size_t> shape
2215 )
2216 {
2217  reshape(shape.begin(), shape.end());
2218 }
2219 
2231 template<class T, bool isConst, class A>
2232 inline View<T, isConst, A>
2233 View<T, isConst, A>::reshapedView
2234 (
2235  std::initializer_list<std::size_t> shape
2236 ) const
2237 {
2238  return reshapedView(shape.begin(), shape.end());
2239 }
2240 #endif
2241 
2252 template<class T, bool isConst, class A>
2253 View<T, isConst, A>
2254 View<T, isConst, A>::boundView
2255 (
2256  const std::size_t dimension,
2257  const std::size_t value
2258 ) const
2259 {
2260  testInvariant();
2261  marray_detail::Assert(MARRAY_NO_ARG_TEST || (dimension < this->dimension()
2262  && value < shape(dimension)));
2263  if(this->dimension() == 1) {
2264  View v(&((*this)(value)));
2265  v.geometry_.coordinateOrder() = coordinateOrder();
2266  return v;
2267  }
2268  else {
2269  View v;
2270  v.geometry_.resize(this->dimension()-1);
2271  v.geometry_.coordinateOrder() = coordinateOrder();
2272  v.geometry_.size() = size() / shape(dimension);
2273  for(std::size_t j=0, k=0; j<this->dimension(); ++j) {
2274  if(j != dimension) {
2275  v.geometry_.shape(k) = shape(j);
2276  v.geometry_.strides(k) = strides(j);
2277  ++k;
2278  }
2279  }
2280  marray_detail::stridesFromShape(v.geometry_.shapeBegin(), v.geometry_.shapeEnd(),
2281  v.geometry_.shapeStridesBegin(), v.geometry_.coordinateOrder());
2282  v.data_ = data_ + strides(dimension) * value;
2283  v.updateSimplicity();
2284  v.testInvariant();
2285  return v;
2286  }
2287 }
2288 
2293 template<class T, bool isConst, class A>
2294 void
2295 View<T, isConst, A>::squeeze()
2296 {
2297  testInvariant();
2298  if(dimension() != 0) {
2299  std::size_t newDimension = dimension();
2300  for(std::size_t j=0; j<dimension(); ++j) {
2301  if(shape(j) == 1) {
2302  --newDimension;
2303  }
2304  }
2305  if(newDimension != dimension()) {
2306  if(newDimension == 0) {
2307  geometry_.resize(0);
2308  geometry_.size() = 1;
2309  }
2310  else {
2311  for(std::size_t j=0, k=0; j<geometry_.dimension(); ++j) {
2312  if(geometry_.shape(j) != 1) {
2313  geometry_.shape(k) = geometry_.shape(j);
2314  geometry_.strides(k) = geometry_.strides(j);
2315  ++k;
2316  }
2317  }
2318  geometry_.resize(newDimension);
2319  marray_detail::stridesFromShape(geometry_.shapeBegin(), geometry_.shapeEnd(),
2320  geometry_.shapeStridesBegin(), geometry_.coordinateOrder());
2321  updateSimplicity();
2322  }
2323  }
2324  }
2325  testInvariant();
2326 }
2327 
2333 template<class T, bool isConst, class A>
2334 inline View<T, isConst, A>
2335 View<T, isConst, A>::squeezedView() const
2336 {
2337  View<T, isConst, A> v = *this;
2338  v.squeeze();
2339  return v;
2340 }
2341 
2342 #ifdef HAVE_CPP11_INITIALIZER_LISTS
2343 
2344 
2345 
2346 
2347 
2348 
2349 
2350 
2351 template<class T, bool isConst, class A>
2352 void
2353 View<T, isConst, A>::permute
2354 (
2355  std::initializer_list<std::size_t> permutation
2356 )
2357 {
2358  permute(permutation.begin());
2359 }
2360 #endif
2361 
2370 template<class T, bool isConst, class A>
2371 template<class CoordinateIterator>
2372 void
2373 View<T, isConst, A>::permute
2374 (
2375  CoordinateIterator begin
2376 )
2377 {
2378  testInvariant();
2379  if(!MARRAY_NO_ARG_TEST) {
2380  marray_detail::Assert(dimension() != 0);
2381  std::set<std::size_t> s1, s2;
2382  CoordinateIterator it = begin;
2383  for(std::size_t j=0; j<dimension(); ++j) {
2384  s1.insert(j);
2385  s2.insert(*it);
2386  ++it;
2387  }
2388  marray_detail::Assert(s1 == s2);
2389  }
2390  // update shape, shape strides, strides, and simplicity
2391  std::vector<std::size_t> newShape = std::vector<std::size_t>(dimension());
2392  std::vector<std::size_t> newStrides = std::vector<std::size_t>(dimension());
2393  for(std::size_t j=0; j<dimension(); ++j) {
2394  newShape[j] = geometry_.shape(static_cast<std::size_t>(*begin));
2395  newStrides[j] = geometry_.strides(static_cast<std::size_t>(*begin));
2396  ++begin;
2397  }
2398  for(std::size_t j=0; j<dimension(); ++j) {
2399  geometry_.shape(j) = newShape[j];
2400  geometry_.strides(j) = newStrides[j];
2401  }
2402  marray_detail::stridesFromShape(geometry_.shapeBegin(), geometry_.shapeEnd(),
2403  geometry_.shapeStridesBegin(), geometry_.coordinateOrder());
2404  updateSimplicity();
2405  testInvariant();
2406 }
2407 
2417 template<class T, bool isConst, class A>
2418 template<class CoordinateIterator>
2419 inline View<T, isConst, A>
2420 View<T, isConst, A>::permutedView
2421 (
2422  CoordinateIterator begin
2423 ) const
2424 {
2425  View<T, isConst, A> out = *this;
2426  out.permute(begin);
2427  return out;
2428 }
2429 
2436 template<class T, bool isConst, class A>
2437 void
2438 View<T, isConst, A>::transpose
2439 (
2440  const std::size_t c1,
2441  const std::size_t c2
2442 )
2443 {
2444  testInvariant();
2445  marray_detail::Assert(MARRAY_NO_ARG_TEST ||
2446  (dimension() != 0 && c1 < dimension() && c2 < dimension()));
2447 
2448  std::size_t j1 = c1;
2449  std::size_t j2 = c2;
2450  std::size_t c;
2451  std::size_t d;
2452 
2453  // transpose shape
2454  c = geometry_.shape(j2);
2455  geometry_.shape(j2) = geometry_.shape(j1);
2456  geometry_.shape(j1) = c;
2457 
2458  // transpose strides
2459  d = geometry_.strides(j2);
2460  geometry_.strides(j2) = geometry_.strides(j1);
2461  geometry_.strides(j1) = d;
2462 
2463  // update shape strides
2464  marray_detail::stridesFromShape(geometry_.shapeBegin(), geometry_.shapeEnd(),
2465  geometry_.shapeStridesBegin(), geometry_.coordinateOrder());
2466 
2467  updateSimplicity();
2468  testInvariant();
2469 }
2470 
2476 template<class T, bool isConst, class A>
2477 void
2478 View<T, isConst, A>::transpose()
2479 {
2480  testInvariant();
2481  for(std::size_t j=0; j<dimension()/2; ++j) {
2482  std::size_t k = dimension()-j-1;
2483 
2484  // transpose shape
2485  std::size_t tmp = geometry_.shape(j);
2486  geometry_.shape(j) = geometry_.shape(k);
2487  geometry_.shape(k) = tmp;
2488 
2489  // transpose strides
2490  tmp = geometry_.strides(j);
2491  geometry_.strides(j) = geometry_.strides(k);
2492  geometry_.strides(k) = tmp;
2493  }
2494  marray_detail::stridesFromShape(geometry_.shapeBegin(), geometry_.shapeEnd(),
2495  geometry_.shapeStridesBegin(), geometry_.coordinateOrder());
2496  updateSimplicity();
2497  testInvariant();
2498 }
2499 
2508 template<class T, bool isConst, class A>
2509 inline View<T, isConst, A>
2510 View<T, isConst, A>::transposedView
2511 (
2512  const std::size_t c1,
2513  const std::size_t c2
2514 ) const
2515 {
2516  View<T, isConst, A> out = *this;
2517  out.transpose(c1, c2);
2518  return out;
2519 }
2520 
2527 template<class T, bool isConst, class A>
2528 inline View<T, isConst, A>
2529 View<T, isConst, A>::transposedView() const
2530 {
2531  View<T, isConst, A> out = *this;
2532  out.transpose();
2533  return out;
2534 }
2535 
2542 template<class T, bool isConst, class A>
2543 inline void
2544 View<T, isConst, A>::shift
2545 (
2546  const int n
2547 )
2548 {
2549  testInvariant();
2550  marray_detail::Assert(MARRAY_NO_DEBUG || dimension() != 0);
2551  if(n <= -static_cast<int>(dimension()) || n >= static_cast<int>(dimension())) {
2552  shift(n % static_cast<int>(dimension()));
2553  }
2554  else {
2555  if(n > 0) {
2556  shift(n - static_cast<int>(dimension()));
2557  }
2558  else {
2559  std::vector<std::size_t> p(dimension());
2560  for(std::size_t j=0; j<dimension(); ++j) {
2561  p[j] = static_cast<std::size_t>((static_cast<int>(j) - n)) % dimension();
2562  }
2563  permute(p.begin());
2564  }
2565  }
2566  testInvariant();
2567 }
2568 
2574 template<class T, bool isConst, class A>
2575 inline View<T, isConst, A>
2576 View<T, isConst, A>::shiftedView
2577 (
2578  const int n
2579 ) const
2580 {
2581  View<T, isConst, A> out = *this;
2582  out.shift(n);
2583  return out;
2584 }
2585 
2591 
2592 template<class T, bool isConst, class A>
2593 inline typename View<T, isConst, A>::iterator
2594 View<T, isConst, A>::begin()
2595 {
2596  testInvariant();
2597  return Iterator<T, isConst, A>(*this, 0);
2598 }
2599 
2605 template<class T, bool isConst, class A>
2606 inline typename View<T, isConst, A>::iterator
2607 View<T, isConst, A>::end()
2608 {
2609  testInvariant();
2610  return Iterator<T, isConst, A>(*this, geometry_.size());
2611 }
2612 
2618 template<class T, bool isConst, class A>
2619 inline typename View<T, isConst, A>::const_iterator
2620 View<T, isConst, A>::begin() const
2621 {
2622  testInvariant();
2623  return Iterator<T, true>(*this, 0);
2624 }
2625 
2631 template<class T, bool isConst, class A>
2632 inline typename View<T, isConst, A>::const_iterator
2633 
2634 View<T, isConst, A>::end() const
2635 {
2636  testInvariant();
2637  return Iterator<T, true>(*this, geometry_.size());
2638 }
2639 
2645 template<class T, bool isConst, class A>
2646 inline typename View<T, isConst, A>::reverse_iterator
2647 View<T, isConst, A>::rbegin()
2648 {
2649  return reverse_iterator(end());
2650 }
2651 
2657 template<class T, bool isConst, class A>
2658 inline typename View<T, isConst, A>::reverse_iterator
2659 View<T, isConst, A>::rend()
2660 {
2661  return reverse_iterator(begin());
2662 }
2663 
2669 template<class T, bool isConst, class A>
2670 inline typename View<T, isConst, A>::const_reverse_iterator
2671 View<T, isConst, A>::rbegin() const
2672 {
2673  return const_reverse_iterator(end());
2674 }
2675 
2681 template<class T, bool isConst, class A>
2682 inline typename View<T, isConst, A>::const_reverse_iterator
2683 View<T, isConst, A>::rend() const
2684 {
2685  return const_reverse_iterator(begin());
2686 }
2687 
2694 template<class T, bool isConst, class A>
2695 inline void
2696 View<T, isConst, A>::updateSimplicity()
2697 {
2698  // no invariant test here because this function
2699  // is called during unsafe updates of a view
2700  geometry_.updateSimplicity();
2701 }
2702 
2711 template<class T, bool isConst, class A>
2712 inline const T&
2713 View<T, isConst, A>::operator[]
2714 (
2715  const std::size_t offset
2716 )
2717 const
2718 {
2719  return data_[offset];
2720 }
2721 
2730 template<class T, bool isConst, class A>
2731 inline T&
2732 View<T, isConst, A>::operator[]
2733 (
2734  const std::size_t offset
2735 )
2736 {
2737  return data_[offset];
2738 }
2739 
2745 template<class T, bool isConst, class A>
2746 inline void
2747 View<T, isConst, A>::testInvariant() const
2748 {
2749  if(!MARRAY_NO_DEBUG) {
2750  if(geometry_.dimension() == 0) {
2751  marray_detail::Assert(geometry_.isSimple() == true);
2752  if(data_ != 0) { // scalar
2753  marray_detail::Assert(geometry_.size() == 1);
2754  }
2755  }
2756  else {
2757  marray_detail::Assert(data_ != 0);
2758 
2759  // test size_ to be consistent with shape_
2760  std::size_t testSize = 1;
2761  for(std::size_t j=0; j<geometry_.dimension(); ++j) {
2762  testSize *= geometry_.shape(j);
2763  }
2764  marray_detail::Assert(geometry_.size() == testSize);
2765 
2766  // test shapeStrides_ to be consistent with shape_
2767  if(geometry_.coordinateOrder() == FirstMajorOrder) {
2768  std::size_t tmp = 1;
2769  for(std::size_t j=0; j<geometry_.dimension(); ++j) {
2770  marray_detail::Assert(geometry_.shapeStrides(geometry_.dimension()-j-1) == tmp);
2771  tmp *= geometry_.shape(geometry_.dimension()-j-1);
2772  }
2773  }
2774  else {
2775  std::size_t tmp = 1;
2776  for(std::size_t j=0; j<geometry_.dimension(); ++j) {
2777  marray_detail::Assert(geometry_.shapeStrides(j) == tmp);
2778  tmp *= geometry_.shape(j);
2779  }
2780  }
2781 
2782  // test the simplicity condition
2783  if(geometry_.isSimple()) {
2784  for(std::size_t j=0; j<geometry_.dimension(); ++j) {
2785  marray_detail::Assert(geometry_.strides(j) == geometry_.shapeStrides(j));
2786  }
2787  }
2788  }
2789  }
2790 }
2791 
2807 template<class T, bool isConst, class A>
2808 template<class TLocal, bool isConstLocal, class ALocal>
2809 inline bool View<T, isConst, A>::overlaps
2810 (
2811  const View<TLocal, isConstLocal, ALocal>& v
2812 ) const
2813 {
2814  testInvariant();
2815  if(!MARRAY_NO_ARG_TEST) {
2816  v.testInvariant();
2817  }
2818  if(data_ == 0 || v.data_ == 0) {
2819  return false;
2820  }
2821  else {
2822  const void* dataPointer_ = data_;
2823  const void* vDataPointer_ = v.data_;
2824  const void* maxPointer = & (*this)(this->size()-1);
2825  const void* maxPointerV = & v(v.size()-1);
2826  if( (dataPointer_ <= vDataPointer_ && vDataPointer_ <= maxPointer)
2827  || (vDataPointer_ <= dataPointer_ && dataPointer_ <= maxPointerV) )
2828  {
2829  return true;
2830  }
2831  }
2832  return false;
2833 }
2834 
2837 template<class T, bool isConst, class A>
2838 std::string
2839 View<T, isConst, A>::asString
2840 (
2841  const StringStyle& style
2842 ) const
2843 {
2844  testInvariant();
2845  std::ostringstream out(std::ostringstream::out);
2846  if(style == MatrixStyle) {
2847  if(dimension() == 0) {
2848  // scalar
2849  out << "A = " << (*this)(0) << std::endl;
2850  }
2851  else if(dimension() == 1) {
2852  // vector
2853  out << "A = (";
2854  for(std::size_t j=0; j<this->size(); ++j) {
2855  out << (*this)(j) << ", ";
2856  }
2857  out << "\b\b)" << std::endl;
2858  }
2859  else if(dimension() == 2) {
2860  // matrix
2861  if(coordinateOrder() == FirstMajorOrder) {
2862  out << "A(r,c) =" << std::endl;
2863  for(std::size_t y=0; y<this->shape(0); ++y) {
2864  for(std::size_t x=0; x<this->shape(1); ++x) {
2865  out << (*this)(y, x) << ' ';
2866  }
2867  out << std::endl;
2868  }
2869  }
2870  else {
2871  out << "A(c,r) =" << std::endl;
2872  for(std::size_t y=0; y<this->shape(1); ++y) {
2873  for(std::size_t x=0; x<this->shape(0); ++x) {
2874  out << (*this)(x, y) << ' ';
2875  }
2876  out << std::endl;
2877  }
2878  }
2879  }
2880  else {
2881  // higher dimensional
2882  std::vector<std::size_t> c1(dimension());
2883  std::vector<std::size_t> c2(dimension());
2884  unsigned short q = 2;
2885  if(coordinateOrder() == FirstMajorOrder) {
2886  q = static_cast<unsigned char>(dimension() - 3);
2887  }
2888  for(const_iterator it = this->begin(); it.hasMore(); ++it) {
2889  it.coordinate(c2.begin());
2890  if(it.index() == 0 || c2[q] != c1[q]) {
2891  if(it.index() != 0) {
2892  out << std::endl << std::endl;
2893  }
2894  if(coordinateOrder() == FirstMajorOrder) {
2895  out << "A(";
2896  for(std::size_t j=0; j<dimension()-2; ++j) {
2897  out << c2[j] << ",";
2898  }
2899  }
2900  else {
2901  out << "A(c,r,";
2902  for(std::size_t j=2; j<dimension(); ++j) {
2903  out << c2[j] << ",";
2904  }
2905  }
2906  out << '\b';
2907  if(coordinateOrder() == FirstMajorOrder) {
2908  out << ",r,c";
2909  }
2910  out << ") =" << std::endl;
2911  }
2912  else if(c2[1] != c1[1]) {
2913  out << std::endl;
2914  }
2915  out << *it << " ";
2916  c1 = c2;
2917  }
2918  out << std::endl;
2919  }
2920  out << std::endl;
2921  }
2922  else if(style == TableStyle) {
2923  if(dimension() == 0) {
2924  // scalar
2925  out << "A = " << (*this)(0) << std::endl;
2926  }
2927  else {
2928  // non-scalar
2929  std::vector<std::size_t> c(dimension());
2930  for(const_iterator it = this->begin(); it.hasMore(); ++it) {
2931  out << "A(";
2932  it.coordinate(c.begin());
2933  for(std::size_t j=0; j<c.size(); ++j) {
2934  out << c[j] << ',';
2935  }
2936  out << "\b) = " << *it << std::endl;
2937  }
2938  }
2939  out << std::endl;
2940  }
2941  return out.str();
2942 }
2943 
2944 // implementation of arithmetic operators of View
2945 
2946 template<class T1, class T2, bool isConst, class A>
2947 inline View<T1, false, A>&
2948 operator+=
2949 (
2950  View<T1, false, A>& v,
2951  const View<T2, isConst, A>& w
2952 )
2953 {
2954  marray_detail::operate(v, w, marray_detail::PlusEqual<T1, T2>());
2955  return v;
2956 }
2957 
2958 // prefix
2959 template<class T, class A>
2960 inline View<T, false, A>&
2961 operator++
2962 (
2963  View<T, false, A>& v
2964 )
2965 {
2966  marray_detail::operate(v, marray_detail::PrefixIncrement<T>());
2967  return v;
2968 }
2969 
2970 // postfix
2971 template<class T, class A>
2972 inline Marray<T, A>
2973 operator++
2974 (
2975  Marray<T, A>& in,
2976  int dummy
2977 )
2978 {
2979  Marray<T, A> out = in;
2980  marray_detail::operate(in, marray_detail::PostfixIncrement<T>());
2981  return out;
2982 }
2983 
2984 template<class T1, class T2, bool isConst, class A>
2985 inline View<T1, false, A>&
2986 operator-=
2987 (
2988  View<T1, false, A>& v,
2989  const View<T2, isConst, A>& w
2990 )
2991 {
2992  marray_detail::operate(v, w, marray_detail::MinusEqual<T1, T2>());
2993  return v;
2994 }
2995 
2996 // prefix
2997 template<class T, class A>
2998 inline View<T, false, A>&
2999 operator--
3000 (
3001  View<T, false, A>& v
3002 )
3003 {
3004  marray_detail::operate(v, marray_detail::PrefixDecrement<T>());
3005  return v;
3006 }
3007 
3008 // postfix
3009 template<class T, class A>
3010 inline Marray<T, A>
3011 operator--
3012 (
3013  Marray<T, A>& in,
3014  int dummy
3015 )
3016 {
3017  Marray<T, A> out = in;
3018  marray_detail::operate(in, marray_detail::PostfixDecrement<T>());
3019  return out;
3020 }
3021 
3022 template<class T1, class T2, bool isConst, class A>
3023 inline View<T1, false, A>&
3024 operator*=
3025 (
3026  View<T1, false, A>& v,
3027  const View<T2, isConst, A>& w
3028 )
3029 {
3030  marray_detail::operate(v, w, marray_detail::TimesEqual<T1, T2>());
3031  return v;
3032 }
3033 
3034 template<class T1, class T2, bool isConst, class A>
3035 inline View<T1, false, A>&
3036 operator/=
3037 (
3038  View<T1, false, A>& v,
3039  const View<T2, isConst, A>& w
3040 )
3041 {
3042  marray_detail::operate(v, w, marray_detail::DividedByEqual<T1, T2>());
3043  return v;
3044 }
3045 
3046 template<class E1, class T1, class E2, class T2>
3047 inline const BinaryViewExpression<E1, T1, E2, T2, marray_detail::Plus<T1, T2, typename marray_detail::PromoteType<T1, T2>::type> >
3048 operator+
3049 (
3050  const ViewExpression<E1, T1>& expression1,
3051  const ViewExpression<E2, T2>& expression2
3052 )
3053 {
3054  typedef typename marray_detail::PromoteType<T1, T2>::type promoted_type;
3055  typedef marray_detail::Plus<T1, T2, promoted_type> Functor;
3056  typedef BinaryViewExpression<E1, T1, E2, T2, Functor> return_type;
3057  return return_type(expression1, expression2);
3058 }
3059 
3060 template<class E, class T>
3061 inline const ViewExpression<E,T>&
3062 operator+
3063 (
3064  const ViewExpression<E,T>& expression
3065 ) // unary
3066 {
3067  return expression;
3068 }
3069 
3070 #define MARRAY_UNARY_OPERATOR(datatype, operation, functorname) \
3071 template<class T, class A> \
3072 inline View<T, false, A>& \
3073 operator operation \
3074 ( \
3075  View<T, false, A>& v, \
3076  const datatype& x \
3077 ) \
3078 { \
3079  marray_detail::operate(v, static_cast<T>(x), marray_detail:: functorname <T, T>()); \
3080  return v; \
3081 } \
3082 
3083 #define MARRAY_UNARY_OPERATOR_ALL_TYPES(op, functorname) \
3084  MARRAY_UNARY_OPERATOR(char, op, functorname) \
3085  MARRAY_UNARY_OPERATOR(unsigned char, op, functorname) \
3086  MARRAY_UNARY_OPERATOR(short, op, functorname) \
3087  MARRAY_UNARY_OPERATOR(unsigned short, op, functorname) \
3088  MARRAY_UNARY_OPERATOR(int, op, functorname) \
3089  MARRAY_UNARY_OPERATOR(unsigned int, op, functorname) \
3090  MARRAY_UNARY_OPERATOR(long, op, functorname) \
3091  MARRAY_UNARY_OPERATOR(unsigned long, op, functorname) \
3092  MARRAY_UNARY_OPERATOR(float, op, functorname) \
3093  MARRAY_UNARY_OPERATOR(double, op, functorname) \
3094  MARRAY_UNARY_OPERATOR(long double, op, functorname) \
3095 
3096 MARRAY_UNARY_OPERATOR_ALL_TYPES(+=, PlusEqual)
3097 MARRAY_UNARY_OPERATOR_ALL_TYPES(-=, MinusEqual)
3098 MARRAY_UNARY_OPERATOR_ALL_TYPES(*=, TimesEqual)
3099 MARRAY_UNARY_OPERATOR_ALL_TYPES(/=, DividedByEqual)
3100 
3101 template<class E1, class T1, class E2, class T2>
3102 inline const BinaryViewExpression<E1, T1, E2, T2,
3103  marray_detail::Minus<T1, T2, typename marray_detail::PromoteType<T1, T2>::type> >
3104 operator-
3105 (
3106  const ViewExpression<E1, T1>& expression1,
3107  const ViewExpression<E2, T2>& expression2
3108 )
3109 {
3110  return BinaryViewExpression<E1, T1, E2, T2,
3111  marray_detail::Minus<T1, T2,
3112  typename marray_detail::PromoteType<T1, T2>::type> >(
3113  expression1, expression2);
3114 }
3115 
3116 template<class E, class T>
3117 inline const UnaryViewExpression<E,T,marray_detail::Negate<T> >
3118 operator-
3119 (
3120  const ViewExpression<E,T>& expression
3121 ) // unary
3122 {
3123  return UnaryViewExpression<E,T,marray_detail::Negate<T> >(
3124  expression);
3125 }
3126 
3127 template<class E1, class T1, class E2, class T2>
3128 inline const BinaryViewExpression<E1, T1, E2, T2,
3129  marray_detail::Times<T1, T2, typename marray_detail::PromoteType<T1, T2>::type> >
3130 operator*
3131 (
3132  const ViewExpression<E1, T1>& expression1,
3133  const ViewExpression<E2, T2>& expression2
3134 )
3135 {
3136  return BinaryViewExpression<E1, T1, E2, T2,
3137  marray_detail::Times<T1, T2,
3138  typename marray_detail::PromoteType<T1, T2>::type > >(
3139  expression1, expression2);
3140 }
3141 
3142 template<class E1, class T1, class E2, class T2>
3143 inline const BinaryViewExpression<E1, T1, E2, T2,
3144  marray_detail::DividedBy<T1, T2, typename marray_detail::PromoteType<T1, T2>::type> >
3145 operator/
3146 (
3147  const ViewExpression<E1, T1>& expression1,
3148  const ViewExpression<E2, T2>& expression2
3149 )
3150 {
3151  return BinaryViewExpression<E1, T1, E2, T2,
3152  marray_detail::DividedBy<T1, T2,
3153  typename marray_detail::PromoteType<T1, T2>::type > >(
3154  expression1, expression2);
3155 }
3156 
3157 #define MARRAY_BINARY_OPERATOR(datatype, operation, functorname) \
3158 template<class E, class T> \
3159 inline const BinaryViewExpressionScalarSecond< \
3160  E, T, datatype, marray_detail:: functorname < \
3161  T, datatype, typename marray_detail::PromoteType<T, datatype>::type \
3162  > \
3163 > \
3164 operator operation \
3165 ( \
3166  const ViewExpression<E, T>& expression, \
3167  const datatype& scalar \
3168 ) \
3169 { \
3170  typedef typename marray_detail::PromoteType<T, datatype>::type \
3171  promoted_type; \
3172  typedef marray_detail:: functorname <T, datatype, promoted_type> Functor; \
3173  typedef BinaryViewExpressionScalarSecond<E, T, datatype, Functor> \
3174  expression_type; \
3175  return expression_type(expression, scalar); \
3176 } \
3177 \
3178 template<class E, class T> \
3179 inline const BinaryViewExpressionScalarFirst \
3180 < \
3181  E, T, datatype, marray_detail:: functorname < \
3182  datatype, T, typename marray_detail::PromoteType<datatype, T>::type \
3183  > \
3184 > \
3185 operator operation \
3186 ( \
3187  const datatype& scalar, \
3188  const ViewExpression<E, T>& expression \
3189 ) \
3190 { \
3191  typedef typename marray_detail::PromoteType<T, datatype>::type \
3192  promoted_type; \
3193  typedef marray_detail:: functorname <datatype, T, promoted_type> Functor; \
3194  typedef BinaryViewExpressionScalarFirst<E, T, datatype, Functor> \
3195  expression_type; \
3196  return expression_type(expression, scalar); \
3197 }
3198 
3199 #define MARRAY_BINARY_OPERATOR_ALL_TYPES(op, functorname) \
3200  MARRAY_BINARY_OPERATOR(char, op, functorname) \
3201  MARRAY_BINARY_OPERATOR(unsigned char, op, functorname) \
3202  MARRAY_BINARY_OPERATOR(short, op, functorname) \
3203  MARRAY_BINARY_OPERATOR(unsigned short, op, functorname) \
3204  MARRAY_BINARY_OPERATOR(int, op, functorname) \
3205  MARRAY_BINARY_OPERATOR(unsigned int, op, functorname) \
3206  MARRAY_BINARY_OPERATOR(long, op, functorname) \
3207  MARRAY_BINARY_OPERATOR(unsigned long, op, functorname) \
3208  MARRAY_BINARY_OPERATOR(float, op, functorname) \
3209  MARRAY_BINARY_OPERATOR(double, op, functorname) \
3210  MARRAY_BINARY_OPERATOR(long double, op, functorname) \
3211 
3212 MARRAY_BINARY_OPERATOR_ALL_TYPES(+, Plus)
3213 MARRAY_BINARY_OPERATOR_ALL_TYPES(-, Minus)
3214 MARRAY_BINARY_OPERATOR_ALL_TYPES(*, Times)
3215 MARRAY_BINARY_OPERATOR_ALL_TYPES(/, DividedBy)
3216 
3217 // implementation of Marray
3218 
3227 template<class T, class A>
3228 inline void
3229 Marray<T, A>::assign
3230 (
3231  const allocator_type& allocator
3232 )
3233 {
3234  if(this->data_ != 0) {
3235  dataAllocator_.deallocate(this->data_, this->size());
3236  this->data_ = 0;
3237  }
3238  dataAllocator_ = allocator;
3239  base::assign();
3240 }
3241 
3246 template<class T, class A>
3247 inline
3248 Marray<T, A>::Marray
3249 (
3250  const allocator_type& allocator
3251 )
3252 : base(allocator),
3253  dataAllocator_(allocator)
3254 {
3255  testInvariant();
3256 }
3257 
3267 template<class T, class A>
3268 inline
3269 Marray<T, A>::Marray
3270 (
3271  const T& value,
3272  const CoordinateOrder& coordinateOrder,
3273  const allocator_type& allocator
3274 )
3275 : dataAllocator_(allocator)
3276 {
3277  this->data_ = dataAllocator_.allocate(1);
3278  this->data_[0] = value;
3279  this->geometry_ = geometry_type(0, coordinateOrder, 1, true, allocator);
3280  testInvariant();
3281 }
3282 
3287 template<class T, class A>
3288 inline
3289 Marray<T, A>::Marray
3290 (
3291  const Marray<T, A>& in
3292 )
3293 : dataAllocator_(in.dataAllocator_)
3294 {
3295  if(!MARRAY_NO_ARG_TEST) {
3296  in.testInvariant();
3297  }
3298  if(in.data_ == 0) {
3299  this->data_ = 0;
3300  }
3301  else {
3302  this->data_ = dataAllocator_.allocate(in.size());
3303  memcpy(this->data_, in.data_, (in.size())*sizeof(T));
3304  }
3305  this->geometry_ = in.geometry_;
3306  testInvariant();
3307 }
3308 
3313 template<class T, class A>
3314 template<class TLocal, bool isConstLocal, class ALocal>
3315 inline
3316 Marray<T, A>::Marray
3317 (
3318  const View<TLocal, isConstLocal, ALocal>& in
3319 )
3320 : dataAllocator_()
3321 {
3322  if(!MARRAY_NO_ARG_TEST) {
3323  in.testInvariant();
3324  }
3325 
3326  // adapt geometry
3327  this->geometry_ = in.geometry_;
3328  for(std::size_t j=0; j<in.dimension(); ++j) {
3329  this->geometry_.strides(j) = in.geometry_.shapeStrides(j); // !
3330  }
3331  this->geometry_.isSimple() = true;
3332 
3333  // copy data
3334  if(in.size() == 0) {
3335  this->data_ = 0;
3336  }
3337  else {
3338  this->data_ = dataAllocator_.allocate(in.size());
3339  }
3340  if(in.isSimple() && marray_detail::IsEqual<T, TLocal>::type) {
3341  memcpy(this->data_, in.data_, (in.size())*sizeof(T));
3342  }
3343  else {
3344  typename View<TLocal, isConstLocal, ALocal>::const_iterator it = in.begin();
3345  for(std::size_t j=0; j<this->size(); ++j, ++it) {
3346  this->data_[j] = static_cast<T>(*it);
3347  }
3348  }
3349 
3350  testInvariant();
3351 }
3352 
3358 template<class T, class A>
3359 template<class E, class Te>
3360 inline
3361 Marray<T, A>::Marray
3362 (
3363  const ViewExpression<E, Te>& expression,
3364  const allocator_type& allocator
3365 )
3366 : dataAllocator_(allocator)
3367 {
3368  this->data_ = dataAllocator_.allocate(expression.size());
3369  if(expression.dimension() == 0) {
3370  this->geometry_ = geometry_type(0,
3371  static_cast<const E&>(expression).coordinateOrder(),
3372  1, true, dataAllocator_);
3373  }
3374  else {
3375  this->geometry_ = geometry_type(
3376  static_cast<const E&>(expression).shapeBegin(),
3377  static_cast<const E&>(expression).shapeEnd(),
3378  static_cast<const E&>(expression).coordinateOrder(),
3379  static_cast<const E&>(expression).coordinateOrder(),
3380  dataAllocator_);
3381 
3382  }
3383  const E& e = static_cast<const E&>(expression);
3384  if(e.dimension() == 0) {
3385  marray_detail::Assert(MARRAY_NO_ARG_TEST || e.size() < 2);
3386  this->data_[0] = static_cast<T>(e(0));
3387  }
3388  else {
3389  marray_detail::Assert(MARRAY_NO_ARG_TEST || e.size() != 0);
3390  marray_detail::operate(*this, e, marray_detail::Assign<T, Te>());
3391  }
3392  testInvariant();
3393 }
3394 
3405 template<class T, class A>
3406 template<class ShapeIterator>
3407 inline
3408 Marray<T, A>::Marray
3409 (
3410  ShapeIterator begin,
3411  ShapeIterator end,
3412  const T& value,
3413  const CoordinateOrder& coordinateOrder,
3414  const allocator_type& allocator
3415 )
3416 : dataAllocator_(allocator)
3417 {
3418  std::size_t size = std::accumulate(begin, end, static_cast<std::size_t>(1),
3419  std::multiplies<std::size_t>());
3420  marray_detail::Assert(MARRAY_NO_ARG_TEST || size != 0);
3421  base::assign(begin, end, dataAllocator_.allocate(size), coordinateOrder,
3422  coordinateOrder, allocator);
3423  for(std::size_t j=0; j<size; ++j) {
3424  this->data_[j] = value;
3425  }
3426  testInvariant();
3427 }
3428 
3439 template<class T, class A>
3440 template<class ShapeIterator>
3441 inline
3442 Marray<T, A>::Marray
3443 (
3444  const InitializationSkipping& is,
3445  ShapeIterator begin,
3446  ShapeIterator end,
3447  const CoordinateOrder& coordinateOrder,
3448  const allocator_type& allocator
3449 )
3450 : dataAllocator_(allocator)
3451 {
3452  std::size_t size = std::accumulate(begin, end, static_cast<std::size_t>(1),
3453  std::multiplies<std::size_t>());
3454  marray_detail::Assert(MARRAY_NO_ARG_TEST || size != 0);
3455  base::assign(begin, end, dataAllocator_.allocate(size), coordinateOrder,
3456  coordinateOrder, allocator);
3457  testInvariant();
3458 }
3459 
3460 #ifdef HAVE_CPP11_INITIALIZER_LISTS
3461 
3462 
3463 
3464 
3465 
3466 
3467 
3468 
3469 template<class T, class A>
3470 inline
3471 Marray<T, A>::Marray
3472 (
3473  std::initializer_list<std::size_t> shape,
3474  const T& value,
3475  const CoordinateOrder& coordinateOrder,
3476  const allocator_type& allocator
3477 
3478 )
3479 : dataAllocator_(allocator)
3480 {
3481  std::size_t size = std::accumulate(shape.begin(), shape.end(),
3482  static_cast<std::size_t>(1), std::multiplies<std::size_t>());
3483  marray_detail::Assert(MARRAY_NO_ARG_TEST || size != 0);
3484  base::assign(shape.begin(), shape.end(), dataAllocator_.allocate(size),
3485  coordinateOrder, coordinateOrder, allocator);
3486  std::fill(this->data_, this->data_+size, value);
3487  testInvariant();
3488 }
3489 #endif
3490 
3493 template<class T, class A>
3494 inline
3495 Marray<T, A>::~Marray()
3496 {
3497  dataAllocator_.deallocate(this->data_, this->size());
3498 }
3499 
3512 template<class T, class A>
3513 
3514 Marray<T, A>&
3515 Marray<T, A>::operator=
3516 (
3517  const Marray<T, A>& in
3518 )
3519 {
3520  testInvariant();
3521  if(!MARRAY_NO_ARG_TEST) {
3522  in.testInvariant();
3523  }
3524  if(this != &in) { // no self-assignment
3525  // copy data
3526  if(in.data_ == 0) { // un-initialized
3527  // free
3528  dataAllocator_.deallocate(this->data_, this->size());
3529  this->data_ = 0;
3530  }
3531  else {
3532  if(this->size() != in.size()) {
3533  // re-alloc
3534  dataAllocator_.deallocate(this->data_, this->size());
3535  this->data_ = dataAllocator_.allocate(in.size());
3536  }
3537  // copy data
3538  memcpy(this->data_, in.data_, in.size() * sizeof(T));
3539  }
3540  this->geometry_ = in.geometry_;
3541  }
3542  testInvariant();
3543  return *this;
3544 }
3545 
3560 template<class T, class A>
3561 template<class TLocal, bool isConstLocal, class ALocal>
3562 Marray<T, A>&
3563 Marray<T, A>::operator=
3564 (
3565  const View<TLocal, isConstLocal, ALocal>& in
3566 )
3567 {
3568  if(!MARRAY_NO_ARG_TEST) {
3569  in.testInvariant();
3570  }
3571  if( (void*)(this) != (void*)(&in) ) { // no self-assignment
3572  if(in.data_ == 0) {
3573  dataAllocator_.deallocate(this->data_, this->size());
3574  this->data_ = 0;
3575  this->geometry_ = in.geometry_;
3576  }
3577  else if(this->overlaps(in)) {
3578  Marray<T, A> m = in; // temporary copy
3579  (*this) = m;
3580  }
3581  else {
3582  // re-alloc memory if necessary
3583  if(this->size() != in.size()) {
3584  dataAllocator_.deallocate(this->data_, this->size());
3585  this->data_ = dataAllocator_.allocate(in.size());
3586  }
3587 
3588  // copy geometry
3589  this->geometry_.resize(in.dimension());
3590  for(std::size_t j=0; j<in.dimension(); ++j) {
3591  this->geometry_.shape(j) = in.geometry_.shape(j);
3592  this->geometry_.shapeStrides(j) = in.geometry_.shapeStrides(j);
3593  this->geometry_.strides(j) = in.geometry_.shapeStrides(j); // !
3594  }
3595  this->geometry_.size() = in.size();
3596  this->geometry_.isSimple() = true;
3597  this->geometry_.coordinateOrder() = in.coordinateOrder();
3598 
3599  // copy data
3600  if(in.isSimple() && marray_detail::IsEqual<T, TLocal>::type) {
3601  memcpy(this->data_, in.data_, (in.size())*sizeof(T));
3602  }
3603  else if(in.dimension() == 1)
3604  marray_detail::OperateHelperBinary<1, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3605  else if(in.dimension() == 2)
3606  marray_detail::OperateHelperBinary<2, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3607  else if(in.dimension() == 3)
3608  marray_detail::OperateHelperBinary<3, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3609  else if(in.dimension() == 4)
3610  marray_detail::OperateHelperBinary<4, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3611  else if(in.dimension() == 5)
3612  marray_detail::OperateHelperBinary<5, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3613  else if(in.dimension() == 6)
3614  marray_detail::OperateHelperBinary<6, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3615  else if(in.dimension() == 7)
3616  marray_detail::OperateHelperBinary<7, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3617  else if(in.dimension() == 8)
3618  marray_detail::OperateHelperBinary<8, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3619  else if(in.dimension() == 9)
3620  marray_detail::OperateHelperBinary<9, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3621  else if(in.dimension() == 10)
3622  marray_detail::OperateHelperBinary<10, marray_detail::Assign<T, TLocal>, T, TLocal, isConstLocal, A, ALocal>::operate(*this, in, marray_detail::Assign<T, TLocal>(), this->data_, &in(0));
3623  else {
3624  typename View<TLocal, isConstLocal, ALocal>::const_iterator it = in.begin();
3625  for(std::size_t j=0; j<this->size(); ++j, ++it) {
3626  this->data_[j] = static_cast<T>(*it);
3627  }
3628  }
3629  }
3630  }
3631  testInvariant();
3632  return *this;
3633 }
3634 
3641 template<class T, class A>
3642 inline Marray<T, A>&
3643 Marray<T, A>::operator=
3644 (
3645  const T& value
3646 )
3647 {
3648  marray_detail::Assert(MARRAY_NO_DEBUG || this->data_ != 0);
3649  for(std::size_t j=0; j<this->size(); ++j) {
3650  this->data_[j] = value;
3651  }
3652  return *this;
3653 }
3654 
3655 template<class T, class A>
3656 template<class E, class Te>
3657 inline Marray<T, A>&
3658 Marray<T, A>::operator=
3659 (
3660  const ViewExpression<E, Te>& expression
3661 )
3662 {
3663  if(expression.overlaps(*this)) {
3664  Marray<T, A> m(expression); // temporary copy
3665  (*this) = m; // recursive call
3666  }
3667  else {
3668  // re-allocate memory (if necessary)
3669  if(this->size() != expression.size()) {
3670  dataAllocator_.deallocate(this->data_, this->size());
3671  this->data_ = dataAllocator_.allocate(expression.size());
3672  }
3673 
3674  // copy geometry
3675  this->geometry_.resize(expression.dimension());
3676  for(std::size_t j=0; j<expression.dimension(); ++j) {
3677  this->geometry_.shape(j) = expression.shape(j);
3678  }
3679  this->geometry_.size() = expression.size();
3680  this->geometry_.isSimple() = true;
3681  this->geometry_.coordinateOrder() = expression.coordinateOrder();
3682  if(this->geometry_.dimension() != 0) {
3683  marray_detail::stridesFromShape(this->geometry_.shapeBegin(), this->geometry_.shapeEnd(),
3684  this->geometry_.shapeStridesBegin(), this->geometry_.coordinateOrder());
3685  marray_detail::stridesFromShape(this->geometry_.shapeBegin(), this->geometry_.shapeEnd(),
3686  this->geometry_.stridesBegin(), this->geometry_.coordinateOrder());
3687  }
3688 
3689  // copy data
3690  marray_detail::operate(*this, expression, marray_detail::Assign<T, Te>());
3691  }
3692  return *this;
3693 }
3694 
3695 template<class T, class A>
3696 template<bool SKIP_INITIALIZATION, class ShapeIterator>
3697 inline void
3698 Marray<T, A>::resizeHelper
3699 (
3700  ShapeIterator begin,
3701  ShapeIterator end,
3702  const T& value
3703 )
3704 {
3705  testInvariant();
3706  // compute size
3707  std::vector<std::size_t> newShape;
3708  std::size_t newSize = 1;
3709  for(ShapeIterator it = begin; it != end; ++it) {
3710  std::size_t x = static_cast<std::size_t>(*it);
3711  marray_detail::Assert(MARRAY_NO_ARG_TEST || x > 0);
3712  newShape.push_back(x);
3713  newSize *= x;
3714  }
3715  // allocate new
3716  value_type* newData = dataAllocator_.allocate(newSize);
3717  if(!SKIP_INITIALIZATION) {
3718  for(std::size_t j=0; j<newSize; ++j) {
3719  newData[j] = value;
3720  }
3721  }
3722  // copy old data in region of overlap
3723  if(this->data_ != 0) {
3724  if(newSize == 1 || this->dimension() == 0) {
3725  newData[0] = this->data_[0];
3726  }
3727  else {
3728  std::vector<std::size_t> base1(this->dimension());
3729  std::vector<std::size_t> base2(newShape.size());
3730  std::vector<std::size_t> shape1(this->dimension(), 1);
3731  std::vector<std::size_t> shape2(newShape.size(), 1);
3732  for(std::size_t j=0; j<std::min(this->dimension(), newShape.size()); ++j) {
3733  shape1[j] = std::min(this->shape(j), newShape[j]);
3734  shape2[j] = shape1[j];
3735  }
3736  View<T, true, A> view1;
3737  this->constView(base1.begin(), shape1.begin(), view1);
3738  View<T, false, A> viewT(newShape.begin(), newShape.end(),
3739  newData, this->coordinateOrder(),
3740  this->coordinateOrder());
3741  View<T, false, A> view2;
3742  viewT.view(base2.begin(), shape2.begin(), view2);
3743  view1.squeeze();
3744  view2.squeeze();
3745  view2 = view1; // copy
3746  }
3747  dataAllocator_.deallocate(this->data_, this->size());
3748  this->data_ = 0;
3749  }
3750  base::assign(begin, end, newData, this->geometry_.coordinateOrder(),
3751  this->geometry_.coordinateOrder());
3752  testInvariant();
3753 }
3754 
3762 template<class T, class A>
3763 template<class ShapeIterator>
3764 void
3765 Marray<T, A>::resize
3766 (
3767  ShapeIterator begin,
3768  ShapeIterator end,
3769  const T& value
3770 )
3771 {
3772  resizeHelper<false>(begin, end, value);
3773 }
3774 
3782 template<class T, class A>
3783 template<class ShapeIterator>
3784 void
3785 Marray<T, A>::resize
3786 (
3787  const InitializationSkipping& is,
3788  ShapeIterator begin,
3789  ShapeIterator end
3790 )
3791 {
3792  resizeHelper<true>(begin, end);
3793 }
3794 
3795 #ifdef HAVE_CPP11_INITIALIZER_LISTS
3796 
3797 
3798 
3799 
3800 
3801 template<class T, class A>
3802 inline void
3803 Marray<T, A>::resize
3804 (
3805  std::initializer_list<std::size_t> shape,
3806  const T& value
3807 )
3808 {
3809  resizeHelper<false>(shape.begin(), shape.end(), value);
3810 }
3811 
3817 template<class T, class A>
3818 inline void
3819 Marray<T, A>::resize
3820 (
3821  const InitializationSkipping& si,
3822  std::initializer_list<std::size_t> shape
3823 )
3824 {
3825  resizeHelper<true>(shape.begin(), shape.end());
3826 }
3827 #endif
3828 
3831 template<class T, class A>
3832 inline void
3833 Marray<T, A>::testInvariant() const
3834 {
3835  View<T, false, A>::testInvariant();
3836  marray_detail::Assert(MARRAY_NO_DEBUG || this->geometry_.isSimple());
3837 }
3838 
3839 // iterator implementation
3840 
3843 template<class T, bool isConst, class A>
3844 inline void
3845 Iterator<T, isConst, A>::testInvariant() const
3846 {
3847  if(!MARRAY_NO_DEBUG) {
3848  if(view_ == 0) {
3849  marray_detail::Assert(coordinates_.size() == 0
3850  && index_ == 0
3851  && pointer_ == 0);
3852  }
3853  else {
3854  if(view_->size() == 0) { // un-initialized view
3855  marray_detail::Assert(coordinates_.size() == 0
3856  && index_ == 0
3857  && pointer_ == 0);
3858  }
3859  else { // initialized view
3860  marray_detail::Assert(index_ >= 0 && index_ <= view_->size());
3861  if(index_ == view_->size()) { // end iterator
3862  marray_detail::Assert(pointer_ == &((*view_)(view_->size()-1)) + 1);
3863  }
3864  else {
3865  marray_detail::Assert(pointer_ == &((*view_)(index_)));
3866  }
3867  if(!view_->isSimple()) {
3868  marray_detail::Assert(coordinates_.size() == view_->dimension());
3869  if(index_ == view_->size()) { // end iterator
3870  if(view_->coordinateOrder() == LastMajorOrder) {
3871  marray_detail::Assert(coordinates_[0] == view_->shape(0));
3872  for(std::size_t j=1; j<coordinates_.size(); ++j) {
3873  marray_detail::Assert(coordinates_[j] == view_->shape(j)-1);
3874  }
3875  }
3876  else { // FirstMajorOrder
3877  std::size_t d = view_->dimension() - 1;
3878  marray_detail::Assert(coordinates_[d] == view_->shape(d));
3879  for(std::size_t j=0; j<d; ++j) {
3880  marray_detail::Assert(coordinates_[j] == view_->shape(j)-1);
3881  }
3882  }
3883  }
3884  else {
3885  std::vector<std::size_t> testCoord(coordinates_.size());
3886  view_->indexToCoordinates(index_, testCoord.begin());
3887  for(std::size_t j=0; j<coordinates_.size(); ++j) {
3888  marray_detail::Assert(coordinates_[j] == testCoord[j]);
3889  }
3890  }
3891  }
3892  }
3893  }
3894  }
3895 }
3896 
3898 template<class T, bool isConst, class A>
3899 inline Iterator<T, isConst, A>::Iterator()
3900 : view_(0),
3901  pointer_(0),
3902  index_(0),
3903  coordinates_(std::vector<std::size_t>())
3904 {
3905  testInvariant();
3906 }
3907 
3913 template<class T, bool isConst, class A>
3914 inline Iterator<T, isConst, A>::Iterator
3915 (
3916  const View<T, true, A>& view,
3917  const std::size_t index
3918 )
3919 : view_(&view),
3920  pointer_(0),
3921  index_(index),
3922  coordinates_(std::vector<std::size_t>(view.dimension()))
3923  // Note for developers: If isConst==false, the construction view_(&view)
3924  // fails due to incompatible types. This is intended because it should
3925  // not be possible to construct a mutable iterator on constant data.
3926 {
3927  if(view.size() == 0) { // un-initialized view
3928  marray_detail::Assert(MARRAY_NO_ARG_TEST || index == 0);
3929  }
3930  else {
3931  if(view.isSimple()) {
3932  marray_detail::Assert(MARRAY_NO_ARG_TEST || index <= view.size());
3933  pointer_ = &view(0) + index;
3934  }
3935  else {
3936  if(index >= view.size()) { // end iterator
3937  if(view_->coordinateOrder() == LastMajorOrder) {
3938  coordinates_[0] = view.shape(0);
3939  for(std::size_t j=1; j<view.dimension(); ++j) {
3940  coordinates_[j] = view.shape(j)-1;
3941  }
3942  }
3943  else { // FirstMajorOrder
3944  std::size_t d = view_->dimension() - 1;
3945  coordinates_[d] = view.shape(d);
3946  for(std::size_t j=0; j<d; ++j) {
3947  coordinates_[j] = view.shape(j)-1;
3948  }
3949  }
3950  pointer_ = &view(view.size()-1) + 1;
3951  }
3952  else {
3953  view.indexToCoordinates(index, coordinates_.begin());
3954  pointer_ = &view(index);
3955  }
3956  }
3957  }
3958  testInvariant();
3959 }
3960 
3966 template<class T, bool isConst, class A>
3967 inline Iterator<T, isConst, A>::Iterator
3968 (
3969  const View<T, false, A>& view,
3970  const std::size_t index
3971 )
3972 : view_(reinterpret_cast<view_pointer>(&view)),
3973  pointer_(0),
3974  index_(index),
3975  coordinates_(std::vector<std::size_t>(view.dimension()))
3976  // Note for developers: If isConst==true, the construction
3977  // view_(reinterpret_cast<view_pointer>(&view)) works as well.
3978  // This is intended because it should be possible to construct
3979  // a constant iterator on mutable data.
3980 {
3981  if(view.size() == 0) { // un-initialized view
3982  marray_detail::Assert(MARRAY_NO_ARG_TEST || index == 0);
3983  }
3984  else {
3985  if(view.isSimple()) {
3986  marray_detail::Assert(MARRAY_NO_ARG_TEST || index <= view.size());
3987  pointer_ = &view(0) + index;
3988  }
3989  else {
3990  if(index >= view.size()) { // end iterator
3991  if(view_->coordinateOrder() == LastMajorOrder) {
3992  coordinates_[0] = view.shape(0);
3993  for(std::size_t j=1; j<view.dimension(); ++j) {
3994  coordinates_[j] = view.shape(j)-1;
3995  }
3996  }
3997  else { // FirstMajorOrder
3998  std::size_t d = view_->dimension() - 1;
3999  coordinates_[d] = view.shape(d);
4000  for(std::size_t j=0; j<d; ++j) {
4001  coordinates_[j] = view.shape(j)-1;
4002  }
4003  }
4004  pointer_ = &view(view.size()-1) + 1;
4005  }
4006  else {
4007  view.indexToCoordinates(index, coordinates_.begin());
4008  pointer_ = &view(index);
4009  }
4010  }
4011  }
4012  testInvariant();
4013 }
4014 
4020 template<class T, bool isConst, class A>
4021 inline Iterator<T, isConst, A>::Iterator
4022 (
4023  View<T, false, A>& view,
4024  const std::size_t index
4025 )
4026 : view_(reinterpret_cast<view_pointer>(&view)),
4027  pointer_(0),
4028  index_(index),
4029  coordinates_(std::vector<std::size_t>(view.dimension()))
4030  // Note for developers: If isConst==true, the construction
4031  // view_(reinterpret_cast<view_pointer>(&view)) works as well.
4032  // This is intended because it should be possible to construct
4033  // a constant iterator on mutable data.
4034 {
4035  if(view.size() == 0) { // un-initialized view
4036  marray_detail::Assert(MARRAY_NO_ARG_TEST || index == 0);
4037  }
4038  else {
4039  if(view.isSimple()) {
4040  marray_detail::Assert(MARRAY_NO_ARG_TEST || index <= view.size());
4041  pointer_ = &view(0) + index;
4042  }
4043  else {
4044  if(index >= view.size()) { // end iterator
4045  if(view_->coordinateOrder() == LastMajorOrder) {
4046  coordinates_[0] = view.shape(0);
4047  for(std::size_t j=1; j<view.dimension(); ++j) {
4048  coordinates_[j] = view.shape(j)-1;
4049  }
4050  }
4051  else { // FirstMajorOrder
4052  std::size_t d = view_->dimension() - 1;
4053  coordinates_[d] = view.shape(d);
4054  for(std::size_t j=0; j<d; ++j) {
4055  coordinates_[j] = view.shape(j)-1;
4056  }
4057  }
4058  pointer_ = &view(view.size()-1) + 1;
4059  }
4060  else {
4061  view.indexToCoordinates(index, coordinates_.begin());
4062  pointer_ = &view(index);
4063  }
4064  }
4065  }
4066  testInvariant();
4067 }
4068 
4071 template<class T, bool isConst, class A>
4072 inline Iterator<T, isConst, A>::Iterator
4073 (
4074  const Iterator<T, false, A>& in
4075 )
4076 : view_(view_pointer(in.view_)),
4077  pointer_(pointer(in.pointer_)),
4078  index_(in.index_),
4079  coordinates_(in.coordinates_)
4080 {
4081  testInvariant();
4082 }
4083 
4086 template<class T, bool isConst, class A>
4087 inline typename Iterator<T, isConst, A>::reference
4088 Iterator<T, isConst, A>::operator*() const
4089 {
4090  marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && index_ < view_->size()));
4091  return *pointer_;
4092 }
4093 
4096 template<class T, bool isConst, class A>
4097 inline typename Iterator<T, isConst, A>::pointer
4098 Iterator<T, isConst, A>::operator->() const
4099 {
4100  marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && index_ < view_->size()));
4101  return pointer_;
4102 }
4103 
4106 template<class T, bool isConst, class A>
4107 inline typename Iterator<T, isConst, A>::reference
4108 Iterator<T, isConst, A>::operator[]
4109 (
4110  const std::size_t x
4111 ) const
4112 {
4113  marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && x+index_ < view_->size()));
4114  return (*view_)(x+index_);
4115 }
4116 
4117 template<class T, bool isConst, class A>
4118 inline Iterator<T, isConst, A>&
4119 Iterator<T, isConst, A>::operator+=
4120 (
4121  const difference_type& x
4122 )
4123 {
4124  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4125  if(index_ < view_->size()) { // view initialized and iterator not at the end
4126  if(index_ + x < view_->size()) {
4127  index_ += x;
4128  if(view_->isSimple()) {
4129  pointer_ += x;
4130  }
4131  else {
4132  pointer_ = &((*view_)(index_));
4133  view_->indexToCoordinates(index_, coordinates_.begin());
4134  }
4135  }
4136  else {
4137  // set to end iterator
4138  index_ = view_->size();
4139  if(view_->isSimple()) {
4140  pointer_ = &(*view_)(0) + view_->size();
4141  }
4142  else {
4143  pointer_ = (&(*view_)(view_->size()-1)) + 1;
4144  view_->indexToCoordinates(view_->size()-1, coordinates_.begin());
4145  if(view_->coordinateOrder() == LastMajorOrder) {
4146  ++coordinates_[0];
4147  }
4148  else { // FirstMajorOrder
4149  ++coordinates_[view_->dimension()-1];
4150  }
4151  }
4152  }
4153  }
4154  testInvariant();
4155  return *this;
4156 }
4157 
4158 template<class T, bool isConst, class A>
4159 inline Iterator<T, isConst, A>&
4160 Iterator<T, isConst, A>::operator-=
4161 (
4162  const difference_type& x
4163 )
4164 {
4165  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4166  marray_detail::Assert(MARRAY_NO_ARG_TEST || static_cast<difference_type>(index_) >= x);
4167  index_ -= x;
4168  if(view_->isSimple()) {
4169  pointer_ -= x;
4170  }
4171  else {
4172  pointer_ = &((*view_)(index_));
4173  view_->indexToCoordinates(index_, coordinates_.begin());
4174  }
4175  testInvariant();
4176  return *this;
4177 }
4178 
4181 template<class T, bool isConst, class A>
4182 inline Iterator<T, isConst, A>&
4183 Iterator<T, isConst, A>::operator++()
4184 {
4185  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4186  if(index_ < view_->size()) { // view initialized and iterator not at the end
4187  ++index_;
4188  if(view_->isSimple()) {
4189  ++pointer_;
4190  }
4191  else {
4192  if(index_ < view_->size()) {
4193  if(view_->coordinateOrder() == LastMajorOrder) {
4194  for(std::size_t j=0; j<coordinates_.size(); ++j) {
4195  if(coordinates_[j] == view_->shape(j)-1) {
4196  pointer_ -= view_->strides(j) * coordinates_[j];
4197  coordinates_[j] = 0;
4198  }
4199  else {
4200  pointer_ += view_->strides(j);
4201  ++coordinates_[j];
4202  break;
4203  }
4204  }
4205  }
4206  else { // FirstMajorOrder
4207  std::size_t j = coordinates_.size() - 1;
4208  for(;;) {
4209  if(coordinates_[j] == view_->shape(j)-1) {
4210  pointer_ -= view_->strides(j) * coordinates_[j];
4211  coordinates_[j] = 0;
4212  }
4213  else {
4214  pointer_ += view_->strides(j);
4215  ++coordinates_[j];
4216  break;
4217  }
4218  if(j == 0) {
4219  break;
4220  }
4221  else {
4222  --j;
4223  }
4224  }
4225  }
4226  }
4227  else {
4228  // set to end iterator
4229  pointer_ = &((*view_)(view_->size()-1)) + 1;
4230  if(view_->coordinateOrder() == LastMajorOrder) {
4231  ++coordinates_[0];
4232  }
4233  else { // FirstMajorOrder
4234  ++coordinates_[view_->dimension()-1];
4235  }
4236  }
4237  }
4238  }
4239  testInvariant();
4240  return *this;
4241 }
4242 
4245 template<class T, bool isConst, class A>
4246 inline Iterator<T, isConst, A>&
4247 Iterator<T, isConst, A>::operator--()
4248 {
4249  marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && index_ > 0));
4250  --index_;
4251  if(view_->isSimple()) {
4252  --pointer_;
4253  }
4254  else {
4255  if(index_ == view_->size()) {
4256  // decrement from end iterator
4257  --pointer_;
4258  if(view_->coordinateOrder() == LastMajorOrder) {
4259  --coordinates_[0];
4260  }
4261  else { // FirstMajorOrder
4262  --coordinates_[view_->dimension() - 1];
4263  }
4264  }
4265  else {
4266  if(view_->coordinateOrder() == LastMajorOrder) {
4267  for(std::size_t j=0; j<coordinates_.size(); ++j) {
4268  if(coordinates_[j] == 0) {
4269  coordinates_[j] = view_->shape(j)-1;
4270  pointer_ += view_->strides(j) * coordinates_[j];
4271  }
4272  else {
4273  pointer_ -= view_->strides(j);
4274  --coordinates_[j];
4275  break;
4276  }
4277  }
4278  }
4279  else { // FirstMajorOrder
4280  std::size_t j = view_->dimension() - 1;
4281  for(;;) {
4282  if(coordinates_[j] == 0) {
4283  coordinates_[j] = view_->shape(j)-1;
4284  pointer_ += view_->strides(j) * coordinates_[j];
4285  }
4286  else {
4287  pointer_ -= view_->strides(j);
4288  --coordinates_[j];
4289  break;
4290  }
4291  if(j == 0) {
4292  break;
4293  }
4294  else {
4295  --j;
4296  }
4297  }
4298  }
4299  }
4300  }
4301  testInvariant();
4302  return *this;
4303 }
4304 
4307 template<class T, bool isConst, class A>
4308 inline Iterator<T, isConst, A>
4309 Iterator<T, isConst, A>::operator++(int)
4310 {
4311  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4312  Iterator<T, isConst, A> copy = *this;
4313  ++(*this);
4314  return copy;
4315 }
4316 
4319 template<class T, bool isConst, class A>
4320 inline Iterator<T, isConst, A>
4321 Iterator<T, isConst, A>::operator--(int)
4322 {
4323  marray_detail::Assert(MARRAY_NO_DEBUG || (view_ != 0 && index_ > 0));
4324  Iterator<T, isConst, A> copy = *this;
4325  --(*this);
4326  return copy;
4327 }
4328 
4329 template<class T, bool isConst, class A>
4330 inline Iterator<T, isConst, A>
4331 Iterator<T, isConst, A>::operator+
4332 (
4333  const difference_type& x
4334 ) const
4335 {
4336  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4337  Iterator<T, isConst, A> tmp = *this;
4338  tmp += x;
4339  return tmp;
4340 }
4341 
4342 template<class T, bool isConst, class A>
4343 inline Iterator<T, isConst, A>
4344 Iterator<T, isConst, A>::operator-
4345 (
4346  const difference_type& x
4347 ) const
4348 {
4349  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4350  Iterator<T, isConst, A> tmp = *this;
4351  tmp -= x;
4352  return tmp;
4353 }
4354 
4355 template<class T, bool isConst, class A>
4356 template<bool isConstLocal>
4357 inline typename Iterator<T, isConst, A>::difference_type
4358 Iterator<T, isConst, A>::operator-
4359 (
4360  const Iterator<T, isConstLocal, A>& it
4361 ) const
4362 {
4363  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4364  marray_detail::Assert(MARRAY_NO_ARG_TEST || it.view_ != 0);
4365  return difference_type(index_)-difference_type(it.index_);
4366 }
4367 
4368 template<class T, bool isConst, class A>
4369 template<bool isConstLocal>
4370 inline bool
4371 Iterator<T, isConst, A>::operator==
4372 (
4373  const Iterator<T, isConstLocal, A>& it
4374 ) const
4375 {
4376  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4377  marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && (void*)it.view_ == (void*)view_));
4378  return index_ == it.index_;
4379 }
4380 
4381 template<class T, bool isConst, class A>
4382 template<bool isConstLocal>
4383 inline bool
4384 Iterator<T, isConst, A>::operator!=
4385 (
4386  const Iterator<T, isConstLocal, A>& it
4387 ) const
4388 {
4389  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4390  marray_detail::Assert(MARRAY_NO_ARG_TEST || it.view_ != 0);
4391  marray_detail::Assert(MARRAY_NO_ARG_TEST ||
4392  static_cast<const void*>(it.view_) == static_cast<const void*>(view_));
4393  return index_ != it.index_;
4394 }
4395 
4396 template<class T, bool isConst, class A>
4397 template<bool isConstLocal>
4398 inline bool
4399 Iterator<T, isConst, A>::operator<
4400 (
4401  const Iterator<T, isConstLocal, A>& it
4402 ) const
4403 {
4404  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4405  marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && it.view_ == view_));
4406  return(index_ < it.index_);
4407 }
4408 
4409 template<class T, bool isConst, class A>
4410 template<bool isConstLocal>
4411 inline bool
4412 Iterator<T, isConst, A>::operator>
4413 (
4414  const Iterator<T, isConstLocal, A>& it
4415 ) const
4416 {
4417  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4418  marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && it.view_ == view_));
4419  return(index_ > it.index_);
4420 }
4421 
4422 template<class T, bool isConst, class A>
4423 template<bool isConstLocal>
4424 inline bool
4425 Iterator<T, isConst, A>::operator<=
4426 (
4427  const Iterator<T, isConstLocal, A>& it
4428 ) const
4429 {
4430  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4431  marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && it.view_ == view_));
4432  return(index_ <= it.index_);
4433 }
4434 
4435 template<class T, bool isConst, class A>
4436 template<bool isConstLocal>
4437 inline bool
4438 Iterator<T, isConst, A>::operator>=
4439 (
4440  const Iterator<T, isConstLocal, A>& it
4441 ) const
4442 {
4443  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4444  marray_detail::Assert(MARRAY_NO_ARG_TEST || (it.view_ != 0 && it.view_ == view_));
4445  return(index_ >= it.index_);
4446 }
4447 
4452 template<class T, bool isConst, class A>
4453 inline bool
4454 Iterator<T, isConst, A>::hasMore() const
4455 {
4456  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4457  return index_ < view_->size();
4458 }
4459 
4464 template<class T, bool isConst, class A>
4465 inline std::size_t
4466 Iterator<T, isConst, A>::index() const
4467 {
4468  return index_;
4469 }
4470 
4476 template<class T, bool isConst, class A>
4477 template<class CoordinateIterator>
4478 inline void
4479 Iterator<T, isConst, A>::coordinate
4480 (
4481  CoordinateIterator it
4482 ) const
4483 {
4484  marray_detail::Assert(MARRAY_NO_DEBUG || view_ != 0);
4485  marray_detail::Assert(MARRAY_NO_ARG_TEST || index_ < view_->size());
4486  if(view_->isSimple()) {
4487  view_->indexToCoordinates(index_, it);
4488  }
4489  else {
4490  for(std::size_t j=0; j<coordinates_.size(); ++j, ++it) {
4491  *it = coordinates_[j];
4492  }
4493  }
4494 }
4495 
4496 // implementation of expression templates
4497 
4499 template<class E, class T>
4500 class ViewExpression {
4501 public:
4502  typedef E expression_type;
4503  typedef T value_type;
4504 
4505  const std::size_t dimension() const
4506  { return static_cast<const E&>(*this).dimension(); }
4507  const std::size_t size() const
4508  { return static_cast<const E&>(*this).size(); }
4509  const std::size_t shape(const std::size_t j) const
4510  { return static_cast<const E&>(*this).shape(j); }
4511  const std::size_t* shapeBegin() const
4512  { return static_cast<const E&>(*this).shapeBegin(); }
4513  const std::size_t* shapeEnd() const
4514  { return static_cast<const E&>(*this).shapeEnd(); }
4515  template<class Tv, bool isConst, class A>
4516  bool overlaps(const View<Tv, isConst, A>& v) const
4517  { return static_cast<const E&>(*this).overlaps(v); }
4518  const CoordinateOrder& coordinateOrder() const
4519  { return static_cast<const E&>(*this).coordinateOrder(); }
4520  const bool isSimple() const
4521  { return static_cast<const E&>(*this).isSimple(); }
4522  template<class Accessor>
4523  const T& operator()(Accessor it) const
4524  { return static_cast<const E&>(*this)(it); }
4525  const T& operator()(const std::size_t c0, const std::size_t c1) const
4526  { return static_cast<const E&>(*this)(c0, c1); }
4527  const T& operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2) const
4528  { return static_cast<const E&>(*this)(c0, c1, c2); }
4529  const T& operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2, const std::size_t c3) const
4530  { return static_cast<const E&>(*this)(c0, c1, c2, c3); }
4531  const T& operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2, const std::size_t c3, const std::size_t c4) const
4532  { return static_cast<const E&>(*this)(c0, c1, c2, c3, c4); }
4533  const T& operator[](const std::size_t offset) const
4534  { return static_cast<const E&>(*this)[offset]; }
4535  operator E&()
4536  { return static_cast<E&>(*this); }
4537  operator E const&() const
4538  { return static_cast<const E&>(*this); }
4539 
4540  // \cond suppress_doxygen
4541  class ExpressionIterator {
4542  public:
4543  ExpressionIterator(const ViewExpression<E, T>& expression)
4544  : expression_(expression), // cast!
4545  data_(&expression_(0)),
4546  offset_(0)
4547  {}
4548  void incrementCoordinate(const std::size_t coordinateIndex)
4549  { offset_ += expression_.strides(coordinateIndex); }
4550  void resetCoordinate(const std::size_t coordinateIndex)
4551  { offset_ -= expression_.strides(coordinateIndex)
4552 
4553  * (expression_.shape(coordinateIndex) - 1); }
4554  const T& operator*() const
4555  { // return expression_[offset_];
4556  // would require making this nested class a friend of View
4557  // which in turn would require a forward declaration of
4558  // this class. work around:
4559  return data_[offset_]; }
4560  private:
4561  const E& expression_;
4562  const T* data_;
4563  std::size_t offset_;
4564  };
4565  // \endcond suppress_doxygen
4566 };
4567 
4568 // \cond suppress_doxygen
4569 template<class E, class T, class UnaryFunctor>
4570 class UnaryViewExpression
4571 : public ViewExpression<UnaryViewExpression<E, T, UnaryFunctor>, T>
4572 {
4573 public:
4574  typedef E expression_type;
4575  typedef T value_type;
4576 
4577  UnaryViewExpression(const ViewExpression<E, T>& e)
4578  : e_(e), // cast!
4579  unaryFunctor_(UnaryFunctor())
4580  {}
4581  const std::size_t dimension() const
4582  { return e_.dimension(); }
4583  const std::size_t size() const
4584  { return e_.size(); }
4585  const std::size_t shape(const std::size_t j) const
4586  { return e_.shape(j); }
4587  const std::size_t* shapeBegin() const
4588  { return e_.shapeBegin(); }
4589  const std::size_t* shapeEnd() const
4590  { return e_.shapeEnd(); }
4591  template<class Tv, bool isConst, class A>
4592  bool overlaps(const View<Tv, isConst, A>& v) const
4593  { return e_.overlaps(v); }
4594  const CoordinateOrder& coordinateOrder() const
4595  { return e_.coordinateOrder(); }
4596  const bool isSimple() const
4597  { return e_.isSimple(); }
4598  template<class Accessor>
4599  const T operator()(Accessor it) const
4600  { return unaryFunctor_(e_(it)); }
4601  const T operator()(const std::size_t c0, const std::size_t c1) const
4602  { return unaryFunctor_(e_(c0, c1)); }
4603  const T operator()(const std::size_t c0, const std::size_t c1,const std::size_t c2) const
4604  { return unaryFunctor_(e_(c0, c1, c2)); }
4605  const T operator()(const std::size_t c0, const std::size_t c1,const std::size_t c2, const std::size_t c3) const
4606  { return unaryFunctor_(e_(c0, c1, c2, c3)); }
4607  const T operator()(const std::size_t c0, const std::size_t c1,const std::size_t c2, const std::size_t c3, const std::size_t c4) const
4608  { return unaryFunctor_(e_(c0, c1, c2, c3, c4)); }
4609  const T operator[](const std::size_t offset) const
4610  { return unaryFunctor_(e_[offset]); }
4611 
4612  class ExpressionIterator {
4613  public:
4614  ExpressionIterator(const UnaryViewExpression<E, T, UnaryFunctor>& expression)
4615  : unaryFunctor_(expression.unaryFunctor_),
4616  iterator_(expression.e_)
4617  {}
4618  void incrementCoordinate(const std::size_t coordinateIndex)
4619  { iterator_.incrementCoordinate(coordinateIndex); }
4620  void resetCoordinate(const std::size_t coordinateIndex)
4621  { iterator_.resetCoordinate(coordinateIndex); }
4622  const T operator*() const
4623  { return unaryFunctor_(*iterator_); }
4624  private:
4625  UnaryFunctor unaryFunctor_;
4626  typename E::ExpressionIterator iterator_;
4627  };
4628 
4629 private:
4630  const E& e_;
4631  UnaryFunctor unaryFunctor_;
4632 };
4633 
4634 template<class E1, class T1, class E2, class T2, class BinaryFunctor>
4635 class BinaryViewExpression
4636 : public ViewExpression<BinaryViewExpression<E1, T1, E2, T2, BinaryFunctor>,
4637  typename marray_detail::PromoteType<T1, T2>::type>
4638 {
4639 public:
4640  typedef E1 expression_type_1;
4641  typedef E2 expression_type_2;
4642  typedef T1 value_type_1;
4643  typedef T2 value_type_2;
4644  typedef typename marray_detail::PromoteType<T1, T2>::type value_type;
4645  typedef BinaryFunctor functor_type;
4646  typedef ViewExpression<BinaryViewExpression<E1, T1, E2, T2, BinaryFunctor>,
4647  value_type> base;
4648 
4649  BinaryViewExpression(const ViewExpression<E1, T1>& e1,
4650  const ViewExpression<E2, T2>& e2)
4651  : e1_(e1), e2_(e2), // cast!
4652  binaryFunctor_(BinaryFunctor())
4653  {
4654  if(!MARRAY_NO_DEBUG) {
4655  marray_detail::Assert(e1_.size() != 0 && e2_.size() != 0);
4656  marray_detail::Assert(e1_.dimension() == e2_.dimension());
4657  for(std::size_t j=0; j<e1_.dimension(); ++j) {
4658  marray_detail::Assert(e1_.shape(j) == e2_.shape(j));
4659  }
4660  }
4661  }
4662  const std::size_t dimension() const
4663  { return e1_.dimension() < e2_.dimension() ? e2_.dimension() : e1_.dimension(); }
4664  const std::size_t size() const
4665  { return e1_.size() < e2_.size() ? e2_.size() : e1_.size(); }
4666  const std::size_t shape(const std::size_t j) const
4667  { return e1_.dimension() < e2_.dimension() ? e2_.shape(j) : e1_.shape(j); }
4668  const std::size_t* shapeBegin() const
4669  { return e1_.dimension() < e2_.dimension() ? e2_.shapeBegin() : e1_.shapeBegin(); }
4670  const std::size_t* shapeEnd() const
4671  { return e1_.dimension() < e2_.dimension() ? e2_.shapeEnd() : e1_.shapeEnd(); }
4672  template<class Tv, bool isConst, class A>
4673  bool overlaps(const View<Tv, isConst, A>& v) const
4674  { return e1_.overlaps(v) || e2_.overlaps(v); }
4675  const CoordinateOrder& coordinateOrder() const
4676  { return e1_.coordinateOrder(); }
4677  const bool isSimple() const
4678  { return e1_.isSimple() && e2_.isSimple()
4679  && e1_.coordinateOrder() == e2_.coordinateOrder(); }
4680  template<class Accessor>
4681  const value_type operator()(Accessor it) const
4682  { return binaryFunctor_(e1_(it), e2_(it)); }
4683  const value_type operator()(const std::size_t c0, const std::size_t c1) const
4684  { return binaryFunctor_(e1_(c0, c1), e2_(c0, c1)); }
4685  const value_type operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2) const
4686  { return binaryFunctor_(e1_(c0, c1, c2), e2_(c0, c1, c2)); }
4687  const value_type operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2, const std::size_t c3) const
4688  { return binaryFunctor_(e1_(c0, c1, c2, c3), e2_(c0, c1, c2, c3)); }
4689  const value_type operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2, const std::size_t c3, const std::size_t c4) const
4690  { return binaryFunctor_(e1_(c0, c1, c2, c3, c4), e2_(c0, c1, c2, c3, c4)); }
4691  const value_type operator[](const std::size_t offset) const
4692  { return binaryFunctor_(e1_[offset], e2_[offset]); }
4693 
4694  class ExpressionIterator {
4695  public:
4696  ExpressionIterator(const BinaryViewExpression<E1, T1, E2, T2, BinaryFunctor>& expression)
4697  : binaryFunctor_(expression.binaryFunctor_),
4698  iterator1_(expression.e1_),
4699  iterator2_(expression.e2_)
4700  {}
4701  void incrementCoordinate(const std::size_t coordinateIndex)
4702  { iterator1_.incrementCoordinate(coordinateIndex);
4703  iterator2_.incrementCoordinate(coordinateIndex); }
4704  void resetCoordinate(const std::size_t coordinateIndex)
4705  { iterator1_.resetCoordinate(coordinateIndex);
4706  iterator2_.resetCoordinate(coordinateIndex); }
4707  const value_type operator*() const
4708  { return binaryFunctor_(*iterator1_, *iterator2_); }
4709  private:
4710  BinaryFunctor binaryFunctor_;
4711  typename E1::ExpressionIterator iterator1_;
4712  typename E2::ExpressionIterator iterator2_;
4713  };
4714 
4715 private:
4716  const expression_type_1& e1_;
4717  const expression_type_2& e2_;
4718  BinaryFunctor binaryFunctor_;
4719 };
4720 
4721 template<class E, class T, class S, class BinaryFunctor>
4722 class BinaryViewExpressionScalarFirst
4723 : public ViewExpression<BinaryViewExpressionScalarFirst<E, T, S, BinaryFunctor>,
4724  typename marray_detail::PromoteType<T, S>::type> {
4725 public:
4726  typedef E expression_type;
4727  typedef T value_type_1;
4728  typedef S scalar_type;
4729  typedef typename marray_detail::PromoteType<T, S>::type value_type;
4730  typedef BinaryFunctor functor_type;
4731  typedef ViewExpression<BinaryViewExpressionScalarFirst<E, T, S, BinaryFunctor>,
4732  value_type> base;
4733 
4734  BinaryViewExpressionScalarFirst(const ViewExpression<E, T>& e,
4735  const scalar_type& scalar)
4736  : e_(e), // cast!
4737  scalar_(scalar), binaryFunctor_(BinaryFunctor())
4738  { }
4739  const std::size_t dimension() const
4740  { return e_.dimension(); }
4741  const std::size_t size() const
4742  { return e_.size(); }
4743  const std::size_t shape(const std::size_t j) const
4744  { return e_.shape(j); }
4745  const std::size_t* shapeBegin() const
4746  { return e_.shapeBegin(); }
4747  const std::size_t* shapeEnd() const
4748  { return e_.shapeEnd(); }
4749  template<class Tv, bool isConst, class A>
4750  bool overlaps(const View<Tv, isConst, A>& v) const
4751  { return e_.overlaps(v); }
4752  const CoordinateOrder& coordinateOrder() const
4753  { return e_.coordinateOrder(); }
4754  const bool isSimple() const
4755  { return e_.isSimple(); }
4756  template<class Accessor>
4757  const value_type operator()(Accessor it) const
4758  { return binaryFunctor_(scalar_, e_(it)); }
4759  const value_type operator()(const std::size_t c0, const std::size_t c1) const
4760  { return binaryFunctor_(scalar_, e_(c0, c1)); }
4761  const value_type operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2) const
4762  { return binaryFunctor_(scalar_, e_(c0, c1, c2)); }
4763  const value_type operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2, const std::size_t c3) const
4764  { return binaryFunctor_(scalar_, e_(c0, c1, c2, c3)); }
4765  const value_type operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2, const std::size_t c3, const std::size_t c4) const
4766  { return binaryFunctor_(scalar_, e_(c0, c1, c2, c3, c4)); }
4767  const value_type operator[](const std::size_t offset) const
4768  { return binaryFunctor_(scalar_, e_[offset]); }
4769 
4770  class ExpressionIterator {
4771  public:
4772  ExpressionIterator(const BinaryViewExpressionScalarFirst<E, T, S, BinaryFunctor>& expression)
4773  : binaryFunctor_(expression.binaryFunctor_),
4774  scalar_(expression.scalar_),
4775  iterator_(expression.e_)
4776  {}
4777  void incrementCoordinate(const std::size_t coordinateIndex)
4778  { iterator_.incrementCoordinate(coordinateIndex); }
4779  void resetCoordinate(const std::size_t coordinateIndex)
4780  { iterator_.resetCoordinate(coordinateIndex); }
4781  const T operator*() const
4782  { return binaryFunctor_(scalar_, *iterator_); }
4783  private:
4784  BinaryFunctor binaryFunctor_;
4785  const typename BinaryViewExpressionScalarFirst<E, T, S, BinaryFunctor>::scalar_type& scalar_;
4786  typename E::ExpressionIterator iterator_;
4787  };
4788 
4789 private:
4790  const expression_type& e_;
4791  const scalar_type scalar_;
4792  BinaryFunctor binaryFunctor_;
4793 };
4794 
4795 template<class E, class T, class S, class BinaryFunctor>
4796 class BinaryViewExpressionScalarSecond
4797 : public ViewExpression<BinaryViewExpressionScalarSecond<E, T, S, BinaryFunctor>,
4798  typename marray_detail::PromoteType<T, S>::type> {
4799 public:
4800  typedef T value_type_1;
4801  typedef E expression_type;
4802  typedef S scalar_type;
4803  typedef typename marray_detail::PromoteType<T, S>::type value_type;
4804  typedef BinaryFunctor functor_type;
4805  typedef ViewExpression<BinaryViewExpressionScalarSecond<E, T, S, BinaryFunctor>,
4806  value_type> base;
4807 
4808  BinaryViewExpressionScalarSecond(const ViewExpression<E, T>& e,
4809  const scalar_type& scalar)
4810  : e_(e), // cast!
4811  scalar_(scalar), binaryFunctor_(BinaryFunctor())
4812  { }
4813  const std::size_t dimension() const
4814  { return e_.dimension(); }
4815  const std::size_t size() const
4816  { return e_.size(); }
4817  const std::size_t shape(const std::size_t j) const
4818  { return e_.shape(j); }
4819  const std::size_t* shapeBegin() const
4820  { return e_.shapeBegin(); }
4821  const std::size_t* shapeEnd() const
4822  { return e_.shapeEnd(); }
4823  template<class Tv, bool isConst, class A>
4824  bool overlaps(const View<Tv, isConst, A>& v) const
4825  { return e_.overlaps(v); }
4826  const CoordinateOrder& coordinateOrder() const
4827  { return e_.coordinateOrder(); }
4828  const bool isSimple() const
4829  { return e_.isSimple(); }
4830  template<class Accessor>
4831  const value_type operator()(Accessor it) const
4832  { return binaryFunctor_(e_(it), scalar_); }
4833  const value_type operator()(const std::size_t c0, const std::size_t c1) const
4834  { return binaryFunctor_(e_(c0, c1), scalar_); }
4835  const value_type operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2) const
4836  { return binaryFunctor_(e_(c0, c1, c2), scalar_); }
4837  const value_type operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2, const std::size_t c3) const
4838  { return binaryFunctor_(e_(c0, c1, c2, c3), scalar_); }
4839  const value_type operator()(const std::size_t c0, const std::size_t c1, const std::size_t c2, const std::size_t c3, const std::size_t c4) const
4840  { return binaryFunctor_(e_(c0, c1, c2, c3, c4), scalar_); }
4841  const value_type operator[](const std::size_t offset) const
4842  { return binaryFunctor_(e_[offset], scalar_); }
4843 
4844  class ExpressionIterator {
4845  public:
4846  ExpressionIterator(const BinaryViewExpressionScalarSecond<E, T, S, BinaryFunctor>& expression)
4847  : binaryFunctor_(expression.binaryFunctor_),
4848  scalar_(expression.scalar_),
4849  iterator_(expression.e_)
4850  {}
4851  void incrementCoordinate(const std::size_t coordinateIndex)
4852  { iterator_.incrementCoordinate(coordinateIndex); }
4853  void resetCoordinate(const std::size_t coordinateIndex)
4854  { iterator_.resetCoordinate(coordinateIndex); }
4855  const T operator*() const
4856  { return binaryFunctor_(*iterator_, scalar_); }
4857  private:
4858  BinaryFunctor binaryFunctor_;
4859  const typename BinaryViewExpressionScalarSecond<E, T, S, BinaryFunctor>::scalar_type& scalar_;
4860  typename E::ExpressionIterator iterator_;
4861  };
4862 
4863 private:
4864  const expression_type& e_;
4865  const scalar_type scalar_;
4866  BinaryFunctor binaryFunctor_;
4867 };
4868 // \endcond suppress_doxygen
4869 
4870 // implementation of marray_detail
4871 
4872 // \cond suppress_doxygen
4873 namespace marray_detail {
4874 
4875 template<class A>
4876 class Geometry
4877 {
4878 public:
4879  typedef typename A::template rebind<std::size_t>::other allocator_type;
4880 
4881  Geometry(const allocator_type& = allocator_type());
4882  Geometry(const std::size_t, const CoordinateOrder& = defaultOrder,
4883  const std::size_t = 0, const bool = true,
4884  const allocator_type& = allocator_type());
4885  template<class ShapeIterator>
4886  Geometry(ShapeIterator, ShapeIterator,
4887  const CoordinateOrder& = defaultOrder,
4888  const CoordinateOrder& = defaultOrder,
4889  const allocator_type& = allocator_type());
4890  template<class ShapeIterator, class StrideIterator>
4891  Geometry(ShapeIterator, ShapeIterator, StrideIterator,
4892  const CoordinateOrder& = defaultOrder,
4893  const allocator_type& = allocator_type());
4894  Geometry(const Geometry<A>&);
4895  ~Geometry();
4896 
4897  Geometry<A>& operator=(const Geometry<A>&);
4898 
4899  void resize(const std::size_t dimension);
4900  const std::size_t dimension() const;
4901  const std::size_t shape(const std::size_t) const;
4902  std::size_t& shape(const std::size_t);
4903  const std::size_t shapeStrides(const std::size_t) const;
4904  std::size_t& shapeStrides(const std::size_t);
4905  const std::size_t strides(const std::size_t) const;
4906  std::size_t& strides(const std::size_t);
4907  const std::size_t* shapeBegin() const;
4908  std::size_t* shapeBegin();
4909  const std::size_t* shapeEnd() const;
4910  std::size_t* shapeEnd();
4911  const std::size_t* shapeStridesBegin() const;
4912  std::size_t* shapeStridesBegin();
4913  const std::size_t* shapeStridesEnd() const;
4914  std::size_t* shapeStridesEnd();
4915  const std::size_t* stridesBegin() const;
4916  std::size_t* stridesBegin();
4917  const std::size_t* stridesEnd() const;
4918  std::size_t* stridesEnd();
4919  const std::size_t size() const;
4920  std::size_t& size();
4921  const CoordinateOrder& coordinateOrder() const;
4922  CoordinateOrder& coordinateOrder();
4923  const bool isSimple() const;
4924  void updateSimplicity();
4925  bool& isSimple();
4926 
4927 private:
4928  allocator_type allocator_;
4929  std::size_t* shape_;
4930  std::size_t* shapeStrides_;
4931  // Intended redundancy: shapeStrides_ could be
4932  // computed from shape_ and coordinateOrder_
4933  std::size_t* strides_;
4934  std::size_t dimension_;
4935  std::size_t size_;
4936  // intended redundancy: size_ could be computed from shape_
4937  CoordinateOrder coordinateOrder_;
4938  bool isSimple_;
4939  // simple array: an array which is unstrided (i.e. the strides
4940  // equal the shape strides), cf. the function testInvariant of
4941  // View for the formal definition.
4942 };
4943 
4944 template<class A>
4945 inline
4946 Geometry<A>::Geometry
4947 (
4948  const typename Geometry<A>::allocator_type& allocator
4949 )
4950 : allocator_(allocator),
4951  shape_(0),
4952  shapeStrides_(0),
4953  strides_(0),
4954  dimension_(0),
4955  size_(0),
4956  coordinateOrder_(defaultOrder),
4957  isSimple_(true)
4958 {
4959 }
4960 
4961 template<class A>
4962 inline
4963 Geometry<A>::Geometry
4964 (
4965  const Geometry<A>& g
4966 )
4967 : allocator_(g.allocator_),
4968  shape_(g.dimension_==0 ? 0 : allocator_.allocate(g.dimension_*3)),
4969  shapeStrides_(shape_ + g.dimension_),
4970  strides_(shapeStrides_ + g.dimension_),
4971  dimension_(g.dimension_),
4972  size_(g.size_),
4973  coordinateOrder_(g.coordinateOrder_),
4974  isSimple_(g.isSimple_)
4975 {
4976  /*
4977  for(std::size_t j=0; j<dimension_; ++j) {
4978  shape_[j] = g.shape_[j];
4979  shapeStrides_[j] = g.shapeStrides_[j];
4980  strides_[j] = g.strides_[j];
4981  }
4982  */
4983  memcpy(shape_, g.shape_, (dimension_*3)*sizeof(std::size_t));
4984 }
4985 
4986 template<class A>
4987 inline
4988 Geometry<A>::Geometry
4989 (
4990  const std::size_t dimension,
4991  const CoordinateOrder& order,
4992  const std::size_t size,
4993  const bool isSimple,
4994  const typename Geometry<A>::allocator_type& allocator
4995 )
4996 : allocator_(allocator),
4997  shape_(allocator_.allocate(dimension*3)),
4998  shapeStrides_(shape_+dimension),
4999  strides_(shapeStrides_+dimension),
5000  dimension_(dimension),
5001  size_(size),
5002  coordinateOrder_(order),
5003  isSimple_(isSimple)
5004 {
5005 }
5006 
5007 template<class A>
5008 template<class ShapeIterator>
5009 inline
5010 Geometry<A>::Geometry
5011 (
5012  ShapeIterator begin,
5013  ShapeIterator end,
5014  const CoordinateOrder& externalCoordinateOrder,
5015  const CoordinateOrder& internalCoordinateOrder,
5016  const typename Geometry<A>::allocator_type& allocator
5017 )
5018 : allocator_(allocator),
5019  shape_(allocator_.allocate(std::distance(begin, end) * 3)),
5020  shapeStrides_(shape_ + std::distance(begin, end)),
5021  strides_(shapeStrides_ + std::distance(begin, end)),
5022  dimension_(std::distance(begin, end)),
5023  size_(1),
5024  coordinateOrder_(internalCoordinateOrder),
5025  isSimple_(true)
5026 {
5027  if(dimension_ != 0) { // if the array is not a scalar
5028  isSimple_ = (externalCoordinateOrder == internalCoordinateOrder);
5029  for(std::size_t j=0; j<dimension(); ++j, ++begin) {
5030  const std::size_t s = static_cast<std::size_t>(*begin);
5031  shape(j) = s;
5032  size() *= s;
5033  }
5034  stridesFromShape(shapeBegin(), shapeEnd(), stridesBegin(),
5035  externalCoordinateOrder);
5036  stridesFromShape(shapeBegin(), shapeEnd(), shapeStridesBegin(),
5037  internalCoordinateOrder);
5038  }
5039 }
5040 
5041 template<class A>
5042 template<class ShapeIterator, class StrideIterator>
5043 inline
5044 Geometry<A>::Geometry
5045 (
5046  ShapeIterator begin,
5047  ShapeIterator end,
5048  StrideIterator it,
5049  const CoordinateOrder& internalCoordinateOrder,
5050  const typename Geometry<A>::allocator_type& allocator
5051 )
5052 : allocator_(allocator),
5053  shape_(allocator_.allocate(std::distance(begin, end) * 3)),
5054  shapeStrides_(shape_ + std::distance(begin, end)),
5055  strides_(shapeStrides_ + std::distance(begin, end)),
5056  dimension_(std::distance(begin, end)),
5057  size_(1),
5058  coordinateOrder_(internalCoordinateOrder),
5059  isSimple_(true)
5060 {
5061  if(dimension() != 0) {
5062  for(std::size_t j=0; j<dimension(); ++j, ++begin, ++it) {
5063  const std::size_t s = static_cast<std::size_t>(*begin);
5064  shape(j) = s;
5065  size() *= s;
5066  strides(j) = *it;
5067  }
5068  stridesFromShape(shapeBegin(), shapeEnd(), shapeStridesBegin(),
5069  internalCoordinateOrder);
5070  updateSimplicity();
5071  }
5072 }
5073 
5074 template<class A>
5075 inline
5076 Geometry<A>::~Geometry()
5077 {
5078  allocator_.deallocate(shape_, dimension_*3);
5079 }
5080 
5081 template<class A>
5082 inline Geometry<A>&
5083 Geometry<A>::operator=
5084 (
5085  const Geometry<A>& g
5086 )
5087 {
5088  if(&g != this) { // no self-assignment
5089  if(g.dimension_ != dimension_) {
5090  allocator_.deallocate(shape_, dimension_*3);
5091  dimension_ = g.dimension_;
5092  shape_ = allocator_.allocate(dimension_*3);
5093  shapeStrides_ = shape_+dimension_;
5094  strides_ = shapeStrides_+dimension_;
5095  dimension_ = g.dimension_;
5096  }
5097  /*
5098  for(std::size_t j=0; j<dimension_; ++j) {
5099  shape_[j] = g.shape_[j];
5100  shapeStrides_[j] = g.shapeStrides_[j];
5101  strides_[j] = g.strides_[j];
5102  }
5103  */
5104  memcpy(shape_, g.shape_, (dimension_*3)*sizeof(std::size_t));
5105  size_ = g.size_;
5106  coordinateOrder_ = g.coordinateOrder_;
5107  isSimple_ = g.isSimple_;
5108  }
5109  return *this;
5110 }
5111 
5112 template<class A>
5113 inline void
5114 Geometry<A>::resize
5115 (
5116  const std::size_t dimension
5117 )
5118 {
5119  if(dimension != dimension_) {
5120  std::size_t* newShape = allocator_.allocate(dimension*3);
5121  std::size_t* newShapeStrides = newShape + dimension;
5122  std::size_t* newStrides = newShapeStrides + dimension;
5123  for(std::size_t j=0; j<( (dimension < dimension_) ? dimension : dimension_); ++j) {
5124  // save existing entries
5125  newShape[j] = shape(j);
5126  newShapeStrides[j] = shapeStrides(j);
5127  newStrides[j] = strides(j);
5128  }
5129  allocator_.deallocate(shape_, dimension_*3);
5130  shape_ = newShape;
5131  shapeStrides_ = newShapeStrides;
5132  strides_ = newStrides;
5133  dimension_ = dimension;
5134  }
5135 }
5136 
5137 template<class A>
5138 inline const std::size_t
5139 Geometry<A>::dimension() const
5140 {
5141  return dimension_;
5142 }
5143 
5144 template<class A>
5145 inline const std::size_t
5146 Geometry<A>::shape(const std::size_t j) const
5147 {
5148  Assert(MARRAY_NO_DEBUG || j<dimension_);
5149  return shape_[j];
5150 }
5151 
5152 template<class A>
5153 inline std::size_t&
5154 Geometry<A>::shape(const std::size_t j)
5155 {
5156  Assert(MARRAY_NO_DEBUG || j<dimension_);
5157  return shape_[j];
5158 }
5159 
5160 template<class A>
5161 inline const std::size_t
5162 Geometry<A>::shapeStrides
5163 (
5164  const std::size_t j
5165 ) const
5166 {
5167  Assert(MARRAY_NO_DEBUG || j<dimension_);
5168  return shapeStrides_[j];
5169 }
5170 
5171 template<class A>
5172 inline std::size_t&
5173 Geometry<A>::shapeStrides
5174 (
5175  const std::size_t j
5176 )
5177 {
5178  Assert(MARRAY_NO_DEBUG || j<dimension_);
5179  return shapeStrides_[j];
5180 }
5181 
5182 template<class A>
5183 inline const std::size_t
5184 Geometry<A>::strides
5185 (
5186  const std::size_t j
5187 ) const
5188 {
5189  Assert(MARRAY_NO_DEBUG || j<dimension_);
5190  return strides_[j];
5191 }
5192 
5193 template<class A>
5194 inline std::size_t&
5195 Geometry<A>::strides
5196 (
5197  const std::size_t j
5198 )
5199 {
5200  Assert(MARRAY_NO_DEBUG || j<dimension_);
5201  return strides_[j];
5202 }
5203 
5204 template<class A>
5205 inline const std::size_t*
5206 Geometry<A>::shapeBegin() const
5207 {
5208  return shape_;
5209 }
5210 
5211 template<class A>
5212 inline std::size_t*
5213 Geometry<A>::shapeBegin()
5214 {
5215  return shape_;
5216 }
5217 
5218 template<class A>
5219 inline const std::size_t*
5220 Geometry<A>::shapeEnd() const
5221 {
5222  return shape_ + dimension_;
5223 }
5224 
5225 template<class A>
5226 inline std::size_t*
5227 Geometry<A>::shapeEnd()
5228 {
5229  return shape_ + dimension_;
5230 }
5231 
5232 template<class A>
5233 inline const std::size_t*
5234 Geometry<A>::shapeStridesBegin() const
5235 {
5236  return shapeStrides_;
5237 }
5238 
5239 template<class A>
5240 inline std::size_t*
5241 Geometry<A>::shapeStridesBegin()
5242 {
5243  return shapeStrides_;
5244 }
5245 
5246 template<class A>
5247 inline const std::size_t*
5248 Geometry<A>::shapeStridesEnd() const
5249 {
5250  return shapeStrides_ + dimension_;
5251 }
5252 
5253 template<class A>
5254 inline std::size_t*
5255 Geometry<A>::shapeStridesEnd()
5256 {
5257  return shapeStrides_ + dimension_;
5258 }
5259 
5260 template<class A>
5261 inline const std::size_t*
5262 Geometry<A>::stridesBegin() const
5263 {
5264  return strides_;
5265 }
5266 
5267 template<class A>
5268 inline std::size_t*
5269 Geometry<A>::stridesBegin()
5270 {
5271  return strides_;
5272 }
5273 
5274 template<class A>
5275 inline const std::size_t*
5276 Geometry<A>::stridesEnd() const
5277 {
5278  return strides_ + dimension_;
5279 }
5280 
5281 template<class A>
5282 inline std::size_t*
5283 Geometry<A>::stridesEnd()
5284 {
5285  return strides_ + dimension_;
5286 }
5287 
5288 template<class A>
5289 inline const std::size_t
5290 Geometry<A>::size() const
5291 {
5292  return size_;
5293 }
5294 
5295 template<class A>
5296 inline std::size_t&
5297 Geometry<A>::size()
5298 {
5299  return size_;
5300 }
5301 
5302 template<class A>
5303 inline const CoordinateOrder&
5304 Geometry<A>::coordinateOrder() const
5305 {
5306  return coordinateOrder_;
5307 }
5308 
5309 template<class A>
5310 inline CoordinateOrder&
5311 Geometry<A>::coordinateOrder()
5312 {
5313  return coordinateOrder_;
5314 }
5315 
5316 template<class A>
5317 inline const bool
5318 Geometry<A>::isSimple() const
5319 {
5320  return isSimple_;
5321 }
5322 
5323 template<class A>
5324 inline bool&
5325 Geometry<A>::isSimple()
5326 {
5327  return isSimple_;
5328 }
5329 
5330 template<class A>
5331 inline void
5332 Geometry<A>::updateSimplicity()
5333 {
5334  for(std::size_t j=0; j<dimension(); ++j) {
5335  if(shapeStrides(j) != strides(j)) {
5336  isSimple_ = false;
5337  return;
5338  }
5339  }
5340  isSimple_ = true;
5341  // a 0-dimensional geometry is simple
5342 }
5343 
5344 template<class ShapeIterator, class StridesIterator>
5345 inline void
5346 stridesFromShape
5347 (
5348  ShapeIterator begin,
5349  ShapeIterator end,
5350  StridesIterator strideBegin,
5351  const CoordinateOrder& coordinateOrder
5352 )
5353 {
5354  Assert(MARRAY_NO_ARG_TEST || std::distance(begin, end) != 0);
5355  std::size_t dimension = std::distance(begin, end);
5356  ShapeIterator shapeIt;
5357  StridesIterator strideIt;
5358  if(coordinateOrder == FirstMajorOrder) {
5359  shapeIt = begin + (dimension-1);
5360  strideIt = strideBegin + (dimension-1);
5361  *strideIt = 1;
5362  for(std::size_t j=1; j<dimension; ++j) {
5363  std::size_t tmp = *strideIt;
5364  --strideIt;
5365  (*strideIt) = tmp * (*shapeIt);
5366  --shapeIt;
5367  }
5368  }
5369  else {
5370  shapeIt = begin;
5371  strideIt = strideBegin;
5372  *strideIt = 1;
5373  for(std::size_t j=1; j<dimension; ++j) {
5374  std::size_t tmp = *strideIt;
5375  ++strideIt;
5376  (*strideIt) = tmp * (*shapeIt);
5377  ++shapeIt;
5378  }
5379  }
5380 }
5381 
5382 template<unsigned short N, class Functor, class T, class A>
5383 struct OperateHelperUnary
5384 {
5385  static inline void operate
5386  (
5387  View<T, false, A>& v,
5388  Functor f,
5389  T* data
5390  )
5391  {
5392  for(std::size_t j=0; j<v.shape(N-1); ++j) {
5393  OperateHelperUnary<N-1, Functor, T, A>::operate(v, f, data);
5394  data += v.strides(N-1);
5395  }
5396  data -= v.shape(N-1) * v.strides(N-1);
5397  }
5398 };
5399 
5400 
5401 template<class Functor, class T, class A>
5402 struct OperateHelperUnary<0, Functor, T, A>
5403 {
5404  static inline void operate
5405  (
5406  View<T, false, A>& v,
5407  Functor f,
5408  T* data
5409  )
5410  {
5411  f(*data);
5412  }
5413 };
5414 
5415 template<unsigned short N, class Functor, class T1, class T2, class A>
5416 struct OperateHelperBinaryScalar
5417 {
5418  static inline void operate
5419  (
5420  View<T1, false, A>& v,
5421  const T2& x,
5422  Functor f,
5423  T1* data
5424  )
5425  {
5426  for(std::size_t j=0; j<v.shape(N-1); ++j) {
5427  OperateHelperBinaryScalar<N-1, Functor, T1, T2, A>::operate(
5428  v, x, f, data);
5429  data += v.strides(N-1);
5430  }
5431  data -= v.shape(N-1) * v.strides(N-1);
5432  }
5433 };
5434 
5435 template<class Functor, class T1, class T2, class A>
5436 struct OperateHelperBinaryScalar<0, Functor, T1, T2, A>
5437 {
5438  static inline void operate
5439  (
5440  View<T1, false, A>& v,
5441  const T2& x,
5442  Functor f,
5443  T1* data
5444  )
5445  {
5446  f(*data, x);
5447  }
5448 };
5449 
5450 template<unsigned short N, class Functor, class T1, class T2,
5451  bool isConst, class A1, class A2>
5452 struct OperateHelperBinary
5453 {
5454  static inline void operate
5455  (
5456  View<T1, false, A1>& v,
5457  const View<T2, isConst, A2>& w,
5458  Functor f,
5459  T1* data1,
5460  const T2* data2
5461  )
5462  {
5463  for(std::size_t j=0; j<v.shape(N-1); ++j) {
5464  OperateHelperBinary<N-1, Functor, T1, T2, isConst, A1, A2>::operate(
5465  v, w, f, data1, data2);
5466  data1 += v.strides(N-1);
5467  data2 += w.strides(N-1);
5468  }
5469  data1 -= v.shape(N-1) * v.strides(N-1);
5470  data2 -= w.shape(N-1) * w.strides(N-1);
5471  }
5472 };
5473 
5474 template<class Functor, class T1, class T2, bool isConst, class A1, class A2>
5475 struct OperateHelperBinary<0, Functor, T1, T2, isConst, A1, A2>
5476 {
5477  static inline void operate
5478  (
5479  View<T1, false, A1>& v,
5480  const View<T2, isConst, A2>& w,
5481  Functor f,
5482  T1* data1,
5483  const T2* data2
5484  )
5485  {
5486  f(*data1, *data2);
5487  }
5488 };
5489 
5490 template<class TFrom, class TTo, class AFrom, class ATo>
5491 struct AssignmentOperatorHelper<false, TFrom, TTo, AFrom, ATo>
5492 {
5493  // from constant to mutable
5494  //
5495  // here, 'to' must be initialized (which is asserted) because
5496  // otherwise, the pointer to.data_ to mutable data would have
5497  // to be initialized with the pointer from.data_ to constant
5498  // data which we don't do.
5499  //
5500  static void execute
5501  (
5502  const View<TFrom, true, AFrom>& from,
5503  View<TTo, false, ATo>& to
5504  )
5505  {
5506  typedef typename View<TFrom, true, AFrom>::const_iterator FromIterator;
5507  typedef typename View<TTo, false, ATo>::iterator ToIterator;
5508  if(!MARRAY_NO_ARG_TEST) {
5509  Assert(from.data_ != 0 && from.dimension() == to.dimension());
5510  for(std::size_t j=0; j<from.dimension(); ++j) {
5511  Assert(from.shape(j) == to.shape(j));
5512  }
5513  }
5514  if(from.overlaps(to)) {
5515  Marray<TFrom, AFrom> m = from; // temporary copy
5516  execute(m, to);
5517  }
5518  else if(from.coordinateOrder() == to.coordinateOrder()
5519  && from.isSimple() && to.isSimple()
5520  && IsEqual<TFrom, TTo>::type) {
5521  memcpy(to.data_, from.data_, (from.size())*sizeof(TFrom));
5522  }
5523  else if(from.dimension() == 1)
5524  OperateHelperBinary<1, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5525  else if(from.dimension() == 2)
5526  OperateHelperBinary<2, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5527  else if(from.dimension() == 3)
5528  OperateHelperBinary<3, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5529  else if(from.dimension() == 4)
5530  OperateHelperBinary<4, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5531  else if(from.dimension() == 5)
5532  OperateHelperBinary<5, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5533  else if(from.dimension() == 6)
5534  OperateHelperBinary<6, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5535  else if(from.dimension() == 7)
5536  OperateHelperBinary<7, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5537  else if(from.dimension() == 8)
5538  OperateHelperBinary<8, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5539  else if(from.dimension() == 9)
5540  OperateHelperBinary<9, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5541  else if(from.dimension() == 10)
5542  OperateHelperBinary<10, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5543  else {
5544  FromIterator itFrom = from.begin();
5545  ToIterator itTo = to.begin();
5546  for(; itFrom.hasMore(); ++itFrom, ++itTo) {
5547  *itTo = static_cast<TTo>(*itFrom);
5548  }
5549  }
5550  }
5551 
5556  static void execute
5557  (
5558  const View<TFrom, false, AFrom>& from,
5559  View<TTo, false, ATo>& to
5560  )
5561  {
5562  typedef typename View<TFrom, false, AFrom>::const_iterator FromIterator;
5563  typedef typename View<TTo, false, ATo>::iterator ToIterator;
5564  if(static_cast<const void*>(&from) != static_cast<const void*>(&to)) { // no self-assignment
5565  if(to.data_ == 0) { // if the view 'to' is not initialized
5566  // initialize the view 'to' with source data
5567  Assert(MARRAY_NO_ARG_TEST || sizeof(TTo) == sizeof(TFrom));
5568  to.data_ = static_cast<TTo*>(static_cast<void*>(from.data_)); // copy pointer
5569  to.geometry_ = from.geometry_;
5570  }
5571  else { // if the view 'to' is initialized
5572  if(!MARRAY_NO_ARG_TEST) {
5573  Assert(from.data_ != 0 && from.dimension() == to.dimension());
5574  for(std::size_t j=0; j<from.dimension(); ++j) {
5575  Assert(from.shape(j) == to.shape(j));
5576  }
5577  }
5578  if(from.overlaps(to)) {
5579  Marray<TFrom, AFrom> m = from; // temporary copy
5580  execute(m, to);
5581  }
5582  else if(from.coordinateOrder() == to.coordinateOrder()
5583  && from.isSimple() && to.isSimple()
5584  && IsEqual<TFrom, TTo>::type) {
5585  memcpy(to.data_, from.data_, (from.size())*sizeof(TFrom));
5586  }
5587  else if(from.dimension() == 1)
5588  OperateHelperBinary<1, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5589  else if(from.dimension() == 2)
5590  OperateHelperBinary<2, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5591  else if(from.dimension() == 3)
5592  OperateHelperBinary<3, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5593  else if(from.dimension() == 4)
5594  OperateHelperBinary<4, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5595  else if(from.dimension() == 5)
5596  OperateHelperBinary<5, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5597  else if(from.dimension() == 6)
5598  OperateHelperBinary<6, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5599  else if(from.dimension() == 7)
5600  OperateHelperBinary<7, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5601  else if(from.dimension() == 8)
5602  OperateHelperBinary<8, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5603  else if(from.dimension() == 9)
5604  OperateHelperBinary<9, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5605  else if(from.dimension() == 10)
5606  OperateHelperBinary<10, Assign<TTo, TFrom>, TTo, TFrom, true, ATo, AFrom>::operate(to, from, Assign<TTo, TFrom>(), &to(0), &from(0));
5607  else {
5608  FromIterator itFrom = from.begin();
5609  ToIterator itTo = to.begin();
5610  for(; itFrom.hasMore(); ++itFrom, ++itTo) {
5611  *itTo = static_cast<TTo>(*itFrom);
5612  }
5613  }
5614  }
5615  }
5616  }
5617 };
5618 
5619 template<class TFrom, class TTo, class AFrom, class ATo>
5620 struct AssignmentOperatorHelper<true, TFrom, TTo, AFrom, ATo>
5621 {
5623  template<bool isConstFrom>
5624  static void execute
5625  (
5626  const View<TFrom, isConstFrom, AFrom>& from,
5627  View<TTo, true, ATo>& to
5628  )
5629  {
5630  Assert(MARRAY_NO_ARG_TEST || sizeof(TFrom) == sizeof(TTo));
5631  to.data_ = static_cast<const TTo*>(
5632  static_cast<const void*>(from.data_)); // copy pointer
5633  to.geometry_ = from.geometry_;
5634  }
5635 };
5636 
5637 template<>
5638 struct AccessOperatorHelper<true>
5639 {
5640  // access by scalar index
5641  template<class T, class U, bool isConst, class A>
5642  static typename View<T, isConst, A>::reference
5643  execute(const View<T, isConst, A>& v, const U& index)
5644  {
5645  v.testInvariant();
5646  Assert(MARRAY_NO_DEBUG || v.data_ != 0);
5647  Assert(MARRAY_NO_DEBUG || v.dimension() != 0 || index == 0);
5648  std::size_t offset;
5649  v.indexToOffset(index, offset);
5650  return v.data_[offset];
5651  }
5652 };
5653 
5654 template<>
5655 struct AccessOperatorHelper<false>
5656 {
5657  // access by iterator
5658  template<class T, class U, bool isConst, class A>
5659  static typename View<T, isConst, A>::reference
5660  execute(const View<T, isConst, A>& v, U it)
5661  {
5662  v.testInvariant();
5663  Assert(MARRAY_NO_DEBUG || v.data_ != 0);
5664  Assert(MARRAY_NO_DEBUG || v.dimension() != 0 || *it == 0);
5665  std::size_t offset;
5666  v.coordinatesToOffset(it, offset);
5667  return v.data_[offset];
5668  }
5669 };
5670 
5671 template<class Functor, class T, class A>
5672 inline void
5673 operate
5674 (
5675  View<T, false, A>& v,
5676  Functor f
5677 )
5678 {
5679  if(v.isSimple()) {
5680  T* data = &v(0);
5681  for(std::size_t j=0; j<v.size(); ++j) {
5682  f(data[j]);
5683  }
5684  }
5685  else if(v.dimension() == 1)
5686  OperateHelperUnary<1, Functor, T, A>::operate(v, f, &v(0));
5687  else if(v.dimension() == 2)
5688  OperateHelperUnary<2, Functor, T, A>::operate(v, f, &v(0));
5689  else if(v.dimension() == 3)
5690  OperateHelperUnary<3, Functor, T, A>::operate(v, f, &v(0));
5691  else if(v.dimension() == 4)
5692  OperateHelperUnary<4, Functor, T, A>::operate(v, f, &v(0));
5693  else if(v.dimension() == 5)
5694  OperateHelperUnary<5, Functor, T, A>::operate(v, f, &v(0));
5695  else if(v.dimension() == 6)
5696  OperateHelperUnary<6, Functor, T, A>::operate(v, f, &v(0));
5697  else if(v.dimension() == 7)
5698  OperateHelperUnary<7, Functor, T, A>::operate(v, f, &v(0));
5699  else if(v.dimension() == 8)
5700  OperateHelperUnary<8, Functor, T, A>::operate(v, f, &v(0));
5701  else if(v.dimension() == 9)
5702  OperateHelperUnary<9, Functor, T, A>::operate(v, f, &v(0));
5703  else if(v.dimension() == 10)
5704  OperateHelperUnary<10, Functor, T, A>::operate(v, f, &v(0));
5705  else {
5706  for(typename View<T, false, A>::iterator it = v.begin(); it.hasMore(); ++it) {
5707  f(*it);
5708  }
5709  }
5710 }
5711 
5712 template<class Functor, class T, class A>
5713 inline void
5714 operate
5715 (
5716  View<T, false, A>& v,
5717  const T& x,
5718  Functor f
5719 )
5720 {
5721  if(v.isSimple()) {
5722  T* data = &v(0);
5723  for(std::size_t j=0; j<v.size(); ++j) {
5724  f(data[j], x);
5725  }
5726  }
5727  else if(v.dimension() == 1)
5728  OperateHelperBinaryScalar<1, Functor, T, T, A>::operate(v, x, f, &v(0));
5729  else if(v.dimension() == 2)
5730  OperateHelperBinaryScalar<2, Functor, T, T, A>::operate(v, x, f, &v(0));
5731  else if(v.dimension() == 3)
5732  OperateHelperBinaryScalar<3, Functor, T, T, A>::operate(v, x, f, &v(0));
5733  else if(v.dimension() == 4)
5734  OperateHelperBinaryScalar<4, Functor, T, T, A>::operate(v, x, f, &v(0));
5735  else if(v.dimension() == 5)
5736  OperateHelperBinaryScalar<5, Functor, T, T, A>::operate(v, x, f, &v(0));
5737  else if(v.dimension() == 6)
5738  OperateHelperBinaryScalar<6, Functor, T, T, A>::operate(v, x, f, &v(0));
5739  else if(v.dimension() == 7)
5740  OperateHelperBinaryScalar<7, Functor, T, T, A>::operate(v, x, f, &v(0));
5741  else if(v.dimension() == 8)
5742  OperateHelperBinaryScalar<8, Functor, T, T, A>::operate(v, x, f, &v(0));
5743  else if(v.dimension() == 9)
5744  OperateHelperBinaryScalar<9, Functor, T, T, A>::operate(v, x, f, &v(0));
5745  else if(v.dimension() == 10)
5746  OperateHelperBinaryScalar<10, Functor, T, T, A>::operate(v, x, f, &v(0));
5747  else {
5748  for(typename View<T, false, A>::iterator it = v.begin(); it.hasMore(); ++it) {
5749  f(*it, x);
5750  }
5751  }
5752 }
5753 
5754 template<class Functor, class T1, class T2, bool isConst, class A>
5755 inline void
5756 operate
5757 (
5758  View<T1, false, A>& v,
5759  const View<T2, isConst, A>& w,
5760  Functor f
5761 )
5762 {
5763  if(!MARRAY_NO_ARG_TEST) {
5764  Assert(v.size() != 0 && w.size() != 0);
5765  Assert(w.dimension() == 0 || v.dimension() == w.dimension());
5766  if(w.dimension() != 0) {
5767  for(std::size_t j=0; j<v.dimension(); ++j) {
5768  Assert(v.shape(j) == w.shape(j));
5769  }
5770  }
5771  }
5772  if(w.dimension() == 0) {
5773  T2 x = w(0);
5774  if(v.isSimple()) {
5775  T1* dataV = &v(0);
5776  for(std::size_t j=0; j<v.size(); ++j) {
5777  f(dataV[j], x);
5778  }
5779  }
5780  else if(v.dimension() == 1)
5781  OperateHelperBinaryScalar<1, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5782  else if(v.dimension() == 2)
5783  OperateHelperBinaryScalar<2, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5784  else if(v.dimension() == 3)
5785  OperateHelperBinaryScalar<3, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5786  else if(v.dimension() == 4)
5787  OperateHelperBinaryScalar<4, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5788  else if(v.dimension() == 5)
5789  OperateHelperBinaryScalar<5, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5790  else if(v.dimension() == 6)
5791  OperateHelperBinaryScalar<6, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5792  else if(v.dimension() == 7)
5793  OperateHelperBinaryScalar<7, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5794  else if(v.dimension() == 8)
5795  OperateHelperBinaryScalar<8, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5796  else if(v.dimension() == 9)
5797  OperateHelperBinaryScalar<9, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5798  else if(v.dimension() == 10)
5799  OperateHelperBinaryScalar<10, Functor, T1, T2, A>::operate(v, x, f, &v(0));
5800  else {
5801  for(typename View<T1, false>::iterator it = v.begin(); it.hasMore(); ++it) {
5802  f(*it, x);
5803  }
5804  }
5805  }
5806  else {
5807  if(v.overlaps(w)) {
5808  Marray<T2> m = w; // temporary copy
5809  operate(v, m, f); // recursive call
5810  }
5811  else {
5812  if(v.coordinateOrder() == w.coordinateOrder()
5813  && v.isSimple() && w.isSimple()) {
5814  T1* dataV = &v(0);
5815  const T2* dataW = &w(0);
5816  for(std::size_t j=0; j<v.size(); ++j) {
5817  f(dataV[j], dataW[j]);
5818  }
5819  }
5820  else if(v.dimension() == 1)
5821  OperateHelperBinary<1, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5822  else if(v.dimension() == 2)
5823  OperateHelperBinary<2, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5824  else if(v.dimension() == 3)
5825  OperateHelperBinary<3, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5826  else if(v.dimension() == 4)
5827  OperateHelperBinary<4, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5828  else if(v.dimension() == 5)
5829  OperateHelperBinary<5, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5830  else if(v.dimension() == 6)
5831  OperateHelperBinary<6, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5832  else if(v.dimension() == 7)
5833  OperateHelperBinary<7, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5834  else if(v.dimension() == 8)
5835  OperateHelperBinary<8, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5836  else if(v.dimension() == 9)
5837  OperateHelperBinary<9, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5838  else if(v.dimension() == 10)
5839  OperateHelperBinary<10, Functor, T1, T2, isConst, A, A>::operate(v, w, f, &v(0), &w(0));
5840  else {
5841  typename View<T1, false>::iterator itV = v.begin();
5842  typename View<T2, isConst>::const_iterator itW = w.begin();
5843  for(; itV.hasMore(); ++itV, ++itW) {
5844  Assert(MARRAY_NO_DEBUG || itW.hasMore());
5845  f(*itV, *itW);
5846  }
5847  Assert(MARRAY_NO_DEBUG || !itW.hasMore());
5848  }
5849  }
5850  }
5851 }
5852 
5853 template<class Functor, class T1, class A, class E, class T2>
5854 inline void operate
5855 (
5856  View<T1, false, A>& v,
5857  const ViewExpression<E, T2>& expression,
5858  Functor f
5859 )
5860 {
5861  const E& e = expression; // cast
5862  if(!MARRAY_NO_DEBUG) {
5863  Assert(v.size() != 0 && e.size() != 0);
5864  Assert(e.dimension() == v.dimension());
5865  if(v.dimension() == 0) {
5866  Assert(v.size() == 1 && e.size() == 1);
5867  }
5868  else {
5869  for(std::size_t j=0; j<v.dimension(); ++j) {
5870  Assert(v.shape(j) == e.shape(j));
5871  }
5872  }
5873  }
5874  if(e.overlaps(v)) {
5875  Marray<T1, A> m(e); // temporary copy
5876  operate(v, m, f);
5877  }
5878  else if(v.dimension() == 0) {
5879  f(v[0], e[0]);
5880  }
5881  else if(v.isSimple() && e.isSimple()
5882  && v.coordinateOrder() == e.coordinateOrder()) {
5883  for(std::size_t j=0; j<v.size(); ++j) {
5884  f(v[j], e[j]);
5885  }
5886  }
5887  else {
5888  // loop unrolling does not improve performance here
5889  typename E::ExpressionIterator itE(e);
5890  std::size_t offsetV = 0;
5891  std::vector<std::size_t> coordinate(v.dimension());
5892  std::size_t maxDimension = v.dimension() - 1;
5893  for(;;) {
5894  f(v[offsetV], *itE);
5895  for(std::size_t j=0; j<v.dimension(); ++j) {
5896  if(coordinate[j]+1 == v.shape(j)) {
5897  if(j == maxDimension) {
5898  return;
5899  }
5900  else {
5901  offsetV -= coordinate[j] * v.strides(j);
5902  itE.resetCoordinate(j);
5903  coordinate[j] = 0;
5904  }
5905  }
5906  else {
5907  offsetV += v.strides(j);
5908  itE.incrementCoordinate(j);
5909  ++coordinate[j];
5910  break;
5911  }
5912  }
5913  }
5914  }
5915 }
5916 
5917 } // namespace marray_detail
5918 // \endcond suppress_doxygen
5919 
5920 } // namespace andres
5921 
5922 #endif
Generated on Thu Feb 18 2016 12:31:34 for Marray by   doxygen 1.8.1.2