blitz  Version 0.9
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
array-impl.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 /***************************************************************************
3  * blitz/array-impl.h Definition of the Array<P_numtype, N_rank> class
4  *
5  * $Id: array-impl.h,v 1.25 2005/10/13 23:46:43 julianc Exp $
6  *
7  * Copyright (C) 1997-2001 Todd Veldhuizen <tveldhui@oonumerics.org>
8  *
9  * This program is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU General Public License
11  * as published by the Free Software Foundation; either version 2
12  * of the License, or (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * Suggestions: blitz-dev@oonumerics.org
20  * Bugs: blitz-bugs@oonumerics.org
21  *
22  * For more information, please see the Blitz++ Home Page:
23  * http://oonumerics.org/blitz/
24  *
25  ***************************************************************************/
26 
27 /*
28  * Wish list for array classes.
29  * - Arrays whose dimensions are unknown at compile time.
30  * - where()/elsewhere()/elsewhere() as in Dan Quinlan's implementation
31  * - block reduction operations
32  * - conversion to/from matrix & vector
33  * - apply(T func(T))
34  * - apply(T func(const T&))
35  * - apply<T func(T)>
36  */
37 
38 #ifndef BZ_ARRAY_H
39 #define BZ_ARRAY_H
40 
41 #include <blitz/blitz.h>
42 #include <blitz/memblock.h>
43 #include <blitz/range.h>
44 #include <blitz/tinyvec.h>
45 
46 #ifdef BZ_ARRAY_SPACE_FILLING_TRAVERSAL
47 #include <blitz/traversal.h>
48 #endif
49 
50 #include <blitz/indexexpr.h>
51 #include <blitz/prettyprint.h>
52 
53 #include <blitz/array/slice.h> // Subarrays and slicing
54 #include <blitz/array/map.h> // Tensor index notation
55 #include <blitz/array/multi.h> // Multicomponent arrays
56 #include <blitz/array/domain.h> // RectDomain class
57 #include <blitz/array/storage.h> // GeneralArrayStorage
58 
59 
60 BZ_NAMESPACE(blitz)
61 
62 /*
63  * Forward declarations
64  */
65 
66 template<typename T_numtype, int N_rank>
67 class ArrayIterator;
68 
69 template<typename T_numtype, int N_rank>
70 class ConstArrayIterator;
71 
72 template<typename T_numtype, int N_rank>
73 class FastArrayIterator;
74 
75 template<typename P_expr>
76 class _bz_ArrayExpr;
77 
78 template<typename T_array, typename T_index>
79 class IndirectArray;
80 
81 template <typename P_numtype,int N_rank>
82 void swap(Array<P_numtype,N_rank>&,Array<P_numtype,N_rank>&);
83 
84 template <typename P_numtype, int N_rank>
85 void find(Array<TinyVector<int,N_rank>,1>&,const Array<P_numtype,N_rank>&);
86 
87 /*
88  * Declaration of class Array
89  */
90 
91 // NEEDS_WORK: Array should inherit protected from MemoryBlockReference.
92 // To make this work, need to expose MemoryBlockReference::numReferences()
93 // and make Array<P,N2> a friend of Array<P,N> for slicing.
94 
95 template<typename P_numtype, int N_rank>
96 class Array : public MemoryBlockReference<P_numtype>
97 #ifdef BZ_NEW_EXPRESSION_TEMPLATES
99 #endif
100 {
101 
102 private:
104  using T_base::data_;
105  using T_base::changeToNullBlock;
106  using T_base::numReferences;
107 
108 public:
110  // Public Types
112 
113  /*
114  * T_numtype is the numeric type stored in the array.
115  * T_index is a vector type which can be used to access elements
116  * of many-dimensional arrays.
117  * T_array is the array type itself -- Array<T_numtype, N_rank>
118  * T_iterator is a a fast iterator for the array, used for expression
119  * templates
120  * iterator is a STL-style iterator
121  * const_iterator is an STL-style const iterator
122  */
123 
124  typedef P_numtype T_numtype;
127  typedef FastArrayIterator<T_numtype, N_rank> T_iterator;
128 
129  typedef ArrayIterator<T_numtype,N_rank> iterator;
130  typedef ConstArrayIterator<T_numtype,N_rank> const_iterator;
131 
132  static const int _bz_rank = N_rank;
133 
135  // Constructors //
137 
138 
139  /*
140  * Construct an array from an array expression.
141  */
142 
143  template<typename T_expr>
144  explicit Array(_bz_ArrayExpr<T_expr> expr);
145 
146  /*
147  * Any missing length arguments will have their value taken from the
148  * last argument. For example,
149  * Array<int,3> A(32,64);
150  * will create a 32x64x64 array. This is handled by setupStorage().
151  */
152 
153  Array(GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
154  : storage_(storage)
155  {
156  length_ = 0;
157  stride_ = 0;
158  zeroOffset_ = 0;
159  }
160 
161  explicit Array(int length0,
162  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
163  : storage_(storage)
164  {
165  length_[0] = length0;
166  setupStorage(0);
167  }
168 
169  Array(int length0, int length1,
170  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
171  : storage_(storage)
172  {
173  BZPRECONDITION(N_rank >= 2);
174  TAU_TYPE_STRING(p1, "Array<T,N>::Array() [T="
175  + CT(T_numtype) + ",N=" + CT(N_rank) + "]");
176  TAU_PROFILE(p1, "void (int,int)", TAU_BLITZ);
177 
178  length_[0] = length0;
179  length_[1] = length1;
180  setupStorage(1);
181  }
182 
183  Array(int length0, int length1, int length2,
184  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
185  : storage_(storage)
186  {
187  BZPRECONDITION(N_rank >= 3);
188  length_[0] = length0;
189  length_[1] = length1;
190  length_[2] = length2;
191  setupStorage(2);
192  }
193 
194  Array(int length0, int length1, int length2, int length3,
195  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
196  : storage_(storage)
197  {
198  BZPRECONDITION(N_rank >= 4);
199  length_[0] = length0;
200  length_[1] = length1;
201  length_[2] = length2;
202  length_[3] = length3;
203  setupStorage(3);
204  }
205 
206  Array(int length0, int length1, int length2, int length3, int length4,
207  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
208  : storage_(storage)
209  {
210  BZPRECONDITION(N_rank >= 5);
211  length_[0] = length0;
212  length_[1] = length1;
213  length_[2] = length2;
214  length_[3] = length3;
215  length_[4] = length4;
216  setupStorage(4);
217  }
218 
219  Array(int length0, int length1, int length2, int length3, int length4,
220  int length5,
221  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
222  : storage_(storage)
223  {
224  BZPRECONDITION(N_rank >= 6);
225  length_[0] = length0;
226  length_[1] = length1;
227  length_[2] = length2;
228  length_[3] = length3;
229  length_[4] = length4;
230  length_[5] = length5;
231  setupStorage(5);
232  }
233 
234  Array(int length0, int length1, int length2, int length3, int length4,
235  int length5, int length6,
236  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
237  : storage_(storage)
238  {
239  BZPRECONDITION(N_rank >= 7);
240  length_[0] = length0;
241  length_[1] = length1;
242  length_[2] = length2;
243  length_[3] = length3;
244  length_[4] = length4;
245  length_[5] = length5;
246  length_[6] = length6;
247  setupStorage(6);
248  }
249 
250  Array(int length0, int length1, int length2, int length3, int length4,
251  int length5, int length6, int length7,
252  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
253  : storage_(storage)
254  {
255  BZPRECONDITION(N_rank >= 8);
256  length_[0] = length0;
257  length_[1] = length1;
258  length_[2] = length2;
259  length_[3] = length3;
260  length_[4] = length4;
261  length_[5] = length5;
262  length_[6] = length6;
263  length_[7] = length7;
264  setupStorage(7);
265  }
266 
267  Array(int length0, int length1, int length2, int length3, int length4,
268  int length5, int length6, int length7, int length8,
269  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
270  : storage_(storage)
271  {
272  BZPRECONDITION(N_rank >= 9);
273  length_[0] = length0;
274  length_[1] = length1;
275  length_[2] = length2;
276  length_[3] = length3;
277  length_[4] = length4;
278  length_[5] = length5;
279  length_[6] = length6;
280  length_[7] = length7;
281  length_[8] = length8;
282  setupStorage(8);
283  }
284 
285  Array(int length0, int length1, int length2, int length3, int length4,
286  int length5, int length6, int length7, int length8, int length9,
287  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
288  : storage_(storage)
289  {
290  BZPRECONDITION(N_rank >= 10);
291  length_[0] = length0;
292  length_[1] = length1;
293  length_[2] = length2;
294  length_[3] = length3;
295  length_[4] = length4;
296  length_[5] = length5;
297  length_[6] = length6;
298  length_[7] = length7;
299  length_[8] = length8;
300  length_[9] = length9;
301  setupStorage(9);
302  }
303 
304  Array(int length0, int length1, int length2, int length3, int length4,
305  int length5, int length6, int length7, int length8, int length9,
306  int length10,
307  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
308  : storage_(storage)
309  {
310  BZPRECONDITION(N_rank >= 11);
311  length_[0] = length0;
312  length_[1] = length1;
313  length_[2] = length2;
314  length_[3] = length3;
315  length_[4] = length4;
316  length_[5] = length5;
317  length_[6] = length6;
318  length_[7] = length7;
319  length_[8] = length8;
320  length_[9] = length9;
321  length_[10] = length10;
322  setupStorage(10);
323  }
324 
325  /*
326  * Construct an array from an existing block of memory. Ownership
327  * is not acquired (this is provided for backwards compatibility).
328  */
329  Array(T_numtype* restrict dataFirst, TinyVector<int, N_rank> shape,
330  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
331  : MemoryBlockReference<T_numtype>(product(shape), dataFirst,
333  storage_(storage)
334  {
335  BZPRECONDITION(dataFirst != 0);
336 
337  length_ = shape;
338  computeStrides();
339  data_ += zeroOffset_;
340  }
341 
342  /*
343  * Construct an array from an existing block of memory, with a
344  * given set of strides. Ownership is not acquired (i.e. the memory
345  * block will not be freed by Blitz++).
346  */
347  Array(T_numtype* restrict dataFirst, TinyVector<int, N_rank> shape,
348  TinyVector<int, N_rank> stride,
349  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
350  : MemoryBlockReference<T_numtype>(product(shape), dataFirst,
352  storage_(storage)
353  {
354  BZPRECONDITION(dataFirst != 0);
355 
356  length_ = shape;
357  stride_ = stride;
358  calculateZeroOffset();
359  data_ += zeroOffset_;
360  }
361 
362  /*
363  * Construct an array from an existing block of memory.
364  */
365  Array(T_numtype* restrict dataFirst, TinyVector<int, N_rank> shape,
366  preexistingMemoryPolicy deletionPolicy,
367  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
368  : MemoryBlockReference<T_numtype>(product(shape), dataFirst,
369  deletionPolicy),
370  storage_(storage)
371  {
372  BZPRECONDITION(dataFirst != 0);
373 
374  length_ = shape;
375  computeStrides();
376  data_ += zeroOffset_;
377 
378  if (deletionPolicy == duplicateData)
379  reference(copy());
380  }
381 
382  /*
383  * Construct an array from an existing block of memory, with a
384  * given set of strides.
385  */
386  Array(T_numtype* restrict dataFirst, TinyVector<int, N_rank> shape,
388  preexistingMemoryPolicy deletionPolicy,
389  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
390  : MemoryBlockReference<T_numtype>(product(shape), dataFirst,
391  deletionPolicy),
392  storage_(storage)
393  {
394  BZPRECONDITION(dataFirst != 0);
395 
396  length_ = shape;
397  stride_ = stride;
398  calculateZeroOffset();
399  data_ += zeroOffset_;
400 
401  if (deletionPolicy == duplicateData)
402  reference(copy());
403  }
404 
405  /*
406  * This constructor takes an extent (length) vector and storage format.
407  */
408 
410  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
411  : storage_(storage)
412  {
413  length_ = extent;
414  setupStorage(N_rank - 1);
415  }
416 
417  /*
418  * This construct takes a vector of bases (lbounds) and a vector of
419  * extents.
420  */
421 
422  Array(const TinyVector<int, N_rank>& lbounds,
423  const TinyVector<int, N_rank>& extent,
424  const GeneralArrayStorage<N_rank>& storage
425  = GeneralArrayStorage<N_rank>());
426 
427  /*
428  * These constructors allow arbitrary bases (starting indices) to be set.
429  * e.g. Array<int,2> A(Range(10,20), Range(20,30))
430  * will create an 11x11 array whose indices are 10..20 and 20..30
431  */
433  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
434  : storage_(storage)
435  {
436  BZPRECONDITION(r0.isAscendingContiguous());
437 
438  length_[0] = r0.length();
439  storage_.setBase(0, r0.first());
440  setupStorage(0);
441  }
442 
443  Array(Range r0, Range r1,
444  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
445  : storage_(storage)
446  {
447  BZPRECONDITION(r0.isAscendingContiguous() &&
448  r1.isAscendingContiguous());
449 
450  length_[0] = r0.length();
451  storage_.setBase(0, r0.first());
452  length_[1] = r1.length();
453  storage_.setBase(1, r1.first());
454 
455  setupStorage(1);
456  }
457 
458  Array(Range r0, Range r1, Range r2,
459  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
460  : storage_(storage)
461  {
462  BZPRECONDITION(r0.isAscendingContiguous() &&
464 
465  length_[0] = r0.length();
466  storage_.setBase(0, r0.first());
467  length_[1] = r1.length();
468  storage_.setBase(1, r1.first());
469  length_[2] = r2.length();
470  storage_.setBase(2, r2.first());
471 
472  setupStorage(2);
473  }
474 
475  Array(Range r0, Range r1, Range r2, Range r3,
476  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
477  : storage_(storage)
478  {
479  BZPRECONDITION(r0.isAscendingContiguous() &&
481  && r3.isAscendingContiguous());
482 
483  length_[0] = r0.length();
484  storage_.setBase(0, r0.first());
485  length_[1] = r1.length();
486  storage_.setBase(1, r1.first());
487  length_[2] = r2.length();
488  storage_.setBase(2, r2.first());
489  length_[3] = r3.length();
490  storage_.setBase(3, r3.first());
491 
492  setupStorage(3);
493  }
494 
495  Array(Range r0, Range r1, Range r2, Range r3, Range r4,
496  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
497  : storage_(storage)
498  {
499  BZPRECONDITION(r0.isAscendingContiguous() &&
502 
503  length_[0] = r0.length();
504  storage_.setBase(0, r0.first());
505  length_[1] = r1.length();
506  storage_.setBase(1, r1.first());
507  length_[2] = r2.length();
508  storage_.setBase(2, r2.first());
509  length_[3] = r3.length();
510  storage_.setBase(3, r3.first());
511  length_[4] = r4.length();
512  storage_.setBase(4, r4.first());
513 
514  setupStorage(4);
515  }
516 
517  Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
518  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
519  : storage_(storage)
520  {
521  BZPRECONDITION(r0.isAscendingContiguous() &&
524  && r5.isAscendingContiguous());
525 
526  length_[0] = r0.length();
527  storage_.setBase(0, r0.first());
528  length_[1] = r1.length();
529  storage_.setBase(1, r1.first());
530  length_[2] = r2.length();
531  storage_.setBase(2, r2.first());
532  length_[3] = r3.length();
533  storage_.setBase(3, r3.first());
534  length_[4] = r4.length();
535  storage_.setBase(4, r4.first());
536  length_[5] = r5.length();
537  storage_.setBase(5, r5.first());
538 
539  setupStorage(5);
540  }
541 
542  Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
543  Range r6,
544  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
545  : storage_(storage)
546  {
547  BZPRECONDITION(r0.isAscendingContiguous() &&
551 
552  length_[0] = r0.length();
553  storage_.setBase(0, r0.first());
554  length_[1] = r1.length();
555  storage_.setBase(1, r1.first());
556  length_[2] = r2.length();
557  storage_.setBase(2, r2.first());
558  length_[3] = r3.length();
559  storage_.setBase(3, r3.first());
560  length_[4] = r4.length();
561  storage_.setBase(4, r4.first());
562  length_[5] = r5.length();
563  storage_.setBase(5, r5.first());
564  length_[6] = r6.length();
565  storage_.setBase(6, r6.first());
566 
567  setupStorage(6);
568  }
569 
570  Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
571  Range r6, Range r7,
572  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
573  : storage_(storage)
574  {
575  BZPRECONDITION(r0.isAscendingContiguous() &&
579  && r7.isAscendingContiguous());
580 
581  length_[0] = r0.length();
582  storage_.setBase(0, r0.first());
583  length_[1] = r1.length();
584  storage_.setBase(1, r1.first());
585  length_[2] = r2.length();
586  storage_.setBase(2, r2.first());
587  length_[3] = r3.length();
588  storage_.setBase(3, r3.first());
589  length_[4] = r4.length();
590  storage_.setBase(4, r4.first());
591  length_[5] = r5.length();
592  storage_.setBase(5, r5.first());
593  length_[6] = r6.length();
594  storage_.setBase(6, r6.first());
595  length_[7] = r7.length();
596  storage_.setBase(7, r7.first());
597 
598  setupStorage(7);
599  }
600 
601  Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
602  Range r6, Range r7, Range r8,
603  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
604  : storage_(storage)
605  {
606  BZPRECONDITION(r0.isAscendingContiguous() &&
611 
612  length_[0] = r0.length();
613  storage_.setBase(0, r0.first());
614  length_[1] = r1.length();
615  storage_.setBase(1, r1.first());
616  length_[2] = r2.length();
617  storage_.setBase(2, r2.first());
618  length_[3] = r3.length();
619  storage_.setBase(3, r3.first());
620  length_[4] = r4.length();
621  storage_.setBase(4, r4.first());
622  length_[5] = r5.length();
623  storage_.setBase(5, r5.first());
624  length_[6] = r6.length();
625  storage_.setBase(6, r6.first());
626  length_[7] = r7.length();
627  storage_.setBase(7, r7.first());
628  length_[8] = r8.length();
629  storage_.setBase(8, r8.first());
630 
631  setupStorage(8);
632  }
633 
634  Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
635  Range r6, Range r7, Range r8, Range r9,
636  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
637  : storage_(storage)
638  {
639  BZPRECONDITION(r0.isAscendingContiguous() &&
644  && r9.isAscendingContiguous());
645 
646  length_[0] = r0.length();
647  storage_.setBase(0, r0.first());
648  length_[1] = r1.length();
649  storage_.setBase(1, r1.first());
650  length_[2] = r2.length();
651  storage_.setBase(2, r2.first());
652  length_[3] = r3.length();
653  storage_.setBase(3, r3.first());
654  length_[4] = r4.length();
655  storage_.setBase(4, r4.first());
656  length_[5] = r5.length();
657  storage_.setBase(5, r5.first());
658  length_[6] = r6.length();
659  storage_.setBase(6, r6.first());
660  length_[7] = r7.length();
661  storage_.setBase(7, r7.first());
662  length_[8] = r8.length();
663  storage_.setBase(8, r8.first());
664  length_[9] = r9.length();
665  storage_.setBase(9, r9.first());
666 
667  setupStorage(9);
668  }
669 
670  Array(Range r0, Range r1, Range r2, Range r3, Range r4, Range r5,
671  Range r6, Range r7, Range r8, Range r9, Range r10,
672  GeneralArrayStorage<N_rank> storage = GeneralArrayStorage<N_rank>())
673  : storage_(storage)
674  {
675  BZPRECONDITION(r0.isAscendingContiguous() &&
681 
682  length_[0] = r0.length();
683  storage_.setBase(0, r0.first());
684  length_[1] = r1.length();
685  storage_.setBase(1, r1.first());
686  length_[2] = r2.length();
687  storage_.setBase(2, r2.first());
688  length_[3] = r3.length();
689  storage_.setBase(3, r3.first());
690  length_[4] = r4.length();
691  storage_.setBase(4, r4.first());
692  length_[5] = r5.length();
693  storage_.setBase(5, r5.first());
694  length_[6] = r6.length();
695  storage_.setBase(6, r6.first());
696  length_[7] = r7.length();
697  storage_.setBase(7, r7.first());
698  length_[8] = r8.length();
699  storage_.setBase(8, r8.first());
700  length_[9] = r9.length();
701  storage_.setBase(9, r9.first());
702  length_[10] = r10.length();
703  storage_.setBase(10, r10.first());
704 
705  setupStorage(10);
706  }
707 
708  /*
709  * Create a reference of another array
710  */
712 #ifdef BZ_NEW_EXPRESSION_TEMPLATES
715 #else
717 #endif
718  {
719  // NEEDS_WORK: this const_cast is a tad ugly.
720  reference(const_cast<T_array&>(array));
721  }
722 
723  /*
724  * These constructors are used for creating interlaced arrays (see
725  * <blitz/arrayshape.h>
726  */
727  Array(const TinyVector<int,N_rank-1>& shape,
728  int lastExtent, const GeneralArrayStorage<N_rank>& storage);
729  //Array(const TinyVector<Range,N_rank-1>& shape,
730  // int lastExtent, const GeneralArrayStorage<N_rank>& storage);
731 
732  /*
733  * These constructors make the array a view of a subportion of another
734  * array. If there fewer than N_rank Range arguments provided, no
735  * slicing is performed in the unspecified ranks.
736  * e.g. Array<int,3> A(20,20,20);
737  * Array<int,3> B(A, Range(5,15));
738  * is equivalent to:
739  * Array<int,3> B(A, Range(5,15), Range::all(), Range::all());
740  */
742  {
743  constructSubarray(array, r0);
744  }
745 
747  {
748  constructSubarray(array, r0, r1);
749  }
750 
752  {
753  constructSubarray(array, r0, r1, r2);
754  }
755 
757  Range r3)
758  {
759  constructSubarray(array, r0, r1, r2, r3);
760  }
761 
763  Range r3, Range r4)
764  {
765  constructSubarray(array, r0, r1, r2, r3, r4);
766  }
767 
769  Range r3, Range r4, Range r5)
770  {
771  constructSubarray(array, r0, r1, r2, r3, r4, r5);
772  }
773 
775  Range r3, Range r4, Range r5, Range r6)
776  {
777  constructSubarray(array, r0, r1, r2, r3, r4, r5, r6);
778  }
779 
781  Range r3, Range r4, Range r5, Range r6, Range r7)
782  {
783  constructSubarray(array, r0, r1, r2, r3, r4, r5, r6, r7);
784  }
785 
787  Range r3, Range r4, Range r5, Range r6, Range r7, Range r8)
788  {
789  constructSubarray(array, r0, r1, r2, r3, r4, r5, r6, r7, r8);
790  }
791 
793  Range r3, Range r4, Range r5, Range r6, Range r7, Range r8, Range r9)
794  {
795  constructSubarray(array, r0, r1, r2, r3, r4, r5, r6, r7, r8, r9);
796  }
797 
799  Range r3, Range r4, Range r5, Range r6, Range r7, Range r8, Range r9,
800  Range r10)
801  {
802  constructSubarray(array, r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10);
803  }
804 
806  const RectDomain<N_rank>& subdomain)
807  {
808  constructSubarray(array, subdomain);
809  }
810 
811  /* Constructor added by Julian Cummings */
813  const StridedDomain<N_rank>& subdomain)
814  {
815  constructSubarray(array, subdomain);
816  }
817 
818  /*
819  * This constructor is invoked by the operator()'s which take
820  * a combination of integer and Range arguments. It's not intended
821  * for end-user use.
822  */
823  template<int N_rank2, typename R0, typename R1, typename R2, typename R3, typename R4,
824  typename R5, typename R6, typename R7, typename R8, typename R9, typename R10>
825  Array(Array<T_numtype,N_rank2>& array, R0 r0, R1 r1, R2 r2,
826  R3 r3, R4 r4, R5 r5, R6 r6, R7 r7, R8 r8, R9 r9, R10 r10)
827  {
828  constructSlice(array, r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10);
829  }
830 
832  // Member functions
834 
835  const TinyVector<int, N_rank>& base() const
836  { return storage_.base(); }
837 
838  int base(int rank) const
839  { return storage_.base(rank); }
840 
841  iterator begin()
842  { return iterator(*this); }
843 
844  const_iterator begin() const
845  { return const_iterator(*this); }
846 
847  T_iterator beginFast() const
848  { return T_iterator(*this); }
849 
850  // Deprecated: now extractComponent(...)
851  template<typename P_numtype2>
852  Array<P_numtype2,N_rank> chopComponent(P_numtype2 a, int compNum,
853  int numComponents) const
854  { return extractComponent(a, compNum, numComponents); }
855 
856  int cols() const
857  { return length_[1]; }
858 
859  int columns() const
860  { return length_[1]; }
861 
862  T_array copy() const;
863 
864  // data_ always refers to the point (0,0,...,0) which may
865  // not be in the array if the base is not zero in each rank.
866  // These data() routines return a pointer to the first
867  // element in the array (but note that it may not be
868  // stored first in memory if some ranks are stored descending).
869 
870  int dataOffset() const
871  {
872  return dot(storage_.base(), stride_);
873  }
874 
875  const T_numtype* restrict data() const
876  { return data_ + dataOffset(); }
877 
878  T_numtype* restrict data()
879  { return data_ + dataOffset(); }
880 
881  // These dataZero() routines refer to the point (0,0,...,0)
882  // which may not be in the array if the bases are nonzero.
883 
884  const T_numtype* restrict dataZero() const
885  { return data_; }
886 
887  T_numtype* restrict dataZero()
888  { return data_; }
889 
890  // These dataFirst() routines refer to the element in the
891  // array which falls first in memory.
892 
893  int dataFirstOffset() const
894  {
895  int pos = 0;
896 
897  // Used to use tinyvector expressions:
898  // return data_ + dot(storage_.base()
899  // + (1 - storage_.ascendingFlag()) * (length_ - 1), stride_);
900 
901  for (int i=0; i < N_rank; ++i)
902  pos += (storage_.base(i) + (1-storage_.isRankStoredAscending(i)) *
903  (length_(i)-1)) * stride_(i);
904 
905  return pos;
906  }
907 
908  const T_numtype* restrict dataFirst() const
909  {
910  return data_ + dataFirstOffset();
911  }
912 
913  T_numtype* restrict dataFirst()
914  {
915  return data_ + dataFirstOffset();
916  }
917 
918  int depth() const
919  { return length_[2]; }
920 
921  int dimensions() const
922  { return N_rank; }
923 
924  RectDomain<N_rank> domain() const
925  {
926  return RectDomain<N_rank>(lbound(), ubound());
927  }
928 
929  void dumpStructureInformation(ostream& os = cout) const;
930 
931  iterator end()
932  {
933  return iterator();
934  }
935 
936  const_iterator end() const
937  {
938  return const_iterator();
939  }
940 
941  int extent(int rank) const
942  { return length_[rank]; }
943 
944  const TinyVector<int,N_rank>& extent() const
945  { return length_; }
946 
947  template<typename P_numtype2>
948  Array<P_numtype2,N_rank> extractComponent(P_numtype2, int compNum,
949  int numComponents) const;
950 
951  void free()
952  {
953  changeToNullBlock();
954  length_ = 0;
955  }
956 
957  bool isMajorRank(int rank) const { return storage_.ordering(rank) == 0; }
958  bool isMinorRank(int rank) const { return storage_.ordering(rank) != 0; }
959  bool isRankStoredAscending(int rank) const {
960  return storage_.isRankStoredAscending(rank);
961  }
962 
963  bool isStorageContiguous() const;
964 
965  int lbound(int rank) const { return base(rank); }
966  TinyVector<int,N_rank> lbound() const { return base(); }
967 
968  int length(int rank) const { return length_[rank]; }
969  const TinyVector<int, N_rank>& length() const { return length_; }
970 
971  void makeUnique();
972 
973  int numElements() const { return product(length_); }
974 
975  // NEEDS_WORK -- Expose the numReferences() method
976  // MemoryBlockReference<T_numtype>::numReferences;
977 
978  // The storage_.ordering_ array is a list of dimensions from
979  // the most minor (stride 1) to major dimension. Generally,
980  // ordering(0) will return the dimension which has the smallest
981  // stride, and ordering(N_rank-1) will return the dimension with
982  // the largest stride.
983  int ordering(int storageRankIndex) const
984  { return storage_.ordering(storageRankIndex); }
985 
986  const TinyVector<int, N_rank>& ordering() const
987  { return storage_.ordering(); }
988 
989  void transposeSelf(int r0, int r1, int r2=0,
990  int r3=0, int r4=0, int r5=0, int r6=0, int r7=0, int r8=0, int
991  r9=0, int r10=0);
992  T_array transpose(int r0, int r1, int r2=0,
993  int r3=0, int r4=0, int r5=0, int r6=0, int r7=0, int r8=0, int
994  r9=0, int r10=0);
995 
996  int rank() const
997  { return N_rank; }
998 
999  void reference(const T_array&);
1000 
1001  // Added by Derrick Bass
1002  T_array reindex(const TinyVector<int,N_rank>&);
1003  void reindexSelf(const
1005 
1006  void resize(int extent);
1007  void resize(int extent1, int extent2);
1008  void resize(int extent1, int extent2,
1009  int extent3);
1010  void resize(int extent1, int extent2,
1011  int extent3, int extent4);
1012  void resize(int extent1, int extent2,
1013  int extent3, int extent4, int extent5);
1014  void resize(int extent1, int extent2,
1015  int extent3, int extent4, int extent5,
1016  int extent6);
1017  void resize(int extent1, int extent2,
1018  int extent3, int extent4, int extent5,
1019  int extent6, int extent7);
1020  void resize(int extent1, int extent2,
1021  int extent3, int extent4, int extent5,
1022  int extent6, int extent7, int extent8);
1023  void resize(int extent1, int extent2,
1024  int extent3, int extent4, int extent5,
1025  int extent6, int extent7, int extent8,
1026  int extent9);
1027  void resize(int extent1, int extent2,
1028  int extent3, int extent4, int extent5,
1029  int extent6, int extent7, int extent8,
1030  int extent9, int extent10);
1031  void resize(int extent1, int extent2,
1032  int extent3, int extent4, int extent5,
1033  int extent6, int extent7, int extent8,
1034  int extent9, int extent10,
1035  int extent11);
1036 
1037 
1038  void resize(Range r1);
1039  void resize(Range r1, Range r2);
1040  void resize(Range r1, Range r2, Range r3);
1041  void resize(Range r1, Range r2, Range r3,
1042  Range r4);
1043  void resize(Range r1, Range r2, Range r3,
1044  Range r4, Range r5);
1045  void resize(Range r1, Range r2, Range r3,
1046  Range r4, Range r5, Range r6);
1047  void resize(Range r1, Range r2, Range r3,
1048  Range r4, Range r5, Range r6,
1049  Range r7);
1050  void resize(Range r1, Range r2, Range r3,
1051  Range r4, Range r5, Range r6,
1052  Range r7, Range r8);
1053  void resize(Range r1, Range r2, Range r3,
1054  Range r4, Range r5, Range r6,
1055  Range r7, Range r8, Range r9);
1056  void resize(Range r1, Range r2, Range r3,
1057  Range r4, Range r5, Range r6,
1058  Range r7, Range r8, Range r9,
1059  Range r10);
1060  void resize(Range r1, Range r2, Range r3,
1061  Range r4, Range r5, Range r6,
1062  Range r7, Range r8, Range r9,
1063  Range r10, Range r11);
1064 
1065  void resize(const TinyVector<int,N_rank>&);
1066 
1067 
1068  void resizeAndPreserve(const TinyVector<int,
1069  N_rank>&);
1070  void resizeAndPreserve(int extent);
1071  void resizeAndPreserve(int extent1,
1072  int extent2);
1073  void resizeAndPreserve(int extent1,
1074  int extent2, int extent3);
1075  void resizeAndPreserve(int extent1,
1076  int extent2, int extent3, int extent4);
1077  void resizeAndPreserve(int extent1,
1078  int extent2, int extent3, int extent4,
1079  int extent5);
1080  void resizeAndPreserve(int extent1,
1081  int extent2, int extent3, int extent4,
1082  int extent5, int extent6);
1083  void resizeAndPreserve(int extent1,
1084  int extent2, int extent3, int extent4,
1085  int extent5, int extent6, int extent7);
1086  void resizeAndPreserve(int extent1,
1087  int extent2, int extent3, int extent4,
1088  int extent5, int extent6, int extent7,
1089  int extent8);
1090  void resizeAndPreserve(int extent1,
1091  int extent2, int extent3, int extent4,
1092  int extent5, int extent6, int extent7,
1093  int extent8, int extent9);
1094  void resizeAndPreserve(int extent1,
1095  int extent2, int extent3, int extent4,
1096  int extent5, int extent6, int extent7,
1097  int extent8, int extent9,
1098  int extent10);
1099  void resizeAndPreserve(int extent1,
1100  int extent2, int extent3, int extent4,
1101  int extent5, int extent6, int extent7,
1102  int extent8, int extent9, int extent10,
1103  int extent11);
1104 
1105  // NEEDS_WORK -- resizeAndPreserve(Range,...)
1106  // NEEDS_WORK -- resizeAndPreserve(const Domain<N_rank>&);
1107 
1108  T_array reverse(int rank);
1109  void reverseSelf(int rank);
1110 
1111  int rows() const
1112  { return length_[0]; }
1113 
1114  void setStorage(GeneralArrayStorage<N_rank>);
1115 
1116  void slice(int rank, Range r);
1117 
1118  const TinyVector<int, N_rank>& shape() const
1119  { return length_; }
1120 
1121  int size() const
1122  { return numElements(); }
1123 
1124  const TinyVector<int, N_rank>& stride() const
1125  { return stride_; }
1126 
1127  int stride(int rank) const
1128  { return stride_[rank]; }
1129 
1130  int ubound(int rank) const
1131  { return base(rank) + length_(rank) - 1; }
1132 
1134  {
1136  for (int i=0; i < N_rank; ++i)
1137  ub(i) = base(i) + extent(i) - 1;
1138  // WAS: ub = base() + extent() - 1;
1139  return ub;
1140  }
1141 
1142  int zeroOffset() const
1143  { return zeroOffset_; }
1144 
1146  // Debugging routines
1148 
1149  bool isInRangeForDim(int i, int d) const {
1150  return i >= base(d) && (i - base(d)) < length_[d];
1151  }
1152 
1153  bool isInRange(int i0) const {
1154  return i0 >= base(0) && (i0 - base(0)) < length_[0];
1155  }
1156 
1157  bool isInRange(int i0, int i1) const {
1158  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1159  && i1 >= base(1) && (i1 - base(1)) < length_[1];
1160  }
1161 
1162  bool isInRange(int i0, int i1, int i2) const {
1163  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1164  && i1 >= base(1) && (i1 - base(1)) < length_[1]
1165  && i2 >= base(2) && (i2 - base(2)) < length_[2];
1166  }
1167 
1168  bool isInRange(int i0, int i1, int i2, int i3) const {
1169  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1170  && i1 >= base(1) && (i1 - base(1)) < length_[1]
1171  && i2 >= base(2) && (i2 - base(2)) < length_[2]
1172  && i3 >= base(3) && (i3 - base(3)) < length_[3];
1173  }
1174 
1175  bool isInRange(int i0, int i1, int i2, int i3, int i4) const {
1176  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1177  && i1 >= base(1) && (i1 - base(1)) < length_[1]
1178  && i2 >= base(2) && (i2 - base(2)) < length_[2]
1179  && i3 >= base(3) && (i3 - base(3)) < length_[3]
1180  && i4 >= base(4) && (i4 - base(4)) < length_[4];
1181  }
1182 
1183  bool isInRange(int i0, int i1, int i2, int i3, int i4, int i5) const {
1184  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1185  && i1 >= base(1) && (i1 - base(1)) < length_[1]
1186  && i2 >= base(2) && (i2 - base(2)) < length_[2]
1187  && i3 >= base(3) && (i3 - base(3)) < length_[3]
1188  && i4 >= base(4) && (i4 - base(4)) < length_[4]
1189  && i5 >= base(5) && (i5 - base(5)) < length_[5];
1190  }
1191 
1192  bool isInRange(int i0, int i1, int i2, int i3, int i4, int i5, int i6) const {
1193  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1194  && i1 >= base(1) && (i1 - base(1)) < length_[1]
1195  && i2 >= base(2) && (i2 - base(2)) < length_[2]
1196  && i3 >= base(3) && (i3 - base(3)) < length_[3]
1197  && i4 >= base(4) && (i4 - base(4)) < length_[4]
1198  && i5 >= base(5) && (i5 - base(5)) < length_[5]
1199  && i6 >= base(6) && (i6 - base(6)) < length_[6];
1200  }
1201 
1202  bool isInRange(int i0, int i1, int i2, int i3, int i4,
1203  int i5, int i6, int i7) const {
1204  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1205  && i1 >= base(1) && (i1 - base(1)) < length_[1]
1206  && i2 >= base(2) && (i2 - base(2)) < length_[2]
1207  && i3 >= base(3) && (i3 - base(3)) < length_[3]
1208  && i4 >= base(4) && (i4 - base(4)) < length_[4]
1209  && i5 >= base(5) && (i5 - base(5)) < length_[5]
1210  && i6 >= base(6) && (i6 - base(6)) < length_[6]
1211  && i7 >= base(7) && (i7 - base(7)) < length_[7];
1212  }
1213 
1214  bool isInRange(int i0, int i1, int i2, int i3, int i4,
1215  int i5, int i6, int i7, int i8) const {
1216  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1217  && i1 >= base(1) && (i1 - base(1)) < length_[1]
1218  && i2 >= base(2) && (i2 - base(2)) < length_[2]
1219  && i3 >= base(3) && (i3 - base(3)) < length_[3]
1220  && i4 >= base(4) && (i4 - base(4)) < length_[4]
1221  && i5 >= base(5) && (i5 - base(5)) < length_[5]
1222  && i6 >= base(6) && (i6 - base(6)) < length_[6]
1223  && i7 >= base(7) && (i7 - base(7)) < length_[7]
1224  && i8 >= base(8) && (i8 - base(8)) < length_[8];
1225  }
1226 
1227  bool isInRange(int i0, int i1, int i2, int i3, int i4,
1228  int i5, int i6, int i7, int i8, int i9) const {
1229  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1230  && i1 >= base(1) && (i1 - base(1)) < length_[1]
1231  && i2 >= base(2) && (i2 - base(2)) < length_[2]
1232  && i3 >= base(3) && (i3 - base(3)) < length_[3]
1233  && i4 >= base(4) && (i4 - base(4)) < length_[4]
1234  && i5 >= base(5) && (i5 - base(5)) < length_[5]
1235  && i6 >= base(6) && (i6 - base(6)) < length_[6]
1236  && i7 >= base(7) && (i7 - base(7)) < length_[7]
1237  && i8 >= base(8) && (i8 - base(8)) < length_[8]
1238  && i9 >= base(9) && (i9 - base(9)) < length_[9];
1239  }
1240 
1241  bool isInRange(int i0, int i1, int i2, int i3, int i4,
1242  int i5, int i6, int i7, int i8, int i9, int i10) const {
1243  return i0 >= base(0) && (i0 - base(0)) < length_[0]
1244  && i1 >= base(1) && (i1 - base(1)) < length_[1]
1245  && i2 >= base(2) && (i2 - base(2)) < length_[2]
1246  && i3 >= base(3) && (i3 - base(3)) < length_[3]
1247  && i4 >= base(4) && (i4 - base(4)) < length_[4]
1248  && i5 >= base(5) && (i5 - base(5)) < length_[5]
1249  && i6 >= base(6) && (i6 - base(6)) < length_[6]
1250  && i7 >= base(7) && (i7 - base(7)) < length_[7]
1251  && i8 >= base(8) && (i8 - base(8)) < length_[8]
1252  && i9 >= base(9) && (i9 - base(9)) < length_[9]
1253  && i10 >= base(10) && (i10 - base(10)) < length_[10];
1254  }
1255 
1256  bool isInRange(const T_index& index) const {
1257  for (int i=0; i < N_rank; ++i)
1258  if (index[i] < base(i) || (index[i] - base(i)) >= length_[i])
1259  return false;
1260 
1261  return true;
1262  }
1263 
1264  bool assertInRange(const T_index& BZ_DEBUG_PARAM(index)) const {
1265  BZPRECHECK(isInRange(index), "Array index out of range: " << index
1266  << endl << "Lower bounds: " << storage_.base() << endl
1267  << "Length: " << length_ << endl);
1268  return true;
1269  }
1270 
1271  bool assertInRange(int BZ_DEBUG_PARAM(i0)) const {
1272  BZPRECHECK(isInRange(i0), "Array index out of range: " << i0
1273  << endl << "Lower bounds: " << storage_.base() << endl
1274  << "Length: " << length_ << endl);
1275  return true;
1276  }
1277 
1278  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1)) const {
1279  BZPRECHECK(isInRange(i0,i1), "Array index out of range: ("
1280  << i0 << ", " << i1 << ")"
1281  << endl << "Lower bounds: " << storage_.base() << endl
1282  << "Length: " << length_ << endl);
1283  return true;
1284  }
1285 
1286  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1),
1287  int BZ_DEBUG_PARAM(i2)) const
1288  {
1289  BZPRECHECK(isInRange(i0,i1,i2), "Array index out of range: ("
1290  << i0 << ", " << i1 << ", " << i2 << ")"
1291  << endl << "Lower bounds: " << storage_.base() << endl
1292  << "Length: " << length_ << endl);
1293  return true;
1294  }
1295 
1296  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1),
1297  int BZ_DEBUG_PARAM(i2), int BZ_DEBUG_PARAM(i3)) const
1298  {
1299  BZPRECHECK(isInRange(i0,i1,i2,i3), "Array index out of range: ("
1300  << i0 << ", " << i1 << ", " << i2 << ", " << i3 << ")"
1301  << endl << "Lower bounds: " << storage_.base() << endl
1302  << "Length: " << length_ << endl);
1303  return true;
1304  }
1305 
1306  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1),
1307  int BZ_DEBUG_PARAM(i2), int BZ_DEBUG_PARAM(i3),
1308  int BZ_DEBUG_PARAM(i4)) const
1309  {
1310  BZPRECHECK(isInRange(i0,i1,i2,i3,i4), "Array index out of range: ("
1311  << i0 << ", " << i1 << ", " << i2 << ", " << i3
1312  << ", " << i4 << ")"
1313  << endl << "Lower bounds: " << storage_.base() << endl
1314  << "Length: " << length_ << endl);
1315  return true;
1316  }
1317 
1318  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1),
1319  int BZ_DEBUG_PARAM(i2), int BZ_DEBUG_PARAM(i3), int BZ_DEBUG_PARAM(i4),
1320  int BZ_DEBUG_PARAM(i5)) const
1321  {
1322  BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5), "Array index out of range: ("
1323  << i0 << ", " << i1 << ", " << i2 << ", " << i3
1324  << ", " << i4 << ", " << i5 << ")"
1325  << endl << "Lower bounds: " << storage_.base() << endl
1326  << "Length: " << length_ << endl);
1327  return true;
1328  }
1329 
1330  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1),
1331  int BZ_DEBUG_PARAM(i2), int BZ_DEBUG_PARAM(i3), int BZ_DEBUG_PARAM(i4),
1332  int BZ_DEBUG_PARAM(i5), int BZ_DEBUG_PARAM(i6)) const
1333  {
1334  BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6),
1335  "Array index out of range: ("
1336  << i0 << ", " << i1 << ", " << i2 << ", " << i3
1337  << ", " << i4 << ", " << i5 << ", " << i6 << ")"
1338  << endl << "Lower bounds: " << storage_.base() << endl
1339  << "Length: " << length_ << endl);
1340  return true;
1341  }
1342 
1343  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1),
1344  int BZ_DEBUG_PARAM(i2), int BZ_DEBUG_PARAM(i3), int BZ_DEBUG_PARAM(i4),
1345  int BZ_DEBUG_PARAM(i5), int BZ_DEBUG_PARAM(i6),
1346  int BZ_DEBUG_PARAM(i7)) const
1347  {
1348  BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6,i7),
1349  "Array index out of range: ("
1350  << i0 << ", " << i1 << ", " << i2 << ", " << i3
1351  << ", " << i4 << ", " << i5 << ", " << i6 << ", " << i7 << ")"
1352  << endl << "Lower bounds: " << storage_.base() << endl
1353  << "Length: " << length_ << endl);
1354  return true;
1355  }
1356 
1357  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1),
1358  int BZ_DEBUG_PARAM(i2), int BZ_DEBUG_PARAM(i3), int BZ_DEBUG_PARAM(i4),
1359  int BZ_DEBUG_PARAM(i5), int BZ_DEBUG_PARAM(i6), int BZ_DEBUG_PARAM(i7),
1360  int BZ_DEBUG_PARAM(i8)) const
1361  {
1362  BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6,i7,i8),
1363  "Array index out of range: ("
1364  << i0 << ", " << i1 << ", " << i2 << ", " << i3
1365  << ", " << i4 << ", " << i5 << ", " << i6 << ", " << i7
1366  << ", " << i8 << ")"
1367  << endl << "Lower bounds: " << storage_.base() << endl
1368  << "Length: " << length_ << endl);
1369  return true;
1370  }
1371 
1372  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1),
1373  int BZ_DEBUG_PARAM(i2), int BZ_DEBUG_PARAM(i3), int BZ_DEBUG_PARAM(i4),
1374  int BZ_DEBUG_PARAM(i5), int BZ_DEBUG_PARAM(i6), int BZ_DEBUG_PARAM(i7),
1375  int BZ_DEBUG_PARAM(i8), int BZ_DEBUG_PARAM(i9)) const
1376  {
1377  BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6,i7,i8,i9),
1378  "Array index out of range: ("
1379  << i0 << ", " << i1 << ", " << i2 << ", " << i3
1380  << ", " << i4 << ", " << i5 << ", " << i6 << ", " << i7
1381  << ", " << i8 << ", " << i9 << ")"
1382  << endl << "Lower bounds: " << storage_.base() << endl
1383  << "Length: " << length_ << endl);
1384  return true;
1385  }
1386 
1387  bool assertInRange(int BZ_DEBUG_PARAM(i0), int BZ_DEBUG_PARAM(i1),
1388  int BZ_DEBUG_PARAM(i2), int BZ_DEBUG_PARAM(i3), int BZ_DEBUG_PARAM(i4),
1389  int BZ_DEBUG_PARAM(i5), int BZ_DEBUG_PARAM(i6), int BZ_DEBUG_PARAM(i7),
1390  int BZ_DEBUG_PARAM(i8), int BZ_DEBUG_PARAM(i9),
1391  int BZ_DEBUG_PARAM(i10)) const
1392  {
1393  BZPRECHECK(isInRange(i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10),
1394  "Array index out of range: ("
1395  << i0 << ", " << i1 << ", " << i2 << ", " << i3
1396  << ", " << i4 << ", " << i5 << ", " << i6 << ", " << i7
1397  << ", " << i8 << ", " << i9 << ", " << i10 << ")"
1398  << endl << "Lower bounds: " << storage_.base() << endl
1399  << "Length: " << length_ << endl);
1400  return true;
1401  }
1402 
1404  // Subscripting operators
1406 
1407  template<int N_rank2>
1408  const T_numtype& restrict operator()(const TinyVector<int,N_rank2>& index) const
1409  {
1410  assertInRange(index);
1411  return data_[dot(index, stride_)];
1412  }
1413 
1414  template<int N_rank2>
1415  T_numtype& restrict operator()(const TinyVector<int,N_rank2>& index)
1416  {
1417  assertInRange(index);
1418  return data_[dot(index, stride_)];
1419  }
1420 
1421  const T_numtype& restrict operator()(TinyVector<int,1> index) const
1422  {
1423  assertInRange(index[0]);
1424  return data_[index[0] * stride_[0]];
1425  }
1426 
1427  T_numtype& operator()(TinyVector<int,1> index)
1428  {
1429  assertInRange(index[0]);
1430  return data_[index[0] * stride_[0]];
1431  }
1432 
1433  const T_numtype& restrict operator()(TinyVector<int,2> index) const
1434  {
1435  assertInRange(index[0], index[1]);
1436  return data_[index[0] * stride_[0] + index[1] * stride_[1]];
1437  }
1438 
1439  T_numtype& operator()(TinyVector<int,2> index)
1440  {
1441  assertInRange(index[0], index[1]);
1442  return data_[index[0] * stride_[0] + index[1] * stride_[1]];
1443  }
1444 
1445  const T_numtype& restrict operator()(TinyVector<int,3> index) const
1446  {
1447  assertInRange(index[0], index[1], index[2]);
1448  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1449  + index[2] * stride_[2]];
1450  }
1451 
1452  T_numtype& operator()(TinyVector<int,3> index)
1453  {
1454  assertInRange(index[0], index[1], index[2]);
1455  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1456  + index[2] * stride_[2]];
1457  }
1458 
1459  const T_numtype& restrict operator()(const TinyVector<int,4>& index) const
1460  {
1461  assertInRange(index[0], index[1], index[2], index[3]);
1462  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1463  + index[2] * stride_[2] + index[3] * stride_[3]];
1464  }
1465 
1466  T_numtype& operator()(const TinyVector<int,4>& index)
1467  {
1468  assertInRange(index[0], index[1], index[2], index[3]);
1469  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1470  + index[2] * stride_[2] + index[3] * stride_[3]];
1471  }
1472 
1473  const T_numtype& restrict operator()(const TinyVector<int,5>& index) const
1474  {
1475  assertInRange(index[0], index[1], index[2], index[3],
1476  index[4]);
1477  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1478  + index[2] * stride_[2] + index[3] * stride_[3]
1479  + index[4] * stride_[4]];
1480  }
1481 
1482  T_numtype& operator()(const TinyVector<int,5>& index)
1483  {
1484  assertInRange(index[0], index[1], index[2], index[3],
1485  index[4]);
1486  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1487  + index[2] * stride_[2] + index[3] * stride_[3]
1488  + index[4] * stride_[4]];
1489  }
1490 
1491  const T_numtype& restrict operator()(const TinyVector<int,6>& index) const
1492  {
1493  assertInRange(index[0], index[1], index[2], index[3],
1494  index[4], index[5]);
1495  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1496  + index[2] * stride_[2] + index[3] * stride_[3]
1497  + index[4] * stride_[4] + index[5] * stride_[5]];
1498  }
1499 
1500  T_numtype& operator()(const TinyVector<int,6>& index)
1501  {
1502  assertInRange(index[0], index[1], index[2], index[3],
1503  index[4], index[5]);
1504  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1505  + index[2] * stride_[2] + index[3] * stride_[3]
1506  + index[4] * stride_[4] + index[5] * stride_[5]];
1507  }
1508 
1509  const T_numtype& restrict operator()(const TinyVector<int,7>& index) const
1510  {
1511  assertInRange(index[0], index[1], index[2], index[3],
1512  index[4], index[5], index[6]);
1513  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1514  + index[2] * stride_[2] + index[3] * stride_[3]
1515  + index[4] * stride_[4] + index[5] * stride_[5]
1516  + index[6] * stride_[6]];
1517  }
1518 
1519  T_numtype& operator()(const TinyVector<int,7>& index)
1520  {
1521  assertInRange(index[0], index[1], index[2], index[3],
1522  index[4], index[5], index[6]);
1523  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1524  + index[2] * stride_[2] + index[3] * stride_[3]
1525  + index[4] * stride_[4] + index[5] * stride_[5]
1526  + index[6] * stride_[6]];
1527  }
1528 
1529  const T_numtype& restrict operator()(const TinyVector<int,8>& index) const
1530  {
1531  assertInRange(index[0], index[1], index[2], index[3],
1532  index[4], index[5], index[6], index[7]);
1533  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1534  + index[2] * stride_[2] + index[3] * stride_[3]
1535  + index[4] * stride_[4] + index[5] * stride_[5]
1536  + index[6] * stride_[6] + index[7] * stride_[7]];
1537  }
1538 
1539  T_numtype& operator()(const TinyVector<int,8>& index)
1540  {
1541  assertInRange(index[0], index[1], index[2], index[3],
1542  index[4], index[5], index[6], index[7]);
1543  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1544  + index[2] * stride_[2] + index[3] * stride_[3]
1545  + index[4] * stride_[4] + index[5] * stride_[5]
1546  + index[6] * stride_[6] + index[7] * stride_[7]];
1547  }
1548 
1549  const T_numtype& restrict operator()(const TinyVector<int,9>& index) const
1550  {
1551  assertInRange(index[0], index[1], index[2], index[3],
1552  index[4], index[5], index[6], index[7], index[8]);
1553  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1554  + index[2] * stride_[2] + index[3] * stride_[3]
1555  + index[4] * stride_[4] + index[5] * stride_[5]
1556  + index[6] * stride_[6] + index[7] * stride_[7]
1557  + index[8] * stride_[8]];
1558  }
1559 
1560  T_numtype& operator()(const TinyVector<int,9>& index)
1561  {
1562  assertInRange(index[0], index[1], index[2], index[3],
1563  index[4], index[5], index[6], index[7], index[8]);
1564  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1565  + index[2] * stride_[2] + index[3] * stride_[3]
1566  + index[4] * stride_[4] + index[5] * stride_[5]
1567  + index[6] * stride_[6] + index[7] * stride_[7]
1568  + index[8] * stride_[8]];
1569  }
1570 
1571  const T_numtype& restrict operator()(const TinyVector<int,10>& index) const
1572  {
1573  assertInRange(index[0], index[1], index[2], index[3],
1574  index[4], index[5], index[6], index[7], index[8], index[9]);
1575  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1576  + index[2] * stride_[2] + index[3] * stride_[3]
1577  + index[4] * stride_[4] + index[5] * stride_[5]
1578  + index[6] * stride_[6] + index[7] * stride_[7]
1579  + index[8] * stride_[8] + index[9] * stride_[9]];
1580  }
1581 
1582  T_numtype& operator()(const TinyVector<int,10>& index)
1583  {
1584  assertInRange(index[0], index[1], index[2], index[3],
1585  index[4], index[5], index[6], index[7], index[8], index[9]);
1586  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1587  + index[2] * stride_[2] + index[3] * stride_[3]
1588  + index[4] * stride_[4] + index[5] * stride_[5]
1589  + index[6] * stride_[6] + index[7] * stride_[7]
1590  + index[8] * stride_[8] + index[9] * stride_[9]];
1591  }
1592 
1593  const T_numtype& restrict operator()(const TinyVector<int,11>& index) const
1594  {
1595  assertInRange(index[0], index[1], index[2], index[3],
1596  index[4], index[5], index[6], index[7], index[8], index[9],
1597  index[10]);
1598  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1599  + index[2] * stride_[2] + index[3] * stride_[3]
1600  + index[4] * stride_[4] + index[5] * stride_[5]
1601  + index[6] * stride_[6] + index[7] * stride_[7]
1602  + index[8] * stride_[8] + index[9] * stride_[9]
1603  + index[10] * stride_[10]];
1604  }
1605 
1606  T_numtype& operator()(const TinyVector<int,11>& index)
1607  {
1608  assertInRange(index[0], index[1], index[2], index[3],
1609  index[4], index[5], index[6], index[7], index[8], index[9],
1610  index[10]);
1611  return data_[index[0] * stride_[0] + index[1] * stride_[1]
1612  + index[2] * stride_[2] + index[3] * stride_[3]
1613  + index[4] * stride_[4] + index[5] * stride_[5]
1614  + index[6] * stride_[6] + index[7] * stride_[7]
1615  + index[8] * stride_[8] + index[9] * stride_[9]
1616  + index[10] * stride_[10]];
1617  }
1618 
1619  const T_numtype& restrict operator()(int i0) const
1620  {
1621  assertInRange(i0);
1622  return data_[i0 * stride_[0]];
1623  }
1624 
1625  T_numtype& restrict operator()(int i0)
1626  {
1627  assertInRange(i0);
1628  return data_[i0 * stride_[0]];
1629  }
1630 
1631  const T_numtype& restrict operator()(int i0, int i1) const
1632  {
1633  assertInRange(i0, i1);
1634  return data_[i0 * stride_[0] + i1 * stride_[1]];
1635  }
1636 
1637  T_numtype& restrict operator()(int i0, int i1)
1638  {
1639  assertInRange(i0, i1);
1640  return data_[i0 * stride_[0] + i1 * stride_[1]];
1641  }
1642 
1643  const T_numtype& restrict operator()(int i0, int i1, int i2) const
1644  {
1645  assertInRange(i0, i1, i2);
1646  return data_[i0 * stride_[0] + i1 * stride_[1]
1647  + i2 * stride_[2]];
1648  }
1649 
1650  T_numtype& restrict operator()(int i0, int i1, int i2)
1651  {
1652  assertInRange(i0, i1, i2);
1653  return data_[i0 * stride_[0] + i1 * stride_[1]
1654  + i2 * stride_[2]];
1655  }
1656 
1657  const T_numtype& restrict operator()(int i0, int i1, int i2, int i3) const
1658  {
1659  assertInRange(i0, i1, i2, i3);
1660  return data_[i0 * stride_[0] + i1 * stride_[1]
1661  + i2 * stride_[2] + i3 * stride_[3]];
1662  }
1663 
1664  T_numtype& restrict operator()(int i0, int i1, int i2, int i3)
1665  {
1666  assertInRange(i0, i1, i2, i3);
1667  return data_[i0 * stride_[0] + i1 * stride_[1]
1668  + i2 * stride_[2] + i3 * stride_[3]];
1669  }
1670 
1671  const T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1672  int i4) const
1673  {
1674  assertInRange(i0, i1, i2, i3, i4);
1675  return data_[i0 * stride_[0] + i1 * stride_[1]
1676  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]];
1677  }
1678 
1679  T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1680  int i4)
1681  {
1682  assertInRange(i0, i1, i2, i3, i4);
1683  return data_[i0 * stride_[0] + i1 * stride_[1]
1684  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]];
1685  }
1686 
1687  const T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1688  int i4, int i5) const
1689  {
1690  assertInRange(i0, i1, i2, i3, i4, i5);
1691  return data_[i0 * stride_[0] + i1 * stride_[1]
1692  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1693  + i5 * stride_[5]];
1694  }
1695 
1696  T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1697  int i4, int i5)
1698  {
1699  assertInRange(i0, i1, i2, i3, i4, i5);
1700  return data_[i0 * stride_[0] + i1 * stride_[1]
1701  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1702  + i5 * stride_[5]];
1703  }
1704 
1705  const T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1706  int i4, int i5, int i6) const
1707  {
1708  assertInRange(i0, i1, i2, i3, i4, i5, i6);
1709  return data_[i0 * stride_[0] + i1 * stride_[1]
1710  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1711  + i5 * stride_[5] + i6 * stride_[6]];
1712  }
1713 
1714  T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1715  int i4, int i5, int i6)
1716  {
1717  assertInRange(i0, i1, i2, i3, i4, i5, i6);
1718  return data_[i0 * stride_[0] + i1 * stride_[1]
1719  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1720  + i5 * stride_[5] + i6 * stride_[6]];
1721  }
1722 
1723  const T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1724  int i4, int i5, int i6, int i7) const
1725  {
1726  assertInRange(i0, i1, i2, i3, i4, i5, i6, i7);
1727  return data_[i0 * stride_[0] + i1 * stride_[1]
1728  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1729  + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]];
1730  }
1731 
1732  T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1733  int i4, int i5, int i6, int i7)
1734  {
1735  assertInRange(i0, i1, i2, i3, i4, i5, i6, i7);
1736  return data_[i0 * stride_[0] + i1 * stride_[1]
1737  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1738  + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]];
1739  }
1740 
1741  const T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1742  int i4, int i5, int i6, int i7, int i8) const
1743  {
1744  assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8);
1745  return data_[i0 * stride_[0] + i1 * stride_[1]
1746  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1747  + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
1748  + i8 * stride_[8]];
1749  }
1750 
1751  T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1752  int i4, int i5, int i6, int i7, int i8)
1753  {
1754  assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8);
1755  return data_[i0 * stride_[0] + i1 * stride_[1]
1756  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1757  + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
1758  + i8 * stride_[8]];
1759  }
1760 
1761  const T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1762  int i4, int i5, int i6, int i7, int i8, int i9) const
1763  {
1764  assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9);
1765  return data_[i0 * stride_[0] + i1 * stride_[1]
1766  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1767  + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
1768  + i8 * stride_[8] + i9 * stride_[9]];
1769  }
1770 
1771  T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1772  int i4, int i5, int i6, int i7, int i8, int i9)
1773  {
1774  assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9);
1775  return data_[i0 * stride_[0] + i1 * stride_[1]
1776  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1777  + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
1778  + i8 * stride_[8] + i9 * stride_[9]];
1779  }
1780 
1781  const T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1782  int i4, int i5, int i6, int i7, int i8, int i9, int i10) const
1783  {
1784  assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8,
1785  i9, i10);
1786  return data_[i0 * stride_[0] + i1 * stride_[1]
1787  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1788  + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
1789  + i8 * stride_[8] + i9 * stride_[9] + i10 * stride_[10]];
1790  }
1791 
1792  T_numtype& restrict operator()(int i0, int i1, int i2, int i3,
1793  int i4, int i5, int i6, int i7, int i8, int i9, int i10)
1794  {
1795  assertInRange(i0, i1, i2, i3, i4, i5, i6, i7, i8,
1796  i9, i10);
1797  return data_[i0 * stride_[0] + i1 * stride_[1]
1798  + i2 * stride_[2] + i3 * stride_[3] + i4 * stride_[4]
1799  + i5 * stride_[5] + i6 * stride_[6] + i7 * stride_[7]
1800  + i8 * stride_[8] + i9 * stride_[9] + i10 * stride_[10]];
1801  }
1802 
1803  /*
1804  * Slicing to produce subarrays. If the number of Range arguments is
1805  * fewer than N_rank, then missing arguments are treated like Range::all().
1806  */
1807 
1808  T_array& noConst() const
1809  { return const_cast<T_array&>(*this); }
1810 
1811  T_array operator()(const RectDomain<N_rank>& subdomain) const
1812  {
1813  return T_array(noConst(), subdomain);
1814  }
1815 
1816  /* Operator added by Julian Cummings */
1817  T_array operator()(const StridedDomain<N_rank>& subdomain) const
1818  {
1819  return T_array(noConst(), subdomain);
1820  }
1821 
1822  T_array operator()(Range r0) const
1823  {
1824  return T_array(noConst(), r0);
1825  }
1826 
1827  T_array operator()(Range r0, Range r1) const
1828  {
1829  return T_array(noConst(), r0, r1);
1830  }
1831 
1832  T_array operator()(Range r0, Range r1, Range r2) const
1833  {
1834  return T_array(noConst(), r0, r1, r2);
1835  }
1836 
1837  T_array operator()(Range r0, Range r1, Range r2, Range r3) const
1838  {
1839  return T_array(noConst(), r0, r1, r2, r3);
1840  }
1841 
1842  T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4) const
1843  {
1844  return T_array(noConst(), r0, r1, r2, r3, r4);
1845  }
1846 
1847  T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
1848  Range r5) const
1849  {
1850  return T_array(noConst(), r0, r1, r2, r3, r4, r5);
1851  }
1852 
1853  T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
1854  Range r5, Range r6) const
1855  {
1856  return T_array(noConst(), r0, r1, r2, r3, r4, r5, r6);
1857  }
1858 
1859  T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
1860  Range r5, Range r6, Range r7) const
1861  {
1862  return T_array(noConst(), r0, r1, r2, r3, r4, r5, r6, r7);
1863  }
1864 
1865  T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
1866  Range r5, Range r6, Range r7, Range r8) const
1867  {
1868  return T_array(noConst(), r0, r1, r2, r3, r4, r5, r6, r7, r8);
1869  }
1870 
1871  T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
1872  Range r5, Range r6, Range r7, Range r8, Range r9) const
1873  {
1874  return T_array(noConst(), r0, r1, r2, r3, r4, r5, r6, r7, r8, r9);
1875  }
1876 
1877  T_array operator()(Range r0, Range r1, Range r2, Range r3, Range r4,
1878  Range r5, Range r6, Range r7, Range r8, Range r9, Range r10) const
1879  {
1880  return T_array(noConst(), r0, r1, r2, r3, r4, r5, r6, r7, r8, r9, r10);
1881  }
1882 
1883  // Allow any mixture of Range, int and Vector<int> objects as
1884  // operands for operator(): A(Range(3,7), 5, Range(2,4))
1885 
1886  /*
1887  * These versions of operator() allow any combination of int
1888  * and Range operands to be used. Each int operand reduces
1889  * the rank of the resulting array by one.
1890  *
1891  * e.g. Array<int,4> A(20,20,20,20);
1892  * Array<int,2> B = A(Range(5,15), 3, 5, Range(8,9));
1893  *
1894  * SliceInfo is a helper class defined in <blitz/arrayslice.h>.
1895  * It counts the number of Range vs. int arguments and does some
1896  * other helpful things.
1897  *
1898  * Once partial specialization becomes widely implemented, these
1899  * operators may be expanded to accept Vector<int> arguments
1900  * and produce ArrayPick<T,N> objects.
1901  *
1902  * This operator() is not provided with a single argument because
1903  * the appropriate cases exist above.
1904  */
1905 
1906 #ifdef BZ_HAVE_PARTIAL_ORDERING
1907 
1908  template<typename T1, typename T2>
1909  typename SliceInfo<T_numtype,T1,T2>::T_slice
1910  operator()(T1 r1, T2 r2) const
1911  {
1912  typedef typename SliceInfo<T_numtype,T1,T2>::T_slice slice;
1913  return slice(noConst(), r1, r2, nilArraySection(), nilArraySection(), nilArraySection(),
1914  nilArraySection(), nilArraySection(), nilArraySection(),
1915  nilArraySection(), nilArraySection(), nilArraySection());
1916  }
1917 
1918  template<typename T1, typename T2, typename T3>
1919  typename SliceInfo<T_numtype,T1,T2,T3>::T_slice
1920  operator()(T1 r1, T2 r2, T3 r3) const
1921  {
1922  typedef typename SliceInfo<T_numtype,T1,T2,T3>::T_slice slice;
1923  return slice(noConst(), r1, r2, r3, nilArraySection(), nilArraySection(), nilArraySection(),
1924  nilArraySection(), nilArraySection(), nilArraySection(),
1925  nilArraySection(), nilArraySection());
1926  }
1927 
1928  template<typename T1, typename T2, typename T3, typename T4>
1929  typename SliceInfo<T_numtype,T1,T2,T3,T4>::T_slice
1930  operator()(T1 r1, T2 r2, T3 r3, T4 r4) const
1931  {
1932  typedef typename SliceInfo<T_numtype,T1,T2,T3,T4>::T_slice slice;
1933  return slice(noConst(), r1, r2, r3, r4, nilArraySection(), nilArraySection(),
1934  nilArraySection(), nilArraySection(), nilArraySection(),
1935  nilArraySection(), nilArraySection());
1936  }
1937 
1938  template<typename T1, typename T2, typename T3, typename T4, typename T5>
1939  typename SliceInfo<T_numtype,T1,T2,T3,T4,T5>::T_slice
1940  operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5) const
1941  {
1942  typedef typename SliceInfo<T_numtype,T1,T2,T3,T4,T5>::T_slice slice;
1943  return slice(noConst(), r1, r2, r3, r4, r5, nilArraySection(),
1944  nilArraySection(), nilArraySection(), nilArraySection(),
1945  nilArraySection(), nilArraySection());
1946  }
1947 
1948  template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6>
1949  typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6>::T_slice
1950  operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6) const
1951  {
1952  typedef typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6>::T_slice slice;
1953  return slice(noConst(), r1, r2, r3, r4, r5, r6, nilArraySection(), nilArraySection(), nilArraySection(),
1954  nilArraySection(), nilArraySection());
1955  }
1956 
1957  template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
1958  typename T7>
1959  typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7>::T_slice
1960  operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7) const
1961  {
1962  typedef typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7>::T_slice slice;
1963  return slice(noConst(), r1, r2, r3, r4, r5, r6, r7, nilArraySection(), nilArraySection(),
1964  nilArraySection(), nilArraySection());
1965  }
1966 
1967  template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
1968  typename T7, typename T8>
1969  typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8>::T_slice
1970  operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8) const
1971  {
1972  typedef typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8>::T_slice slice;
1973  return slice(noConst(), r1, r2, r3, r4, r5, r6, r7, r8,
1974  nilArraySection(), nilArraySection(), nilArraySection());
1975  }
1976 
1977  template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
1978  typename T7, typename T8, typename T9>
1979  typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9>::T_slice
1980  operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9) const
1981  {
1982  typedef typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9>::T_slice slice;
1983  return slice(noConst(), r1, r2, r3, r4, r5, r6, r7, r8, r9, nilArraySection(), nilArraySection());
1984  }
1985 
1986  template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
1987  typename T7, typename T8, typename T9, typename T10>
1988  typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10>::T_slice
1989  operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9, T10 r10) const
1990  {
1991  typedef typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10>::T_slice slice;
1992  return slice(noConst(), r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, nilArraySection());
1993  }
1994 
1995  template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
1996  typename T7, typename T8, typename T9, typename T10, typename T11>
1997  typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
1998  operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9, T10 r10, T11 r11) const
1999  {
2000  typedef typename SliceInfo<T_numtype,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice slice;
2001  return slice(noConst(), r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11);
2002  }
2003 
2004 #endif // BZ_HAVE_PARTIAL_ORDERING
2005 
2006  /*
2007  * These versions of operator() are provided to support tensor-style
2008  * array notation, e.g.
2009  *
2010  * Array<float, 2> A, B;
2011  * firstIndex i;
2012  * secondIndex j;
2013  * thirdIndex k;
2014  * Array<float, 3> C = A(i,j) * B(j,k);
2015  */
2016 
2017  template<int N0>
2018  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0> >
2019  operator()(IndexPlaceholder<N0>) const
2020  {
2021  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0> >
2022  (noConst());
2023  }
2024 
2025  template<int N0, int N1>
2026  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1> >
2028  {
2029  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2030  N1> >(noConst());
2031  }
2032 
2033  template<int N0, int N1, int N2>
2034  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2> >
2036  IndexPlaceholder<N2>) const
2037  {
2038  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2039  N1, N2> >(noConst());
2040  }
2041 
2042  template<int N0, int N1, int N2, int N3>
2043  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3> >
2046  {
2047  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2048  N1, N2, N3> >(noConst());
2049  }
2050 
2051  template<int N0, int N1, int N2, int N3, int N4>
2052  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3, N4> >
2055  IndexPlaceholder<N4>) const
2056  {
2057  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2058  N1, N2, N3, N4> >(noConst());
2059  }
2060 
2061  template<int N0, int N1, int N2, int N3, int N4, int N5>
2062  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
2063  N4, N5> >
2066  IndexPlaceholder<N5>) const
2067  {
2068  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2069  N1, N2, N3, N4, N5> >(noConst());
2070  }
2071 
2072  template<int N0, int N1, int N2, int N3, int N4, int N5, int N6>
2073  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
2074  N4, N5, N6> >
2078  {
2079  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2080  N1, N2, N3, N4, N5, N6> >(noConst());
2081  }
2082 
2083  template<int N0, int N1, int N2, int N3, int N4, int N5, int N6,
2084  int N7>
2085  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
2086  N4, N5, N6, N7> >
2090  IndexPlaceholder<N7>) const
2091  {
2092  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2093  N1, N2, N3, N4, N5, N6, N7> >(noConst());
2094  }
2095 
2096  template<int N0, int N1, int N2, int N3, int N4, int N5, int N6,
2097  int N7, int N8>
2098  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
2099  N4, N5, N6, N7, N8> >
2103  IndexPlaceholder<N8>) const
2104  {
2105  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2106  N1, N2, N3, N4, N5, N6, N7, N8> >(noConst());
2107  }
2108 
2109  template<int N0, int N1, int N2, int N3, int N4, int N5, int N6,
2110  int N7, int N8, int N9>
2111  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
2112  N4, N5, N6, N7, N8, N9> >
2117  {
2118  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2119  N1, N2, N3, N4, N5, N6, N7, N8, N9> >(noConst());
2120  }
2121 
2122  template<int N0, int N1, int N2, int N3, int N4, int N5, int N6,
2123  int N7, int N8, int N9, int N10>
2124  _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0, N1, N2, N3,
2125  N4, N5, N6, N7, N8, N9, N10> >
2130  IndexPlaceholder<N10>) const
2131  {
2132  return _bz_ArrayExpr<ArrayIndexMapping<T_numtype, N_rank, N0,
2133  N1, N2, N3, N4, N5, N6, N7, N8, N9, N10> >(noConst());
2134  }
2135 
2137  // Support for multicomponent arrays
2139 
2140  /*
2141  * See <blitz/array/multi.h> for an explanation of the traits class
2142  * multicomponent_traits.
2143  */
2144 
2146  operator[](const unsigned component) {
2147  typedef typename multicomponent_traits<T_numtype>::T_element T_compType;
2148 
2149  return extractComponent(T_compType(),component,
2150  multicomponent_traits<T_numtype>::numComponents);
2151  }
2152 
2154  operator[](const unsigned component) const {
2155  typedef typename multicomponent_traits<T_numtype>::T_element T_compType;
2156 
2157  return extractComponent(T_compType(),component,
2158  multicomponent_traits<T_numtype>::numComponents);
2159  }
2160 
2162  operator[](const int component) {
2163  return operator[](static_cast<unsigned>(component));
2164  }
2165 
2167  operator[](const int component) const {
2168  return operator[](static_cast<unsigned>(component));
2169  }
2170 
2172  // Indirection
2174 
2175  template<typename T_indexContainer>
2176  IndirectArray<T_array, T_indexContainer>
2177  operator[](const T_indexContainer& index)
2178  {
2179  return IndirectArray<T_array, T_indexContainer>(*this,
2180  const_cast<T_indexContainer&>(index));
2181  }
2182 
2184  // Assignment Operators
2186 
2187  // Scalar operand
2188  // NEEDS_WORK : need a precondition check on
2189  // isStorageContiguous when operator, is used.
2191  {
2193  }
2194 
2195  T_array& initialize(T_numtype);
2196 
2197  // Was:
2198  // T_array& operator=(T_numtype);
2199 
2200 #ifdef BZ_NEW_EXPRESSION_TEMPLATES
2201  template<typename T_expr>
2202  T_array& operator=(const ETBase<T_expr>&);
2203  T_array& operator=(const Array<T_numtype,N_rank>&);
2204 
2205  template<typename T> T_array& operator+=(const T&);
2206  template<typename T> T_array& operator-=(const T&);
2207  template<typename T> T_array& operator*=(const T&);
2208  template<typename T> T_array& operator/=(const T&);
2209  template<typename T> T_array& operator%=(const T&);
2210  template<typename T> T_array& operator^=(const T&);
2211  template<typename T> T_array& operator&=(const T&);
2212  template<typename T> T_array& operator|=(const T&);
2213  template<typename T> T_array& operator>>=(const T&);
2214  template<typename T> T_array& operator<<=(const T&);
2215 
2216 #else
2217  T_array& operator+=(T_numtype);
2218  T_array& operator-=(T_numtype);
2219  T_array& operator*=(T_numtype);
2220  T_array& operator/=(T_numtype);
2221  T_array& operator%=(T_numtype);
2222  T_array& operator^=(T_numtype);
2223  T_array& operator&=(T_numtype);
2224  T_array& operator|=(T_numtype);
2225  T_array& operator>>=(T_numtype);
2226  T_array& operator<<=(T_numtype);
2227 
2228  // Array operands
2229  T_array& operator=(const Array<T_numtype,N_rank>&);
2230 
2231  template<typename P_numtype2>
2232  T_array& operator=(const Array<P_numtype2,N_rank>&);
2233  template<typename P_numtype2>
2234  T_array& operator+=(const Array<P_numtype2,N_rank>&);
2235  template<typename P_numtype2>
2236  T_array& operator-=(const Array<P_numtype2,N_rank>&);
2237  template<typename P_numtype2>
2238  T_array& operator*=(const Array<P_numtype2,N_rank>&);
2239  template<typename P_numtype2>
2240  T_array& operator/=(const Array<P_numtype2,N_rank>&);
2241  template<typename P_numtype2>
2242  T_array& operator%=(const Array<P_numtype2,N_rank>&);
2243  template<typename P_numtype2>
2244  T_array& operator^=(const Array<P_numtype2,N_rank>&);
2245  template<typename P_numtype2>
2246  T_array& operator&=(const Array<P_numtype2,N_rank>&);
2247  template<typename P_numtype2>
2248  T_array& operator|=(const Array<P_numtype2,N_rank>&);
2249  template<typename P_numtype2>
2250  T_array& operator>>=(const Array<P_numtype2,N_rank>&);
2251  template<typename P_numtype2>
2252  T_array& operator<<=(const Array<P_numtype2,N_rank>&);
2253 
2254  // Array expression operands
2255  template<typename T_expr>
2256  inline T_array& operator=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2257  template<typename T_expr>
2258  inline T_array& operator+=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2259  template<typename T_expr>
2260  inline T_array& operator-=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2261  template<typename T_expr>
2262  inline T_array& operator*=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2263  template<typename T_expr>
2264  inline T_array& operator/=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2265  template<typename T_expr>
2266  inline T_array& operator%=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2267  template<typename T_expr>
2268  inline T_array& operator^=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2269  template<typename T_expr>
2270  inline T_array& operator&=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2271  template<typename T_expr>
2272  inline T_array& operator|=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2273  template<typename T_expr>
2274  inline T_array& operator>>=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2275  template<typename T_expr>
2276  inline T_array& operator<<=(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr);
2277 
2278  // NEEDS_WORK -- Index placeholder operand
2279 
2280  // NEEDS_WORK -- Random operand
2281 #endif
2282 
2283 public:
2284  // Undocumented implementation routines
2285 
2286  template<typename T_expr, typename T_update>
2287  inline T_array& evaluate(T_expr expr, T_update);
2288 
2289 #ifdef BZ_HAVE_STD
2290 #ifdef BZ_ARRAY_SPACE_FILLING_TRAVERSAL
2291  template<typename T_expr, typename T_update>
2292  inline T_array& evaluateWithFastTraversal(
2293  const TraversalOrder<N_rank - 1>& order,
2294  T_expr expr, T_update);
2295 #endif // BZ_ARRAY_SPACE_FILLING_TRAVERSAL
2296 #endif
2297 
2298 #ifdef BZ_ARRAY_2D_STENCIL_TILING
2299  template<typename T_expr, typename T_update>
2300  inline T_array& evaluateWithTiled2DTraversal(
2301  T_expr expr, T_update);
2302 #endif
2303 
2304  template<typename T_expr, typename T_update>
2305  inline T_array& evaluateWithIndexTraversal1(
2306  T_expr expr, T_update);
2307 
2308  template<typename T_expr, typename T_update>
2309  inline T_array& evaluateWithIndexTraversalN(
2310  T_expr expr, T_update);
2311 
2312  template<typename T_expr, typename T_update>
2313  inline T_array& evaluateWithStackTraversal1(
2314  T_expr expr, T_update);
2315 
2316  template<typename T_expr, typename T_update>
2317  inline T_array& evaluateWithStackTraversalN(
2318  T_expr expr, T_update);
2319 
2320 
2321  T_numtype* restrict getInitializationIterator() { return dataFirst(); }
2322 
2323  bool canCollapse(int outerRank, int innerRank) const {
2324 #ifdef BZ_DEBUG_TRAVERSE
2325  BZ_DEBUG_MESSAGE("stride(" << innerRank << ")=" << stride(innerRank)
2326  << ", extent()=" << extent(innerRank) << ", stride(outerRank)="
2327  << stride(outerRank));
2328 #endif
2329  return (stride(innerRank) * extent(innerRank) == stride(outerRank));
2330  }
2331 
2332 protected:
2334  // Implementation routines
2336 
2337  _bz_inline2 void computeStrides();
2338  _bz_inline2 void setupStorage(int rank);
2339  void constructSubarray(Array<T_numtype, N_rank>& array,
2340  const RectDomain<N_rank>&);
2341  void constructSubarray(Array<T_numtype, N_rank>& array,
2342  const StridedDomain<N_rank>&);
2343  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0);
2344  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0, Range r1);
2345  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
2346  Range r1, Range r2);
2347  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
2348  Range r1, Range r2, Range r3);
2349  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
2350  Range r1, Range r2, Range r3, Range r4);
2351  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
2352  Range r1, Range r2, Range r3, Range r4, Range r5);
2353  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
2354  Range r1, Range r2, Range r3, Range r4, Range r5, Range r6);
2355  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
2356  Range r1, Range r2, Range r3, Range r4, Range r5, Range r6,
2357  Range r7);
2358  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
2359  Range r1, Range r2, Range r3, Range r4, Range r5, Range r6,
2360  Range r7, Range r8);
2361  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
2362  Range r1, Range r2, Range r3, Range r4, Range r5, Range r6,
2363  Range r7, Range r8, Range r9);
2364  void constructSubarray(Array<T_numtype, N_rank>& array, Range r0,
2365  Range r1, Range r2, Range r3, Range r4, Range r5, Range r6,
2366  Range r7, Range r8, Range r9, Range r10);
2367 
2368  void calculateZeroOffset();
2369 
2370  template<int N_rank2, typename R0, typename R1, typename R2, typename R3, typename R4,
2371  typename R5, typename R6, typename R7, typename R8, typename R9, typename R10>
2372  void constructSlice(Array<T_numtype, N_rank2>& array, R0 r0, R1 r1, R2 r2,
2373  R3 r3, R4 r4, R5 r5, R6 r6, R7 r7, R8 r8, R9 r9, R10 r10);
2374 
2375  template<int N_rank2>
2376  void slice(int& setRank, Range r, Array<T_numtype,N_rank2>& array,
2377  TinyVector<int,N_rank2>& rankMap, int sourceRank);
2378 
2379  template<int N_rank2>
2380  void slice(int& setRank, int i, Array<T_numtype,N_rank2>& array,
2381  TinyVector<int,N_rank2>& rankMap, int sourceRank);
2382 
2383  template<int N_rank2>
2384  void slice(int&, nilArraySection, Array<T_numtype,N_rank2>&,
2386  { }
2387 
2388  void doTranspose(int destRank, int sourceRank, T_array& array);
2389 
2390 protected:
2392  // Data members
2394 
2395  // NB: adding new data members may require changes to ctors, reference()
2396 
2397  /*
2398  * For a description of the storage_ members, see the comments for class
2399  * GeneralArrayStorage<N_rank> above.
2400  *
2401  * length_[] contains the extent of each rank. E.g. a 10x20x30 array
2402  * would have length_ = { 10, 20, 30}.
2403  * stride_[] contains the stride to move to the next element along each
2404  * rank.
2405  * zeroOffset_ is the distance from the first element in the array
2406  * to the point (0,0,...,0). If base_ is zero and all ranks are
2407  * stored ascending, then zeroOffset_ is zero. This value
2408  * is needed because to speed up indexing, the data_ member
2409  * (inherited from MemoryBlockReference) always refers to
2410  * (0,0,...,0).
2411  */
2412  GeneralArrayStorage<N_rank> storage_;
2416 };
2417 
2418 /*
2419  * Rank numbers start with zero, which may be confusing to users coming
2420  * from Fortran. To make code more readable, the following constants
2421  * may help. Example: instead of
2422  *
2423  * int firstRankExtent = A.extent(0);
2424  *
2425  * One can write:
2426  *
2427  * int firstRankExtent = A.extent(firstRank);
2428  */
2429 
2430 const int firstRank = 0;
2431 const int secondRank = 1;
2432 const int thirdRank = 2;
2433 const int fourthRank = 3;
2434 const int fifthRank = 4;
2435 const int sixthRank = 5;
2436 const int seventhRank = 6;
2437 const int eighthRank = 7;
2438 const int ninthRank = 8;
2439 const int tenthRank = 9;
2440 const int eleventhRank = 10;
2441 
2442 const int firstDim = 0;
2443 const int secondDim = 1;
2444 const int thirdDim = 2;
2445 const int fourthDim = 3;
2446 const int fifthDim = 4;
2447 const int sixthDim = 5;
2448 const int seventhDim = 6;
2449 const int eighthDim = 7;
2450 const int ninthDim = 8;
2451 const int tenthDim = 9;
2452 const int eleventhDim = 10;
2453 
2454 /*
2455  * Global Functions
2456  */
2457 
2458 template<typename T_numtype>
2459 ostream& operator<<(ostream&, const Array<T_numtype,1>&);
2460 
2461 template<typename T_numtype>
2462 ostream& operator<<(ostream&, const Array<T_numtype,2>&);
2463 
2464 template<typename T_numtype, int N_rank>
2465 ostream& operator<<(ostream&, const Array<T_numtype,N_rank>&);
2466 
2467 template<typename T_numtype, int N_rank>
2468 istream& operator>>(istream& is, Array<T_numtype,N_rank>& x);
2469 
2470 template <typename P_numtype,int N_rank>
2473  a.reference(b);
2474  b.reference(c);
2475 }
2476 
2477 template <typename P_expr>
2479  const _bz_ArrayExpr<P_expr>& expr) {
2480  find(indices,
2481  static_cast< Array<typename P_expr::T_numtype,P_expr::rank> >(expr));
2482 }
2483 
2484 template <typename P_numtype, int N_rank>
2486  const Array<P_numtype,N_rank>& exprVals) {
2487  indices.resize(exprVals.size());
2488  typename Array<P_numtype,N_rank>::const_iterator it, end = exprVals.end();
2489  int j=0;
2490  for (it = exprVals.begin(); it != end; ++it)
2491  if (*it)
2492  indices(j++) = it.position();
2493  if (j)
2494  indices.resizeAndPreserve(j);
2495  else
2496  indices.free();
2497  return;
2498 }
2499 
2500 
2502 
2503 /*
2504  * Include implementations of the member functions and some additional
2505  * global functions.
2506  */
2507 
2508 #include <blitz/array/iter.h> // Array iterators
2509 #include <blitz/array/fastiter.h> // Fast Array iterators (for et)
2510 #include <blitz/array/expr.h> // Array expression objects
2511 #include <blitz/array/methods.cc> // Member functions
2512 #include <blitz/array/eval.cc> // Array expression evaluation
2513 #include <blitz/array/ops.cc> // Assignment operators
2514 #include <blitz/array/io.cc> // Output formatting
2515 #include <blitz/array/et.h> // Expression templates
2516 #include <blitz/array/reduce.h> // Array reduction expression templates
2517 #include <blitz/array/interlace.cc> // Allocation of interlaced arrays
2518 #include <blitz/array/resize.cc> // Array resize, resizeAndPreserve
2519 #include <blitz/array/slicing.cc> // Slicing and subarrays
2520 #include <blitz/array/cycle.cc> // Cycling arrays
2521 #include <blitz/array/complex.cc> // Special support for complex arrays
2522 #include <blitz/array/zip.h> // Zipping multicomponent types
2523 #include <blitz/array/where.h> // where(X,Y,Z)
2524 #include <blitz/array/indirect.h> // Indirection
2525 #include <blitz/array/stencils.h> // Stencil objects
2526 
2527 #endif // BZ_ARRAY_H