source: XIOS/dev/dev_olga/src/extern/blitz/include/blitz/array/expr.h @ 1022

Last change on this file since 1022 was 1022, checked in by mhnguyen, 7 years ago
File size: 71.6 KB
Line 
1// -*- C++ -*-
2/***************************************************************************
3 * blitz/array/expr.h     Array<T,N> expression templates
4 *
5 * $Id$
6 *
7 * Copyright (C) 1997-2011 Todd Veldhuizen <tveldhui@acm.org>
8 *
9 * This file is a part of Blitz.
10 *
11 * Blitz is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation, either version 3
14 * of the License, or (at your option) any later version.
15 *
16 * Blitz is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19 * GNU Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with Blitz.  If not, see <http://www.gnu.org/licenses/>.
23 *
24 * Suggestions:          blitz-devel@lists.sourceforge.net
25 * Bugs:                 blitz-support@lists.sourceforge.net   
26 *
27 * For more information, please see the Blitz++ Home Page:
28 *    https://sourceforge.net/projects/blitz/
29 *
30 ****************************************************************************/
31#ifndef BZ_EXPR_H
32#define BZ_EXPR_H
33
34#include <blitz/prettyprint.h>
35#include <blitz/shapecheck.h>
36#include <blitz/numinquire.h>
37#include <blitz/array/domain.h>
38#include <blitz/array/slice.h>
39#include <blitz/bounds.h>
40
41/*
42 * The array expression templates iterator interface is followed by
43 * these classes:
44 *
45 * FastArrayIterator          <blitz/array/fastiter.h>
46 * _bz_ArrayExpr              <blitz/array/expr.h>
47 * _bz_ArrayExprUnaryOp               "
48 * _bz_ArrayExprBinaryOp              "
49 * _bz_ArrayExprTernaryOp             "
50 * _bz_ArrayExprConstant              "
51 * ArrayIndexMapping          <blitz/array/map.h>
52 * _bz_ArrayExprReduce        <blitz/array/reduce.h>
53 * _bz_StencilExpr            <blitz/array/stencil-et.h>
54 *   ... and derived types            "
55 * IndexPlaceholder           <blitz/indexexpr.h>
56 * _bz_ArrayWhere             <blitz/array/where.h>
57 * _bz_FunctorExpr            <blitz/array/functorExpr.h>
58 * _bz_FunctorExpr2                   "
59 * _bz_FunctorExpr3                   "
60 */
61
62BZ_NAMESPACE(blitz)
63
64#define BZ_MAX(a,b) (a)>(b) ? (a) : (b)
65#define BZ_MIN(a,b) (a)<(b) ? (a) : (b)
66
67template<typename T1, typename T2>
68class _bz_ExprPair {
69public:
70    _bz_ExprPair(const T1& a, const T2& b)
71        : first_(a), second_(b)
72    { }
73
74    const T1& first() const
75    { return first_; }
76
77    const T2& second() const
78    { return second_; }
79
80protected:
81    T1 first_;
82    T2 second_;
83};
84
85template<typename T1, typename T2>
86inline _bz_ExprPair<T1,T2> makeExprPair(const T1& a, const T2& b)
87{
88    return _bz_ExprPair<T1,T2>(a,b);
89}
90
91
92/** The main expression template class which is used to wrap all other
93    elements. This makes an expression easy to recognize.  A container
94    that wants to be usable through the expression template machinery
95    must must implement the same interface as this class. Once it
96    does, it will be able to interoperate with the other containers.
97    Note that the expression components (_bz_ArrayExprUnaryOp,
98    _bz_ArrayExprBinaryOp, _bz_ArrayExprReduce, etc.) do not inherit
99    from ETBase. Since all the functions are duplicated in all ET
100    classes, they are documented only in this class. \todo Actually
101    document the interface.  */
102template<typename P_expr>
103class _bz_ArrayExpr
104#ifdef BZ_NEW_EXPRESSION_TEMPLATES
105    : public ETBase<_bz_ArrayExpr<P_expr> >
106#endif
107{
108
109public:
110    typedef P_expr T_expr;
111    typedef _bz_typename T_expr::T_numtype T_numtype;
112  // select return type
113  typedef typename unwrapET<typename T_expr::T_result>::T_unwrapped test;
114  typedef typename selectET<typename T_expr::T_typeprop, 
115                            T_numtype,
116                            _bz_ArrayExpr<test> >::T_selected T_typeprop;
117  typedef typename unwrapET<T_typeprop>::T_unwrapped T_result;
118
119  typedef typename T_expr::T_optype T_optype;
120    typedef T_expr T_ctorArg1;
121    typedef int    T_ctorArg2;    // dummy
122  typedef _bz_ArrayExpr<_bz_typename P_expr::T_range_result> T_range_result;
123
124    static const int 
125        numArrayOperands = T_expr::numArrayOperands,
126        numTVOperands = T_expr::numTVOperands,
127        numTMOperands = T_expr::numTMOperands,
128        numIndexPlaceholders = T_expr::numIndexPlaceholders,
129      minWidth = T_expr::minWidth,
130      maxWidth = T_expr::maxWidth,
131        rank_ = T_expr::rank_;
132
133  // traits class that computes the vectorized return type for vector
134  // width N.
135  template<int N> struct tvresult {
136    typedef _bz_ArrayExpr<typename T_expr::template tvresult<N>::Type> Type;
137  };
138
139    _bz_ArrayExpr(const _bz_ArrayExpr<T_expr>& a)
140#ifdef BZ_NEW_EXPRESSION_TEMPLATES
141        : ETBase< _bz_ArrayExpr<T_expr> >(a), iter_(a.iter_)
142#else
143        : iter_(a.iter_)
144#endif
145    { }
146
147#if defined(BZ_NEW_EXPRESSION_TEMPLATES) && ! defined(__MWERKS__)
148    template<typename T>
149    _bz_ArrayExpr(const T& a)
150        : iter_(a)
151    { }
152#else
153
154    _bz_ArrayExpr(BZ_ETPARM(T_expr) a)
155        : iter_(a)
156    { }
157#if !defined(__MWERKS__)
158    _bz_ArrayExpr(BZ_ETPARM(_bz_typename T_expr::T_ctorArg1) a)
159        : iter_(a)
160    { }
161#endif
162#endif
163
164    template<typename T1, typename T2>
165    _bz_ArrayExpr(BZ_ETPARM(T1) a, BZ_ETPARM(T2) b)
166        : iter_(a, b)
167    { }
168
169    template<typename T1, typename T2, typename T3>
170    _bz_ArrayExpr(BZ_ETPARM(T1) a, BZ_ETPARM(T2) b, BZ_ETPARM(T3) c)
171        : iter_(a, b, c)
172    { }
173
174    template<typename T1, typename T2, typename T3, typename T4>
175    _bz_ArrayExpr(BZ_ETPARM(T1) a, BZ_ETPARM(T2) b, BZ_ETPARM(T3) c,
176        BZ_ETPARM(T4) d) : iter_(a, b, c, d)
177    { }
178
179    template<typename T1, typename T2>
180    _bz_ArrayExpr(const _bz_ExprPair<T1,T2>& exprpair)
181        : iter_(exprpair.first(), exprpair.second())
182    { }
183
184  T_result operator*() const { return *iter_; }
185
186  T_result first_value() const { return iter_.first_value(); }
187
188#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
189    template<int N_rank>
190    T_result operator()(const TinyVector<int, N_rank> i) const { return iter_(i); }
191#else
192    template<int N_rank>
193    T_result operator()(const TinyVector<int, N_rank>& i) const { return iter_(i); }
194#endif
195
196  template<int N>
197  T_range_result operator()(const RectDomain<N>& d) const
198  {
199    return T_range_result(iter_(d));
200  }
201
202    int ascending(const int rank) const { return iter_.ascending(rank); }
203    int ordering(const int rank)  const { return iter_.ordering(rank);  }
204    int lbound(const int rank)    const { return iter_.lbound(rank);    }
205    int ubound(const int rank)    const { return iter_.ubound(rank);    }
206    RectDomain<rank_> domain() const { return iter_.domain(); }
207
208    void push(int position)
209    { iter_.push(position); }
210
211    void pop(int position)
212    { iter_.pop(position); }
213
214    void advance()
215    { iter_.advance(); }
216
217    void advance(int n)
218    { iter_.advance(n); }
219
220    void loadStride(int rank)
221    { iter_.loadStride(rank); }
222
223    bool isUnitStride(int rank) const
224    { return iter_.isUnitStride(rank); }
225
226    bool isUnitStride() const
227    { return iter_.isUnitStride(); }
228
229    void advanceUnitStride()
230    { iter_.advanceUnitStride(); }
231
232    bool canCollapse(int outerLoopRank, int innerLoopRank) const
233    { 
234        // BZ_DEBUG_MESSAGE("_bz_ArrayExpr<>::canCollapse()");
235        return iter_.canCollapse(outerLoopRank, innerLoopRank); 
236    }
237
238    T_result operator[](int i) const
239    { return iter_[i]; }
240
241    T_result fastRead(diffType i) const
242    { return iter_.fastRead(i); }
243
244  template<int N>
245  typename tvresult<N>::Type fastRead_tv(diffType i) const
246  { return iter_.template fastRead_tv<N>(i); }
247
248  bool isVectorAligned(diffType offset) const 
249  { return iter_.isVectorAligned(offset); }
250
251    // this is needed for the stencil expression fastRead to work
252    void _bz_offsetData(sizeType i)
253    { iter_._bz_offsetData(i); }
254
255    // and these are needed for stencil expression shift to work
256    void _bz_offsetData(sizeType offset, int dim)
257    { iter_._bz_offsetData(offset, dim);}
258 
259    void _bz_offsetData(sizeType offset1, int dim1, sizeType offset2, int dim2)
260    { iter_._bz_offsetData(offset1, dim1, offset2, dim2);}
261
262    diffType suggestStride(int rank) const
263    { return iter_.suggestStride(rank); }
264
265    bool isStride(int rank, diffType stride) const
266    { return iter_.isStride(rank,stride); }
267
268    void prettyPrint(BZ_STD_SCOPE(string) &str) const
269    {
270        prettyPrintFormat format(true);  // Terse formatting by default
271        iter_.prettyPrint(str, format);
272    }
273
274    void prettyPrint(BZ_STD_SCOPE(string) &str, 
275        prettyPrintFormat& format) const
276    { iter_.prettyPrint(str, format); }
277
278    template<typename T_shape>
279    bool shapeCheck(const T_shape& shape) const
280    { return iter_.shapeCheck(shape); }
281
282    template<int N>
283    void moveTo(const TinyVector<int, N>& i)
284    {
285        iter_.moveTo(i);
286    }
287
288    T_result shift(int offset, int dim) const
289    {
290      return iter_.shift(offset, dim);
291    }
292
293    T_result shift(int offset1, int dim1,int offset2, int dim2) const
294    {
295      return iter_.shift(offset1, dim1, offset2, dim2);
296    }
297
298  // sliceinfo for expressions
299  template<typename T1, typename T2 = nilArraySection, 
300           class T3 = nilArraySection, typename T4 = nilArraySection, 
301           class T5 = nilArraySection, typename T6 = nilArraySection, 
302           class T7 = nilArraySection, typename T8 = nilArraySection, 
303           class T9 = nilArraySection, typename T10 = nilArraySection, 
304           class T11 = nilArraySection>
305  class SliceInfo {
306  public:
307    typedef typename T_expr::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_subexpr;
308    typedef _bz_ArrayExpr<T_subexpr> T_slice;
309};
310
311  // slicing (experimental)
312
313  // because _bz_ArrayExpr is the "top-level" expression class, it has
314  // the burden of supplying enough variants of the operator to make
315  // it user-friendly. Hence the large numbers that follow
316  template<typename T1>
317  typename SliceInfo<T1>::T_slice
318  operator()(T1 r1) const
319  {
320    return typename SliceInfo<T1>::T_slice
321      (iter_
322       (r1, 
323        nilArraySection(), nilArraySection(), nilArraySection(), nilArraySection(),
324        nilArraySection(), nilArraySection(), nilArraySection(),
325        nilArraySection(), nilArraySection(), nilArraySection()));
326  }
327
328  template<typename T1, typename T2>
329  typename SliceInfo<T1,T2>::T_slice
330  operator()(T1 r1, T2 r2) const
331  {
332    typedef typename SliceInfo<T1,T2>::T_slice slice;
333    return slice(iter_
334                 (r1, r2, nilArraySection(), nilArraySection(), nilArraySection(),
335                  nilArraySection(), nilArraySection(), nilArraySection(),
336                  nilArraySection(), nilArraySection(), nilArraySection()));
337  }
338
339
340    template<typename T1, typename T2, typename T3>
341    typename SliceInfo<T1,T2,T3>::T_slice
342    operator()(T1 r1, T2 r2, T3 r3) const
343    {
344        typedef typename SliceInfo<T1,T2,T3>::T_slice slice;
345        return slice(iter_(r1, r2, r3, nilArraySection(), nilArraySection(), 
346                           nilArraySection(), nilArraySection(), nilArraySection(), 
347                           nilArraySection(), nilArraySection(), nilArraySection()));
348    }
349
350    template<typename T1, typename T2, typename T3, typename T4>
351    typename SliceInfo<T1,T2,T3,T4>::T_slice
352    operator()(T1 r1, T2 r2, T3 r3, T4 r4) const
353    {
354        typedef typename SliceInfo<T1,T2,T3,T4>::T_slice slice;
355        return slice(iter_(r1, r2, r3, r4, nilArraySection(), nilArraySection(),
356                           nilArraySection(), nilArraySection(), nilArraySection(),
357                           nilArraySection(), nilArraySection()));
358    }
359
360    template<typename T1, typename T2, typename T3, typename T4, typename T5>
361    typename SliceInfo<T1,T2,T3,T4,T5>::T_slice
362    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5) const
363    {
364        typedef typename SliceInfo<T1,T2,T3,T4,T5>::T_slice slice;
365        return slice(iter_(r1, r2, r3, r4, r5, nilArraySection(),
366                           nilArraySection(), nilArraySection(), nilArraySection(),
367                           nilArraySection(), nilArraySection()));
368    }
369
370    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6>
371    typename SliceInfo<T1,T2,T3,T4,T5,T6>::T_slice
372    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6) const
373    {
374        typedef typename SliceInfo<T1,T2,T3,T4,T5,T6>::T_slice slice;
375        return slice(iter_(r1, r2, r3, r4, r5, r6, nilArraySection(), 
376                           nilArraySection(), nilArraySection(),
377                           nilArraySection(), nilArraySection()));
378    }
379
380    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
381        typename T7>
382    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7>::T_slice
383    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7) const
384    {
385        typedef typename SliceInfo<T1,T2,T3,T4,T5,T6,T7>::T_slice slice;
386        return slice(iter_(r1, r2, r3, r4, r5, r6, r7, 
387                           nilArraySection(), nilArraySection(),
388                           nilArraySection(), nilArraySection()));
389    }
390
391    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
392        typename T7, typename T8>
393    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8>::T_slice
394    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8) const
395    {
396        typedef typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8>::T_slice slice;
397        return slice(iter_(r1, r2, r3, r4, r5, r6, r7, r8,
398                           nilArraySection(), nilArraySection(), nilArraySection()));
399    }
400
401    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
402        typename T7, typename T8, typename T9>
403    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9>::T_slice
404    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9) const
405    {
406        typedef typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9>::T_slice slice;
407        return slice(iter_(r1, r2, r3, r4, r5, r6, r7, r8, r9, 
408                           nilArraySection(), nilArraySection()));
409    }
410
411    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
412        typename T7, typename T8, typename T9, typename T10>
413    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10>::T_slice
414    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9, T10 r10) const
415    {
416        typedef typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10>::T_slice slice;
417        return slice(iter_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, 
418                           nilArraySection()));
419    }
420
421    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
422        typename T7, typename T8, typename T9, typename T10, typename T11>
423    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
424    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9, T10 r10, T11 r11) const
425    {
426        typedef typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice slice;
427        return slice(iter_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11));
428    }
429
430  // now complete slicings to a scalar, which can't be expressed with the above
431  // unlike for arrays, these expressions are rvalues so we return by value
432    T_numtype operator()(int i0) const
433    { 
434        return iter_(TinyVector<int, 1>(i0));
435    }
436
437  T_numtype operator()(int i0, int i1) const
438    { 
439        return iter_(TinyVector<int, 2>(i0, i1));
440    }
441
442  T_numtype operator()(int i0, int i1, int i2) const
443    { 
444      return iter_(TinyVector<int, 3>(i0, i1, i2));
445    }
446
447  T_numtype operator()(int i0, int i1, int i2, int i3) const
448    { 
449      return iter_(TinyVector<int, 4>(i0, i1, i2, i3));
450    }
451
452  T_numtype operator()(int i0, int i1, int i2, int i3, int i4) const
453    { 
454      return iter_(TinyVector<int, 5>(i0, i1, i2, i3, i4));
455    }
456
457  T_numtype operator()(int i0, int i1, int i2, int i3, int i4, int i5) const
458    { 
459      return iter_(TinyVector<int, 6>(i0, i1, i2, i3, i4, i5));
460    }
461
462  T_numtype operator()(int i0, int i1, int i2, int i3, int i4, int i5,
463                       int i6) const
464    { 
465      return iter_(TinyVector<int, 7>(i0, i1, i2, i3, i4, i5, i6));
466    }
467
468  T_numtype operator()(int i0, int i1, int i2, int i3, int i4, int i5,
469                       int i6, int i7) const
470    { 
471      return iter_(TinyVector<int, 8>(i0, i1, i2, i3, i4, i5, i6, i7));
472    }
473
474  T_numtype operator()(int i0, int i1, int i2, int i3, int i4, int i5,
475                       int i6, int i7, int i8) const
476    { 
477      return iter_(TinyVector<int, 9>(i0, i1, i2, i3, i4, i5, i6, i7, i8));
478    }
479
480  T_numtype operator()(int i0, int i1, int i2, int i3, int i4, int i5,
481                       int i6, int i7, int i8, int i9) const
482    { 
483      return iter_(TinyVector<int, 10>(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9));
484    }
485
486  T_numtype operator()(int i0, int i1, int i2, int i3, int i4, int i5,
487                       int i6, int i7, int i8, int i9, int i10) const
488    { 
489      return iter_(TinyVector<int, 11>(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10));
490    }
491
492protected:
493    _bz_ArrayExpr() { }
494
495    T_expr iter_;
496};
497
498
499template<typename P_expr, typename P_op>
500class _bz_ArrayExprUnaryOp {
501public:
502    typedef P_expr T_expr;
503    typedef P_op T_op;
504    typedef _bz_typename T_expr::T_numtype T_numtype1;
505    typedef _bz_typename T_op::T_numtype T_numtype;
506
507  // select return type
508  typedef typename unwrapET<typename T_expr::T_result>::T_unwrapped test;
509  typedef typename selectET<typename T_expr::T_typeprop, 
510                            T_numtype, 
511                            _bz_ArrayExprUnaryOp<test, T_op> >::T_selected T_typeprop;
512  typedef typename unwrapET<T_typeprop>::T_unwrapped T_result;
513
514  typedef T_numtype T_optype;
515
516    typedef T_expr T_ctorArg1;
517    typedef int    T_ctorArg2;    // dummy
518  typedef _bz_ArrayExprUnaryOp<_bz_typename P_expr::T_range_result,
519                               P_op> T_range_result;
520
521    static const int 
522        numArrayOperands = T_expr::numArrayOperands,
523        numTVOperands = T_expr::numTVOperands,
524        numTMOperands = T_expr::numTMOperands,
525        numIndexPlaceholders = T_expr::numIndexPlaceholders,
526      minWidth = T_expr::minWidth,
527      maxWidth = T_expr::maxWidth,
528        rank_ = T_expr::rank_;
529
530  template<int N> struct tvresult {
531    typedef _bz_ArrayExprUnaryOp<
532      typename T_expr::template tvresult<N>::Type,
533      T_op> Type; 
534  };
535
536    _bz_ArrayExprUnaryOp(const _bz_ArrayExprUnaryOp<T_expr, T_op>& a)
537        : iter_(a.iter_)
538    { }
539
540    _bz_ArrayExprUnaryOp(BZ_ETPARM(T_expr) a)
541        : iter_(a)
542    { }
543
544  /*
545    _bz_ArrayExprUnaryOp(_bz_typename T_expr::T_ctorArg1 a)
546        : iter_(a)
547    { }
548  */
549
550#if BZ_TEMPLATE_CTOR_DOESNT_CAUSE_HAVOC
551    template<typename T1>
552    explicit _bz_ArrayExprUnaryOp(BZ_ETPARM(T1) a)
553        : iter_(a)
554    { }
555#endif
556
557    int ascending(const int rank) const { return iter_.ascending(rank); }
558    int ordering(const int rank)  const { return iter_.ordering(rank);  }
559    int lbound(const int rank)    const { return iter_.lbound(rank);    }
560    int ubound(const int rank)    const { return iter_.ubound(rank);    }
561    RectDomain<rank_> domain() const { return iter_.domain(); }
562
563
564  /* Functions for reading data. Because they must depend on the
565   * result type, they utilize a helper class.
566   */
567
568  // For numtypes, apply operator
569  template<typename T> struct readHelper {
570    static T_result fastRead(const T_expr& iter, diffType i) {
571      return (T_op::apply(iter.fastRead(i))); };
572    static T_result indexop(const T_expr& iter, int i) {
573      return (T_op::apply(iter[i])); };
574    static T_result deref(const T_expr& iter) {
575      return T_op::apply(*iter); }
576    static T_result first_value(const T_expr& iter)  {
577      return T_op::apply(iter.first_value()); }
578    static T_result shift(const T_expr& iter,
579                          int offset, int dim) {
580      return T_op::apply(iter.shift(offset, dim)); }
581    static T_result shift(const T_expr& iter,
582                          int offset1, int dim1, int offset2, int dim2) {
583      return T_op::apply(iter.shift(offset1, dim1, offset2, dim2)); }
584    template<int N_rank>
585#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
586    static T_result indexop(const T_expr& iter, 
587                            const TinyVector<int, N_rank> i) {
588#else
589      static T_result indexop(const T_expr& iter,
590                              const TinyVector<int, N_rank>& i) {
591#endif
592      return T_op::apply(iter(i)); }
593  };
594
595  // For ET types, bypass operator and create expression
596    template<typename T> struct readHelper<ETBase<T> > {
597      static T_result fastRead(const T_expr& iter, diffType i) {
598        return iter.fastRead(i); };
599      static T_result indexop(const T_expr& iter, int i) {
600        return iter[i]; };
601      static T_result deref(const T_expr& iter) {
602        return *iter; }
603    static T_result first_value(const T_expr& iter)  {
604      return iter.first_value(); }
605    static T_result shift(const T_expr& iter, 
606                          int offset, int dim) {
607      return T_result(iter.shift(offset, dim)); }
608    static T_result shift(const T_expr& iter, 
609                          int offset1, int dim1, int offset2, int dim2) {
610      return T_result(iter.shift(offset1, dim1, offset2, dim2)); }
611      template<int N_rank>
612#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
613      static T_result indexop(const T_expr& iter,
614                              const TinyVector<int, N_rank> i) {
615#else
616    static T_result indexop(const T_expr& iter,
617                            const TinyVector<int, N_rank>& i) {
618#endif
619      return iter(i); }
620    };
621
622    T_result fastRead(diffType i) const { 
623      return readHelper<T_typeprop>::fastRead(iter_, i); }
624
625  template<int N>
626  typename tvresult<N>::Type fastRead_tv(diffType i) const
627  { return iter_.template fastRead_tv<N>(i); }
628
629
630    T_result operator[](int i) const { 
631      return readHelper<T_typeprop>::indexop(iter_, i); }
632
633    template<int N_rank>
634#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
635    T_result operator()(const TinyVector<int, N_rank> i) const {
636#else
637      T_result operator()(const TinyVector<int, N_rank>& i) const {
638#endif
639      return readHelper<T_typeprop>::indexop(iter_,i); }
640
641      T_result operator*() const {
642        return readHelper<T_typeprop>::deref(iter_); }
643
644      T_result first_value() const { 
645        return readHelper<T_typeprop>::first_value(iter_); }
646
647    T_result shift(int offset, int dim) const
648    {
649      return readHelper<T_typeprop>::shift(iter_, offset, dim);
650    }
651
652    T_result shift(int offset1, int dim1,int offset2, int dim2) const
653    {
654      return readHelper<T_typeprop>::shift(iter_, 
655                                           offset1, dim1, offset2, dim2);
656    }
657
658      // ****** end reading
659
660  bool isVectorAligned(diffType offset) const 
661  { return iter_.isVectorAligned(offset); }
662
663  template<int N>
664  T_range_result operator()(const RectDomain<N>& d) const
665  {
666    return T_range_result(iter_(d));
667  }
668
669    void push(int position)
670    {
671        iter_.push(position);
672    }
673
674    void pop(int position)
675    {
676        iter_.pop(position);
677    }
678
679    void advance()
680    {
681        iter_.advance();
682    }
683
684    void advance(int n)
685    {
686        iter_.advance(n);
687    }
688
689    void loadStride(int rank)
690    {
691        iter_.loadStride(rank);
692    }
693
694    bool isUnitStride(int rank) const
695    { return iter_.isUnitStride(rank); }
696
697    bool isUnitStride() const
698    { return iter_.isUnitStride(); }
699
700    void advanceUnitStride()
701    {
702        iter_.advanceUnitStride();
703    }
704
705  template<int N>
706  void moveTo(const TinyVector<int,N>& i)
707    {
708        iter_.moveTo(i);
709    }
710
711    bool canCollapse(int outerLoopRank, int innerLoopRank) const
712    { 
713        // BZ_DEBUG_MESSAGE("_bz_ArrayExprUnaryOp<>::canCollapse");
714        return iter_.canCollapse(outerLoopRank, innerLoopRank); 
715    }
716
717  // this is needed for the stencil expression fastRead to work
718  void _bz_offsetData(sizeType i)
719  {
720    iter_._bz_offsetData(i);
721  }
722
723    // and these are needed for stencil expression shift to work
724    void _bz_offsetData(sizeType offset, int dim)
725    { iter_._bz_offsetData(offset, dim);}
726 
727    void _bz_offsetData(sizeType offset1, int dim1, sizeType offset2, int dim2)
728    { iter_._bz_offsetData(offset1, dim1, offset2, dim2);}
729
730    diffType suggestStride(int rank) const
731    { return iter_.suggestStride(rank); }
732
733    bool isStride(int rank, diffType stride) const
734    { return iter_.isStride(rank,stride); }
735
736    void prettyPrint(BZ_STD_SCOPE(string) &str, 
737        prettyPrintFormat& format) const
738    { T_op::prettyPrint(str, format, iter_); }
739
740    template<typename T_shape>
741    bool shapeCheck(const T_shape& shape) const
742    { return iter_.shapeCheck(shape); }
743
744
745  // sliceinfo for expressions
746  template<typename T1, typename T2 = nilArraySection, 
747           class T3 = nilArraySection, typename T4 = nilArraySection, 
748           class T5 = nilArraySection, typename T6 = nilArraySection, 
749           class T7 = nilArraySection, typename T8 = nilArraySection, 
750           class T9 = nilArraySection, typename T10 = nilArraySection, 
751           class T11 = nilArraySection>
752  class SliceInfo {
753  public:
754    typedef typename T_expr::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice1;
755    typedef _bz_ArrayExprUnaryOp<T_slice1, T_op> T_slice;
756};
757
758    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
759        typename T7, typename T8, typename T9, typename T10, typename T11>
760    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
761    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9, T10 r10, T11 r11) const
762    {
763      return typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
764        (iter_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11));
765    }
766
767protected:
768    _bz_ArrayExprUnaryOp() { }
769
770    T_expr iter_;
771};
772
773
774template<typename P_expr1, typename P_expr2, typename P_op>
775class _bz_ArrayExprBinaryOp {
776public:
777    typedef P_expr1 T_expr1;
778    typedef P_expr2 T_expr2;
779    typedef P_op T_op;
780    typedef _bz_typename T_expr1::T_numtype T_numtype1;
781    typedef _bz_typename T_expr2::T_numtype T_numtype2;
782    typedef _bz_typename T_op::T_numtype T_numtype;
783
784  // select return type
785  typedef typename unwrapET<typename T_expr1::T_result>::T_unwrapped T_unwrapped1;
786  typedef typename unwrapET<typename T_expr2::T_result>::T_unwrapped T_unwrapped2;
787  typedef typename selectET2<typename T_expr1::T_typeprop, 
788                             typename T_expr2::T_typeprop, 
789                             T_numtype, 
790                             _bz_ArrayExprBinaryOp<typename asExpr<T_unwrapped1>::T_expr, 
791                                                   typename asExpr<T_unwrapped2>::T_expr, 
792                                                   T_op> >::T_selected T_typeprop;
793  typedef typename unwrapET<T_typeprop>::T_unwrapped T_result;
794  typedef typename T_op::T_numtype T_optype;
795
796    typedef T_expr1 T_ctorArg1;
797    typedef T_expr2 T_ctorArg2;
798  typedef _bz_ArrayExprBinaryOp<_bz_typename P_expr1::T_range_result, 
799                                _bz_typename P_expr2::T_range_result, 
800                                P_op> T_range_result;
801
802    static const int 
803        numArrayOperands = T_expr1::numArrayOperands
804                         + T_expr2::numArrayOperands,
805        numTVOperands = T_expr1::numTVOperands +
806      T_expr2::numTVOperands,
807        numTMOperands = T_expr1::numTMOperands +
808      T_expr2::numTMOperands,
809        numIndexPlaceholders = T_expr1::numIndexPlaceholders
810                             + T_expr2::numIndexPlaceholders,
811      minWidth = BZ_MIN(T_expr1::minWidth, T_expr2::minWidth),
812      maxWidth = BZ_MAX(T_expr1::maxWidth, T_expr2::maxWidth),
813      rank_ = BZ_MAX(T_expr1::rank_, T_expr2::rank_);
814
815  template<int N> struct tvresult {
816    typedef _bz_ArrayExprBinaryOp<
817      typename T_expr1::template tvresult<N>::Type,
818      typename T_expr2::template tvresult<N>::Type,
819      T_op> Type; 
820  };
821
822    _bz_ArrayExprBinaryOp(
823        const _bz_ArrayExprBinaryOp<T_expr1, T_expr2, T_op>& a)
824        : iter1_(a.iter1_), iter2_(a.iter2_)
825    { }
826
827    template<typename T1, typename T2>
828    _bz_ArrayExprBinaryOp(BZ_ETPARM(T1) a, BZ_ETPARM(T2) b)
829        : iter1_(a), iter2_(b)
830    { }
831
832
833  /* Functions for reading. Because they must depend on the result
834   * type, they utilize a helper class.
835   */
836
837  // For numtypes, apply operator
838  template<typename T> struct readHelper {
839    static T_result fastRead(const T_expr1& iter1, const T_expr2& iter2, 
840                             diffType i) {
841      return T_op::apply(iter1.fastRead(i), iter2.fastRead(i)); }
842    static T_result indexop(const T_expr1& iter1, const T_expr2& iter2, int i) {
843      return T_op::apply(iter1[i], iter2[i]); };
844    static T_result deref(const T_expr1& iter1, const T_expr2& iter2) {
845      return T_op::apply(*iter1, *iter2); }
846    static T_result first_value(const T_expr1& iter1, const T_expr2& iter2)  {
847      return T_op::apply(iter1.first_value(), iter2.first_value()); }
848    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
849                          int offset, int dim) {
850      return T_op::apply(iter1.shift(offset, dim),iter2.shift(offset, dim)); }
851    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
852                          int offset1, int dim1, int offset2, int dim2) {
853      return T_op::apply(iter1.shift(offset1, dim1, offset2, dim2),
854                         iter2.shift(offset1, dim1, offset2, dim2)); }
855    template<int N_rank>
856#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
857    static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
858                            const TinyVector<int, N_rank> i) {
859#else
860      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
861                              const TinyVector<int, N_rank>& i) {
862#endif
863        return T_op::apply(iter1(i), iter2(i)); };
864    };
865   
866  // For ET types, bypass operator and create expression
867    template<typename T> struct readHelper<ETBase<T> > {
868      static T_result fastRead(const T_expr1& iter1, const T_expr2& iter2, 
869                               diffType i) {
870        return T_result(iter1.fastRead(i), iter2.fastRead(i)); }
871    static T_result indexop(const T_expr1& iter1, const T_expr2& iter2, int i) {
872      return T_result(iter1[i], iter2[i]); };
873    static T_result deref(const T_expr1& iter1, const T_expr2& iter2) {
874      return T_result(*iter1, *iter2); }
875    static T_result first_value(const T_expr1& iter1, const T_expr2& iter2)  {
876      return T_result(iter1.first_value(), iter2.first_value()); }
877    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
878                          int offset, int dim) {
879      return T_result(iter1.shift(offset, dim),iter2.shift(offset, dim)); }
880    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
881                          int offset1, int dim1, int offset2, int dim2) {
882      return T_result(iter1.shift(offset1, dim1, offset2, dim2),
883                      iter2.shift(offset1, dim1, offset2, dim2)); }
884      template<int N_rank>
885#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
886      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
887                              const TinyVector<int, N_rank> i) {
888#else
889        static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
890                                const TinyVector<int, N_rank>& i) {
891#endif
892          return T_result(iter1(i), iter2(i)); }
893      };
894
895    T_result fastRead(diffType i) const { 
896      return readHelper<T_typeprop>::fastRead(iter1_, iter2_, i); }
897
898  template<int N>
899  typename tvresult<N>::Type fastRead_tv(diffType i) const
900      { return typename tvresult<N>::Type(iter1_.template fastRead_tv<N>(i),
901                                          iter2_.template fastRead_tv<N>(i)); }
902
903    T_result operator[](int i) const { 
904      return readHelper<T_typeprop>::indexop(iter1_, iter2_, i); }
905
906    template<int N_rank>
907#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
908    T_result operator()(const TinyVector<int, N_rank> i) const {
909#else
910      T_result operator()(const TinyVector<int, N_rank>& i) const {
911#endif
912        return readHelper<T_typeprop>::indexop(iter1_, iter2_, i); }
913
914      T_result operator*() const {
915        return readHelper<T_typeprop>::deref(iter1_, iter2_); }
916
917      T_result first_value() const { 
918        return readHelper<T_typeprop>::first_value(iter1_, iter2_); }
919     
920    T_result shift(int offset, int dim) const {
921      return readHelper<T_typeprop>::shift(iter1_, iter2_, offset, dim); }
922
923    T_result shift(int offset1, int dim1,int offset2, int dim2) const {
924      return readHelper<T_typeprop>::shift(iter1_, iter2_, 
925                                           offset1, dim1, offset2, dim2); }
926   
927
928      // ****** end reading
929
930  bool isVectorAligned(diffType offset) const 
931  { return iter1_.isVectorAligned(offset) && 
932      iter2_.isVectorAligned(offset); }
933
934  template<int N>
935  T_range_result operator()(const RectDomain<N>& d) const
936  {
937    return T_range_result(iter1_(d), iter2_(d));
938  }
939
940    int ascending(const int rank) const {
941        return bounds::compute_ascending(rank, iter1_.ascending(rank), iter2_.ascending(rank));
942    }
943
944    int ordering(const int rank) const {
945        return bounds::compute_ordering(rank, iter1_.ordering(rank), iter2_.ordering(rank));
946    }
947
948    int lbound(const int rank) const { 
949        return bounds::compute_lbound(rank, iter1_.lbound(rank), iter2_.lbound(rank));
950    }
951
952    int ubound(const int rank) const {
953        return bounds::compute_ubound(rank, iter1_.ubound(rank), iter2_.ubound(rank));
954    }
955
956  // defer calculation to lbound/ubound
957  RectDomain<rank_> domain() const 
958  { 
959    TinyVector<int, rank_> lb, ub;
960    for(int r=0; r<rank_; ++r) {
961      lb[r]=lbound(r); ub[r]=ubound(r); 
962    }
963    return RectDomain<rank_>(lb,ub);
964  }
965
966    void push(int position)
967    { 
968        iter1_.push(position); 
969        iter2_.push(position);
970    }
971
972    void pop(int position)
973    { 
974        iter1_.pop(position); 
975        iter2_.pop(position);
976    }
977
978    void advance()
979    { 
980        iter1_.advance(); 
981        iter2_.advance();
982    }
983
984    void advance(int n)
985    {
986        iter1_.advance(n);
987        iter2_.advance(n);
988    }
989
990    void loadStride(int rank)
991    { 
992        iter1_.loadStride(rank); 
993        iter2_.loadStride(rank);
994    }
995   
996    bool isUnitStride(int rank) const
997    { return iter1_.isUnitStride(rank) && iter2_.isUnitStride(rank); }
998
999    bool isUnitStride() const
1000    { return iter1_.isUnitStride() && iter2_.isUnitStride(); }
1001
1002    void advanceUnitStride()
1003    { 
1004        iter1_.advanceUnitStride(); 
1005        iter2_.advanceUnitStride();
1006    }
1007
1008    bool canCollapse(int outerLoopRank, int innerLoopRank) const
1009    { 
1010        // BZ_DEBUG_MESSAGE("_bz_ArrayExprBinaryOp<>::canCollapse");
1011        return iter1_.canCollapse(outerLoopRank, innerLoopRank)
1012            && iter2_.canCollapse(outerLoopRank, innerLoopRank);
1013    } 
1014
1015    // this is needed for the stencil expression fastRead to work
1016    void _bz_offsetData(sizeType i)
1017  {
1018    iter1_._bz_offsetData(i);
1019    iter2_._bz_offsetData(i);
1020  }
1021
1022    // and these are needed for stencil expression shift to work
1023    void _bz_offsetData(sizeType offset, int dim)
1024    { 
1025      iter1_._bz_offsetData(offset, dim);
1026      iter2_._bz_offsetData(offset, dim);
1027    }
1028 
1029    void _bz_offsetData(sizeType offset1, int dim1, sizeType offset2, int dim2)
1030    { 
1031      iter1_._bz_offsetData(offset1, dim1, offset2, dim2);
1032      iter2_._bz_offsetData(offset1, dim1, offset2, dim2);
1033    }
1034
1035    diffType suggestStride(int rank) const
1036    {
1037        diffType stride1 = iter1_.suggestStride(rank);
1038        diffType stride2 = iter2_.suggestStride(rank);
1039        return (stride1 > stride2) ? stride1 : stride2;
1040    }
1041
1042    bool isStride(int rank, diffType stride) const
1043    {
1044        return iter1_.isStride(rank,stride) && iter2_.isStride(rank,stride);
1045    }
1046
1047  template<int N>
1048  void moveTo(const TinyVector<int,N>& i)
1049    {
1050        iter1_.moveTo(i);
1051        iter2_.moveTo(i);
1052    }
1053
1054    void prettyPrint(BZ_STD_SCOPE(string) &str, 
1055        prettyPrintFormat& format) const
1056    {
1057        T_op::prettyPrint(str, format, iter1_, iter2_);
1058    }
1059
1060    template<typename T_shape>
1061    bool shapeCheck(const T_shape& shape) const
1062    { return iter1_.shapeCheck(shape) && iter2_.shapeCheck(shape); }
1063
1064
1065  // sliceinfo for expressions
1066  template<typename T1, typename T2 = nilArraySection, 
1067           class T3 = nilArraySection, typename T4 = nilArraySection, 
1068           class T5 = nilArraySection, typename T6 = nilArraySection, 
1069           class T7 = nilArraySection, typename T8 = nilArraySection, 
1070           class T9 = nilArraySection, typename T10 = nilArraySection, 
1071           class T11 = nilArraySection>
1072  class SliceInfo {
1073  public:
1074    typedef typename T_expr1::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice1;
1075    typedef typename T_expr2::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice2;
1076    typedef _bz_ArrayExprBinaryOp<T_slice1, T_slice2, T_op> T_slice;
1077};
1078
1079    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
1080        typename T7, typename T8, typename T9, typename T10, typename T11>
1081    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
1082    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9, T10 r10, T11 r11) const
1083    {
1084      return typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
1085        (iter1_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11),
1086         iter2_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11));
1087    }
1088
1089protected:
1090    _bz_ArrayExprBinaryOp() { }
1091
1092    T_expr1 iter1_;
1093    T_expr2 iter2_; 
1094};
1095
1096template<typename P_expr1, typename P_expr2, typename P_expr3, typename P_op>
1097class _bz_ArrayExprTernaryOp {
1098public:
1099    typedef P_expr1 T_expr1;
1100    typedef P_expr2 T_expr2;
1101    typedef P_expr3 T_expr3;
1102    typedef P_op T_op;
1103    typedef _bz_typename T_expr1::T_numtype T_numtype1;
1104    typedef _bz_typename T_expr2::T_numtype T_numtype2;
1105    typedef _bz_typename T_expr3::T_numtype T_numtype3;
1106    typedef _bz_typename T_op::T_numtype T_numtype;
1107
1108  // select return type
1109  typedef typename unwrapET<
1110    typename T_expr1::T_result>::T_unwrapped T_unwrapped1;
1111  typedef typename unwrapET<
1112    typename T_expr2::T_result>::T_unwrapped T_unwrapped2;
1113  typedef typename unwrapET<
1114    typename T_expr3::T_result>::T_unwrapped T_unwrapped3;
1115  typedef typename selectET2<typename T_expr1::T_typeprop, 
1116                             typename T_expr2::T_typeprop, 
1117                             T_numtype, 
1118                             char>::T_selected T_intermediary;
1119
1120  typedef typename selectET2<
1121    T_intermediary,
1122    typename T_expr3::T_typeprop, 
1123    T_numtype, 
1124    _bz_ArrayExprTernaryOp<typename asExpr<T_unwrapped1>::T_expr, 
1125                           typename asExpr<T_unwrapped2>::T_expr, 
1126                           typename asExpr<T_unwrapped3>::T_expr, 
1127                           T_op> >::T_selected T_typeprop;
1128  typedef typename unwrapET<T_typeprop>::T_unwrapped T_result;
1129  typedef typename T_op::T_numtype T_optype;
1130
1131    typedef T_expr1 T_ctorArg1;
1132    typedef T_expr2 T_ctorArg2;
1133    typedef T_expr3 T_ctorArg3;
1134  typedef _bz_ArrayExprTernaryOp<_bz_typename P_expr1::T_range_result, 
1135                                 _bz_typename P_expr2::T_range_result, 
1136                                 _bz_typename P_expr3::T_range_result, P_op> T_range_result;
1137
1138    static const int 
1139        numArrayOperands = T_expr1::numArrayOperands
1140                         + T_expr2::numArrayOperands
1141                         + T_expr3::numArrayOperands,
1142        numTVOperands = T_expr1::numTVOperands +
1143      T_expr2::numTVOperands +
1144      T_expr3::numTVOperands,
1145        numTMOperands = T_expr1::numTMOperands +
1146      T_expr2::numTMOperands +
1147      T_expr3::numTMOperands,
1148        numIndexPlaceholders = T_expr1::numIndexPlaceholders
1149                             + T_expr2::numIndexPlaceholders
1150                             + T_expr3::numIndexPlaceholders,
1151      minWidth = BZ_MIN(BZ_MIN(T_expr1::minWidth, T_expr2::minWidth),
1152                        T_expr3::minWidth),
1153      maxWidth = BZ_MAX(BZ_MAX(T_expr1::maxWidth, T_expr2::maxWidth), 
1154                        T_expr3::maxWidth),
1155      rank_ = BZ_MAX(BZ_MAX(T_expr1::rank_, T_expr2::rank_),
1156                     T_expr3::rank_);
1157
1158  template<int N> struct tvresult {
1159    typedef _bz_ArrayExprTernaryOp<
1160      typename T_expr1::template tvresult<N>::Type,
1161      typename T_expr2::template tvresult<N>::Type,
1162      typename T_expr3::template tvresult<N>::Type,
1163      T_op> Type; 
1164  };
1165
1166    _bz_ArrayExprTernaryOp(
1167        const _bz_ArrayExprTernaryOp<T_expr1, T_expr2, T_expr3, T_op>& a)
1168        : iter1_(a.iter1_), iter2_(a.iter2_), iter3_(a.iter3_)
1169    { }
1170
1171    template<typename T1, typename T2, typename T3>
1172    _bz_ArrayExprTernaryOp(BZ_ETPARM(T1) a, BZ_ETPARM(T2) b, BZ_ETPARM(T3) c)
1173        : iter1_(a), iter2_(b), iter3_(c)
1174    { }
1175
1176  /* Functions for reading. Because they must depend on the result
1177   * type, they utilize a helper class.
1178   */
1179
1180  // For numtypes, apply operator
1181  template<typename T> struct readHelper {
1182    static T_result fastRead(const T_expr1& iter1, const T_expr2& iter2, 
1183                             const T_expr3& iter3, diffType i) {
1184      return T_op::apply(iter1.fastRead(i), iter2.fastRead(i), 
1185                         iter3.fastRead(i)); }
1186    static T_result indexop(const T_expr1& iter1, const T_expr2& iter2, 
1187                            const T_expr3& iter3, int i) {
1188      return T_op::apply(iter1[i], iter2[i], iter3[i]); };
1189    static T_result deref(const T_expr1& iter1, const T_expr2& iter2, 
1190                          const T_expr3& iter3) {
1191      return T_op::apply(*iter1, *iter2, *iter3); };
1192    static T_result first_value(const T_expr1& iter1, const T_expr2& iter2,
1193                                const T_expr3& iter3)  {
1194      return T_op::apply(iter1.first_value(), iter2.first_value(),
1195                         iter3.first_value()); }
1196    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
1197                          const T_expr3& iter3, int offset, int dim) {
1198      return T_op::apply(iter1.shift(offset, dim),iter2.shift(offset, dim),
1199                         iter3.shift(offset, dim)); }
1200    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
1201                          const T_expr3& iter3, int offset1, 
1202                          int dim1, int offset2, int dim2) {
1203      return T_op::apply(iter1.shift(offset1, dim1, offset2, dim2),
1204                         iter2.shift(offset1, dim1, offset2, dim2),
1205                         iter3.shift(offset1, dim1, offset2, dim2)); }
1206    template<int N_rank>
1207#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
1208      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
1209                              const T_expr3& iter3, 
1210                              const TinyVector<int, N_rank> i) {
1211#else
1212      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
1213                              const T_expr3& iter3, 
1214                              const TinyVector<int, N_rank>& i) {
1215#endif
1216        return T_op::apply(iter1(i), iter2(i), iter3(i) ); };
1217    };
1218   
1219    // For ET types, bypass operator and create expression
1220    template<typename T> struct readHelper<ETBase<T> > {
1221    static T_result fastRead(const T_expr1& iter1, const T_expr2& iter2, 
1222                             const T_expr3& iter3, diffType i) {
1223      return T_result(iter1.fastRead(i), iter2.fastRead(i), 
1224                      iter3.fastRead(i)); }
1225    static T_result indexop(const T_expr1& iter1, const T_expr2& iter2, 
1226                            const T_expr3& iter3, int i) {
1227      return T_result(iter1[i], iter2[i], iter3[i]); };
1228    static T_result deref(const T_expr1& iter1, const T_expr2& iter2, 
1229                            const T_expr3& iter3) {
1230      return T_result(*iter1, *iter2, *iter3); };
1231    static T_result first_value(const T_expr1& iter1, const T_expr2& iter2,
1232                                const T_expr3& iter3)  {
1233      return T_result(iter1.first_value(), iter2.first_value(),
1234                      iter3.first_value()); }
1235    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
1236                          const T_expr3& iter3, int offset, int dim) {
1237      return T_result(iter1.shift(offset, dim),iter2.shift(offset, dim),
1238                      iter3.shift(offset, dim)); }
1239    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
1240                          const T_expr3& iter3, int offset1, 
1241                          int dim1, int offset2, int dim2) {
1242      return T_result(iter1.shift(offset1, dim1, offset2, dim2),
1243                      iter2.shift(offset1, dim1, offset2, dim2),
1244                      iter3.shift(offset1, dim1, offset2, dim2)); }
1245      template<int N_rank>
1246#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
1247      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
1248                              const T_expr3& iter3, 
1249                              const TinyVector<int, N_rank> i) {
1250#else
1251      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
1252                              const T_expr3& iter3, 
1253                              const TinyVector<int, N_rank>& i) {
1254#endif
1255        return T_result(iter1(i), iter2(i), iter3(i) ); }
1256      };
1257
1258    T_result fastRead(diffType i) const { 
1259      return readHelper<T_typeprop>::fastRead(iter1_, iter2_, iter3_, i); }
1260
1261      template<int N>
1262      typename tvresult<N>::Type fastRead_tv(diffType i) const
1263      { return typename tvresult<N>::Type(iter1_.template fastRead_tv<N>(i),
1264                                          iter2_.template fastRead_tv<N>(i),
1265                                          iter3_.template fastRead_tv<N>(i)); }
1266     
1267    T_result operator[](int i) const { 
1268      return readHelper<T_typeprop>::indexop(iter1_, iter2_, iter3_, i); }
1269
1270    template<int N_rank>
1271#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
1272    T_result operator()(const TinyVector<int, N_rank> i) const {
1273#else
1274      T_result operator()(const TinyVector<int, N_rank>& i) const {
1275#endif
1276        return readHelper<T_typeprop>::indexop(iter1_, iter2_, iter3_, i); }
1277
1278      T_result operator*() const {
1279        return readHelper<T_typeprop>::deref(iter1_, iter2_, iter3_); }
1280   
1281      T_result first_value() const { 
1282        return readHelper<T_typeprop>::first_value(iter1_, iter2_, iter3_); }
1283
1284    T_result shift(int offset, int dim) const {
1285      return readHelper<T_typeprop>::shift(iter1_, iter2_, iter3_, 
1286                                           offset, dim); }
1287
1288    T_result shift(int offset1, int dim1,int offset2, int dim2) const {
1289      return readHelper<T_typeprop>::shift(iter1_, iter2_, iter3_,
1290                                           offset1, dim1, offset2, dim2); }
1291
1292
1293      // ****** end reading
1294
1295  bool isVectorAligned(diffType offset) const 
1296  { return iter1_.isVectorAligned(offset) && 
1297      iter2_.isVectorAligned(offset) &&
1298      iter3_.isVectorAligned(offset); }
1299
1300  template<int N>
1301  T_range_result operator()(const RectDomain<N>& d) const
1302  {
1303    return T_range_result(iter1_(d), iter2_(d), iter3_(d));
1304  }
1305
1306    int ascending(const int rank) const {
1307        return bounds::compute_ascending(rank, bounds::compute_ascending(
1308            rank, iter1_.ascending(rank), iter2_.ascending(rank)),
1309            iter3_.ascending(rank));
1310    }
1311
1312    int ordering(const int rank) const {
1313        return bounds::compute_ordering(rank, bounds::compute_ordering(
1314            rank, iter1_.ordering(rank), iter2_.ordering(rank)),
1315            iter3_.ordering(rank));
1316    }
1317
1318    int lbound(const int rank) const { 
1319        return bounds::compute_lbound(rank, bounds::compute_lbound(
1320            rank, iter1_.lbound(rank), iter2_.lbound(rank)), 
1321            iter3_.lbound(rank));
1322    }
1323
1324    int ubound(const int rank) const {
1325        return bounds::compute_ubound(rank, bounds::compute_ubound(
1326            rank, iter1_.ubound(rank), iter2_.ubound(rank)), 
1327            iter3_.ubound(rank));
1328    }
1329
1330  // defer calculation to lbound/ubound
1331  RectDomain<rank_> domain() const 
1332  { 
1333    TinyVector<int, rank_> lb, ub;
1334    for(int r=0; r<rank_; ++r) {
1335      lb[r]=lbound(r); ub[r]=ubound(r); 
1336    }
1337    return RectDomain<rank_>(lb,ub);
1338  }
1339
1340    void push(int position)
1341    { 
1342        iter1_.push(position); 
1343        iter2_.push(position);
1344        iter3_.push(position);
1345    }
1346
1347    void pop(int position)
1348    { 
1349        iter1_.pop(position); 
1350        iter2_.pop(position);
1351        iter3_.pop(position);
1352    }
1353
1354    void advance()
1355    { 
1356        iter1_.advance(); 
1357        iter2_.advance();
1358        iter3_.advance();
1359    }
1360
1361    void advance(int n)
1362    {
1363        iter1_.advance(n);
1364        iter2_.advance(n);
1365        iter3_.advance(n);
1366    }
1367
1368    void loadStride(int rank)
1369    { 
1370        iter1_.loadStride(rank); 
1371        iter2_.loadStride(rank);
1372        iter3_.loadStride(rank);
1373    }
1374   
1375    bool isUnitStride(int rank) const
1376    {
1377        return iter1_.isUnitStride(rank)
1378            && iter2_.isUnitStride(rank)
1379            && iter3_.isUnitStride(rank);
1380    }
1381
1382    bool isUnitStride() const
1383    {
1384        return iter1_.isUnitStride()
1385            && iter2_.isUnitStride()
1386            && iter3_.isUnitStride();
1387    }
1388
1389    void advanceUnitStride()
1390    { 
1391        iter1_.advanceUnitStride(); 
1392        iter2_.advanceUnitStride();
1393        iter3_.advanceUnitStride();
1394    }
1395
1396    bool canCollapse(int outerLoopRank, int innerLoopRank) const
1397    { 
1398        // BZ_DEBUG_MESSAGE("_bz_ArrayExprTernaryOp<>::canCollapse");
1399        return iter1_.canCollapse(outerLoopRank, innerLoopRank)
1400            && iter2_.canCollapse(outerLoopRank, innerLoopRank)
1401            && iter3_.canCollapse(outerLoopRank, innerLoopRank);
1402    }
1403
1404    // this is needed for the stencil expression fastRead to work
1405    void _bz_offsetData(sizeType i)
1406    {
1407      iter1_._bz_offsetData(i);
1408      iter2_._bz_offsetData(i);
1409      iter3_._bz_offsetData(i);
1410    }
1411
1412    // and these are needed for stencil expression shift to work
1413    void _bz_offsetData(sizeType offset, int dim)
1414    { 
1415      iter1_._bz_offsetData(offset, dim);
1416      iter2_._bz_offsetData(offset, dim);
1417      iter3_._bz_offsetData(offset, dim);
1418    }
1419 
1420    void _bz_offsetData(sizeType offset1, int dim1, sizeType offset2, int dim2)
1421    {
1422      iter1_._bz_offsetData(offset1, dim1, offset2, dim2);
1423      iter2_._bz_offsetData(offset1, dim1, offset2, dim2);
1424      iter3_._bz_offsetData(offset1, dim1, offset2, dim2);
1425    }
1426
1427    diffType suggestStride(int rank) const
1428    {
1429        diffType stride1 = iter1_.suggestStride(rank);
1430        diffType stride2 = iter2_.suggestStride(rank);
1431        diffType stride3 = iter3_.suggestStride(rank);
1432        return stride1 > ( stride2 = (stride2>stride3 ? stride2 : stride3) ) ?
1433            stride1 : stride2;
1434    }
1435
1436    bool isStride(int rank, diffType stride) const
1437    {
1438        return iter1_.isStride(rank,stride)
1439            && iter2_.isStride(rank,stride)
1440            && iter3_.isStride(rank,stride);
1441    }
1442
1443    template<int N>
1444    void moveTo(const TinyVector<int,N>& i)
1445    {
1446        iter1_.moveTo(i);
1447        iter2_.moveTo(i);
1448        iter3_.moveTo(i);
1449    }
1450
1451    void prettyPrint(BZ_STD_SCOPE(string) &str, 
1452        prettyPrintFormat& format) const
1453    {
1454        T_op::prettyPrint(str, format, iter1_, iter2_, iter3_);
1455    }
1456
1457    template<typename T_shape>
1458    bool shapeCheck(const T_shape& shape) const
1459    {
1460        return iter1_.shapeCheck(shape)
1461            && iter2_.shapeCheck(shape)
1462            && iter3_.shapeCheck(shape);
1463    }
1464
1465
1466  // sliceinfo for expressions
1467  template<typename T1, typename T2 = nilArraySection, 
1468           class T3 = nilArraySection, typename T4 = nilArraySection, 
1469           class T5 = nilArraySection, typename T6 = nilArraySection, 
1470           class T7 = nilArraySection, typename T8 = nilArraySection, 
1471           class T9 = nilArraySection, typename T10 = nilArraySection, 
1472           class T11 = nilArraySection>
1473  class SliceInfo {
1474  public:
1475    typedef typename T_expr1::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice1;
1476    typedef typename T_expr2::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice2;
1477    typedef typename T_expr3::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice3;
1478    typedef _bz_ArrayExprTernaryOp<T_slice1, T_slice2, T_slice3, T_op> T_slice;
1479};
1480
1481    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
1482        typename T7, typename T8, typename T9, typename T10, typename T11>
1483    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
1484    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9, T10 r10, T11 r11) const
1485    {
1486      return typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
1487        (iter1_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11),
1488         iter2_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11),
1489         iter3_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11));
1490    }
1491
1492protected:
1493    _bz_ArrayExprTernaryOp() { }
1494
1495    T_expr1 iter1_;
1496    T_expr2 iter2_; 
1497    T_expr3 iter3_; 
1498};
1499
1500template<typename P_expr1, typename P_expr2, typename P_expr3,
1501         typename P_expr4, typename P_op>
1502class _bz_ArrayExprQuaternaryOp {
1503public:
1504    typedef P_expr1 T_expr1;
1505    typedef P_expr2 T_expr2;
1506    typedef P_expr3 T_expr3;
1507    typedef P_expr4 T_expr4;
1508    typedef P_op T_op;
1509    typedef _bz_typename T_expr1::T_numtype T_numtype1;
1510    typedef _bz_typename T_expr2::T_numtype T_numtype2;
1511    typedef _bz_typename T_expr3::T_numtype T_numtype3;
1512    typedef _bz_typename T_expr4::T_numtype T_numtype4;
1513    typedef _bz_typename T_op::T_numtype T_numtype;
1514
1515  // select return type
1516  typedef typename unwrapET<typename T_expr1::T_result>::T_unwrapped T_unwrapped1;
1517  typedef typename unwrapET<typename T_expr2::T_result>::T_unwrapped T_unwrapped2;
1518  typedef typename unwrapET<typename T_expr3::T_result>::T_unwrapped T_unwrapped3;
1519  typedef typename unwrapET<typename T_expr4::T_result>::T_unwrapped T_unwrapped4;
1520  typedef typename selectET2<typename T_expr1::T_typeprop, 
1521                             typename T_expr2::T_typeprop, 
1522                             T_numtype, 
1523                             char>::T_selected T_intermediary1;
1524  typedef typename selectET2<T_intermediary1,
1525                             typename T_expr3::T_typeprop, 
1526                             T_numtype, 
1527                             char>::T_selected T_intermediary2;
1528  typedef typename selectET2<
1529    T_intermediary2,
1530    typename T_expr4::T_typeprop, 
1531    T_numtype, 
1532    _bz_ArrayExprQuaternaryOp<typename asExpr<T_unwrapped1>::T_expr, 
1533                              typename asExpr<T_unwrapped2>::T_expr, 
1534                              typename asExpr<T_unwrapped3>::T_expr, 
1535                              typename asExpr<T_unwrapped4>::T_expr, 
1536                              T_op> >::T_selected T_typeprop;
1537  typedef typename unwrapET<T_typeprop>::T_unwrapped T_result;
1538  typedef typename T_op::T_numtype T_optype;
1539
1540    typedef T_expr1 T_ctorArg1;
1541    typedef T_expr2 T_ctorArg2;
1542    typedef T_expr3 T_ctorArg3;
1543    typedef T_expr4 T_ctorArg4;
1544  typedef _bz_ArrayExprQuaternaryOp<_bz_typename P_expr1::T_range_result, 
1545                                    _bz_typename P_expr2::T_range_result, 
1546                                    _bz_typename P_expr3::T_range_result, 
1547                                    _bz_typename P_expr4::T_range_result, 
1548                                    P_op> T_range_result;
1549
1550  static const int 
1551  numArrayOperands = T_expr1::numArrayOperands
1552    + T_expr2::numArrayOperands
1553    + T_expr3::numArrayOperands
1554    + T_expr4::numArrayOperands,
1555
1556    numTVOperands = T_expr1::numTVOperands +
1557    T_expr2::numTVOperands +
1558    T_expr3::numTVOperands +
1559    T_expr4::numTVOperands,
1560
1561    numTMOperands = T_expr1::numTMOperands +
1562    T_expr2::numTMOperands +
1563    T_expr3::numTMOperands +
1564    T_expr4::numTMOperands,
1565
1566    numIndexPlaceholders = T_expr1::numIndexPlaceholders
1567    + T_expr2::numIndexPlaceholders
1568    + T_expr3::numIndexPlaceholders
1569    + T_expr4::numIndexPlaceholders,
1570
1571    minWidth = BZ_MIN(BZ_MIN(T_expr1::minWidth, T_expr2::minWidth),
1572                      BZ_MIN(T_expr3::minWidth, T_expr4::minWidth)),
1573    maxWidth = BZ_MAX(BZ_MAX(T_expr1::maxWidth, T_expr2::maxWidth),
1574                      BZ_MAX(T_expr3::maxWidth, T_expr4::maxWidth)),
1575
1576    rank_ = BZ_MAX(BZ_MAX(T_expr1::rank_, T_expr2::rank_),
1577                  BZ_MAX(T_expr3::rank_, T_expr4::rank_));
1578
1579  template<int N> struct tvresult {
1580    typedef _bz_ArrayExprQuaternaryOp<
1581      typename T_expr1::template tvresult<N>::Type,
1582      typename T_expr2::template tvresult<N>::Type,
1583      typename T_expr3::template tvresult<N>::Type,
1584      typename T_expr4::template tvresult<N>::Type,
1585      T_op> Type; 
1586  };
1587
1588    _bz_ArrayExprQuaternaryOp(
1589        const _bz_ArrayExprQuaternaryOp<T_expr1, T_expr2, T_expr3, T_expr4, T_op>& a)
1590        : iter1_(a.iter1_), iter2_(a.iter2_), iter3_(a.iter3_), iter4_(a.iter4_)
1591    { }
1592
1593    template<typename T1, typename T2, typename T3, typename T4>
1594    _bz_ArrayExprQuaternaryOp(BZ_ETPARM(T1) a, BZ_ETPARM(T2) b, 
1595                              BZ_ETPARM(T3) c, BZ_ETPARM(T4) d)
1596        : iter1_(a), iter2_(b), iter3_(c), iter4_(d)
1597    { }
1598
1599  /* Functions for reading. Because they must depend on the result
1600   * type, they utilize a helper class.
1601   */
1602
1603  // For numtypes, apply operator
1604  template<typename T> struct readHelper {
1605    static T_result fastRead(const T_expr1& iter1, const T_expr2& iter2, 
1606                             const T_expr3& iter3, const T_expr4& iter4, 
1607                             diffType i) {
1608      return T_op::apply(iter1.fastRead(i), iter2.fastRead(i), 
1609                         iter3.fastRead(i), iter4.fastRead(i)); }
1610    static T_result indexop(const T_expr1& iter1, const T_expr2& iter2, 
1611                            const T_expr3& iter3, const T_expr4& iter4, 
1612                            int i) {
1613      return T_op::apply(iter1[i], iter2[i], iter3[i], iter4[i]); };
1614    static T_result deref(const T_expr1& iter1, const T_expr2& iter2, 
1615                          const T_expr3& iter3, const T_expr4& iter4) {
1616      return T_op::apply(*iter1, *iter2, *iter3, *iter4); };
1617    static T_result first_value(const T_expr1& iter1, const T_expr2& iter2,
1618                                const T_expr3& iter3, const T_expr4& iter4)  {
1619      return T_op::apply(iter1.first_value(), iter2.first_value(),
1620                         iter3.first_value(), iter4.first_value()); }
1621    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
1622                          const T_expr3& iter3, const T_expr4& iter4,
1623                          int offset, int dim) {
1624      return T_op::apply(iter1.shift(offset, dim),iter2.shift(offset, dim),
1625                         iter3.shift(offset, dim), iter4.shift(offset, dim)); }
1626    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
1627                          const T_expr3& iter3, const T_expr4& iter4,
1628                          int offset1, int dim1, int offset2, int dim2) {
1629      return T_op::apply(iter1.shift(offset1, dim1, offset2, dim2),
1630                         iter2.shift(offset1, dim1, offset2, dim2),
1631                         iter3.shift(offset1, dim1, offset2, dim2),
1632                         iter4.shift(offset1, dim1, offset2, dim2)); }
1633    template<int N_rank>
1634#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
1635      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
1636                              const T_expr3& iter3, const T_expr4& iter4, 
1637                              const TinyVector<int, N_rank> i) {
1638#else
1639      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
1640                              const T_expr3& iter3, const T_expr4& iter4, 
1641                              const TinyVector<int, N_rank>& i) {
1642#endif
1643        return T_op::apply(iter1(i), iter2(i), iter3(i), iter4(i) ); };
1644    };
1645   
1646    // For ET types, bypass operator and create expression
1647    template<typename T> struct readHelper<ETBase<T> > {
1648    static T_result fastRead(const T_expr1& iter1, const T_expr2& iter2, 
1649                             const T_expr3& iter3, const T_expr4& iter4, 
1650                             diffType i) {
1651      return T_result(iter1.fastRead(i), iter2.fastRead(i),
1652                      iter3.fastRead(i), iter4.fastRead(i)); }
1653    static T_result indexop(const T_expr1& iter1, const T_expr2& iter2, 
1654                            const T_expr3& iter3, const T_expr4& iter4,
1655                            int i) {
1656      return T_result(iter1[i], iter2[i], iter3[i], iter4[i]); };
1657    static T_result deref(const T_expr1& iter1, const T_expr2& iter2, 
1658                          const T_expr3& iter3, const T_expr4& iter4) {
1659      return T_result(*iter1, *iter2, *iter3, *iter4); };
1660    static T_result first_value(const T_expr1& iter1, const T_expr2& iter2,
1661                                const T_expr3& iter3, const T_expr3& iter4)  {
1662      return T_result(iter1.first_value(), iter2.first_value(),
1663                      iter3.first_value(), iter4.first_value()); }
1664    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
1665                          const T_expr3& iter3, const T_expr4& iter4,
1666                          int offset, int dim) {
1667      return T_result(iter1.shift(offset, dim),iter2.shift(offset, dim),
1668                      iter3.shift(offset, dim), iter4.shift(offset, dim)); }
1669    static T_result shift(const T_expr1& iter1, const T_expr2& iter2, 
1670                          const T_expr3& iter3, const T_expr4& iter4,
1671                          int offset1, int dim1, int offset2, int dim2) {
1672      return T_result(iter1.shift(offset1, dim1, offset2, dim2),
1673                      iter2.shift(offset1, dim1, offset2, dim2),
1674                      iter3.shift(offset1, dim1, offset2, dim2),
1675                      iter4.shift(offset1, dim1, offset2, dim2)); }
1676      template<int N_rank>
1677#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
1678      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
1679                              const T_expr3& iter3, const T_expr4& iter4, 
1680                              const TinyVector<int, N_rank> i) {
1681#else
1682      static T_result indexop(const T_expr1& iter1, const T_expr2& iter2,
1683                              const T_expr3& iter3, const T_expr4& iter4, 
1684                              const TinyVector<int, N_rank>& i) {
1685#endif
1686        return T_result(iter1(i), iter2(i), iter3(i), iter4(i) ); }
1687      };
1688
1689    T_result fastRead(diffType i) const { 
1690      return readHelper<T_typeprop>::fastRead(iter1_, iter2_, 
1691                                              iter3_, iter4_, i); }
1692
1693      template<int N>
1694      typename tvresult<N>::Type fastRead_tv(diffType i) const
1695      { return typename tvresult<N>::Type(iter1_.template fastRead_tv<N>(i),
1696                                          iter2_.template fastRead_tv<N>(i),
1697                                          iter3_.template fastRead_tv<N>(i),
1698                                          iter4_.template fastRead_tv<N>(i)); }
1699
1700    T_result operator[](int i) const { 
1701      return readHelper<T_typeprop>::indexop(iter1_, iter2_, 
1702                                             iter3_, iter4_, i); }
1703
1704    template<int N_rank>
1705#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
1706    T_result operator()(const TinyVector<int, N_rank> i) const {
1707#else
1708      T_result operator()(const TinyVector<int, N_rank>& i) const {
1709#endif
1710        return readHelper<T_typeprop>::indexop(iter1_, iter2_, 
1711                                               iter3_, iter4_, i); }
1712
1713      T_result operator*() const
1714      { return readHelper<T_typeprop>::deref(*iter1_, *iter2_, 
1715                                             *iter3_, *iter4_); }
1716   
1717      T_result first_value() const { 
1718        return readHelper<T_typeprop>::first_value(iter1_, iter2_, 
1719                                                   iter3_, iter4_); }
1720
1721    T_result shift(int offset, int dim) const {
1722      return readHelper<T_typeprop>::shift(iter1_, iter2_, iter3_, 
1723                                           iter4_, offset, dim); }
1724
1725    T_result shift(int offset1, int dim1,int offset2, int dim2) const {
1726      return readHelper<T_typeprop>::shift(iter1_, iter2_, iter3_, iter4_, 
1727                                           offset1, dim1, offset2, dim2); }
1728
1729      // ****** end reading
1730
1731  bool isVectorAligned(diffType offset) const 
1732  { return iter1_.isVectorAligned(offset) &&
1733      iter2_.isVectorAligned(offset) &&
1734      iter3_.isVectorAligned(offset) &&
1735      iter3_.isVectorAligned(offset); }
1736
1737  template<int N>
1738  T_range_result operator()(const RectDomain<rank_>& d) const
1739  {
1740    return T_range_result(iter1_(d), iter2_(d), iter3_(d), iter4_(d));
1741  }
1742
1743    int ascending(const int rank) const {
1744        return bounds::compute_ascending(rank, 
1745                                         bounds::compute_ascending(rank, 
1746                                                                   iter1_.ascending(rank), 
1747                                                                   iter2_.ascending(rank)),
1748                                         bounds::compute_ascending(rank, 
1749                                                                   iter3_.ascending(rank),
1750                                                                   iter4_.ascending(rank)));
1751    }
1752
1753    int ordering(const int rank) const {
1754        return bounds::compute_ordering(rank, 
1755                                        bounds::compute_ordering(rank, 
1756                                                                 iter1_.ordering(rank), 
1757                                                                 iter2_.ordering(rank)),
1758                                        bounds::compute_ordering(rank, 
1759                                                                 iter3_.ordering(rank),
1760                                                                 iter4_.ordering(rank)));
1761    }
1762
1763    int lbound(const int rank) const { 
1764        return bounds::compute_lbound(rank, 
1765                                      bounds::compute_lbound(rank, 
1766                                                             iter1_.lbound(rank), 
1767                                                             iter2_.lbound(rank)),
1768                                      bounds::compute_lbound(rank, 
1769                                                             iter3_.lbound(rank),
1770                                                             iter4_.lbound(rank)));
1771    }
1772
1773    int ubound(const int rank) const {
1774        return bounds::compute_ubound(rank, 
1775                                      bounds::compute_ubound(rank, 
1776                                                             iter1_.ubound(rank), 
1777                                                             iter2_.ubound(rank)),
1778                                      bounds::compute_ubound(rank, 
1779                                                             iter3_.ubound(rank),
1780                                                             iter4_.ubound(rank)));
1781    }
1782
1783  // defer calculation to lbound/ubound
1784  RectDomain<rank_> domain() const 
1785  { 
1786    TinyVector<int, rank_> lb, ub;
1787    for(int r=0; r<rank_; ++r) {
1788      lb[r]=lbound(r); ub[r]=ubound(r); 
1789    }
1790    return RectDomain<rank_>(lb,ub);
1791  }
1792
1793    void push(int position)
1794    { 
1795        iter1_.push(position); 
1796        iter2_.push(position);
1797        iter3_.push(position);
1798        iter4_.push(position);
1799    }
1800
1801    void pop(int position)
1802    { 
1803        iter1_.pop(position); 
1804        iter2_.pop(position);
1805        iter3_.pop(position);
1806        iter4_.pop(position);
1807    }
1808
1809    void advance()
1810    { 
1811        iter1_.advance(); 
1812        iter2_.advance();
1813        iter3_.advance();
1814        iter4_.advance();
1815    }
1816
1817    void advance(int n)
1818    {
1819        iter1_.advance(n);
1820        iter2_.advance(n);
1821        iter3_.advance(n);
1822        iter4_.advance(n);
1823    }
1824
1825    void loadStride(int rank)
1826    { 
1827        iter1_.loadStride(rank); 
1828        iter2_.loadStride(rank);
1829        iter3_.loadStride(rank);
1830        iter4_.loadStride(rank);
1831    }
1832   
1833    bool isUnitStride(int rank) const
1834    {
1835        return iter1_.isUnitStride(rank)
1836            && iter2_.isUnitStride(rank)
1837            && iter3_.isUnitStride(rank)
1838            && iter4_.isUnitStride(rank);
1839    }
1840
1841    bool isUnitStride() const
1842    {
1843        return iter1_.isUnitStride()
1844            && iter2_.isUnitStride()
1845            && iter3_.isUnitStride()
1846            && iter4_.isUnitStride();
1847    }
1848
1849    void advanceUnitStride()
1850    { 
1851        iter1_.advanceUnitStride(); 
1852        iter2_.advanceUnitStride();
1853        iter3_.advanceUnitStride();
1854        iter4_.advanceUnitStride();
1855    }
1856
1857    bool canCollapse(int outerLoopRank, int innerLoopRank) const
1858    { 
1859        // BZ_DEBUG_MESSAGE("_bz_ArrayExprQuaternaryOp<>::canCollapse");
1860        return iter1_.canCollapse(outerLoopRank, innerLoopRank)
1861            && iter2_.canCollapse(outerLoopRank, innerLoopRank)
1862            && iter3_.canCollapse(outerLoopRank, innerLoopRank)
1863            && iter4_.canCollapse(outerLoopRank, innerLoopRank);
1864    }
1865
1866  // this is needed for the stencil expression fastRead to work
1867  void _bz_offsetData(sizeType i)
1868  {
1869    iter1_._bz_offsetData(i);
1870    iter2_._bz_offsetData(i);
1871    iter3_._bz_offsetData(i);
1872    iter4_._bz_offsetData(i);
1873  }
1874
1875    // and these are needed for stencil expression shift to work
1876    void _bz_offsetData(sizeType offset, int dim)
1877    { 
1878      iter1_._bz_offsetData(offset, dim);
1879      iter2_._bz_offsetData(offset, dim);
1880      iter3_._bz_offsetData(offset, dim);
1881      iter4_._bz_offsetData(offset, dim);
1882    }
1883 
1884    void _bz_offsetData(sizeType offset1, int dim1, sizeType offset2, int dim2)
1885    {
1886      iter1_._bz_offsetData(offset1, dim1, offset2, dim2);
1887      iter2_._bz_offsetData(offset1, dim1, offset2, dim2);
1888      iter3_._bz_offsetData(offset1, dim1, offset2, dim2);
1889      iter4_._bz_offsetData(offset1, dim1, offset2, dim2);
1890    }
1891
1892    diffType suggestStride(int rank) const
1893    {
1894        diffType stride1 = iter1_.suggestStride(rank);
1895        diffType stride2 = iter2_.suggestStride(rank);
1896        diffType stride3 = iter3_.suggestStride(rank);
1897        diffType stride4 = iter4_.suggestStride(rank);
1898
1899        //return stride1 > ( stride2 = (stride2>stride3 ? stride2 : stride3) ) ?
1900        //  stride1 : stride2;
1901        return std::max(std::max(stride1, stride2),
1902                        std::max(stride3, stride4));
1903    }
1904
1905    bool isStride(int rank, diffType stride) const
1906    {
1907        return iter1_.isStride(rank,stride)
1908            && iter2_.isStride(rank,stride)
1909            && iter3_.isStride(rank,stride)
1910            && iter4_.isStride(rank,stride);
1911    }
1912
1913    template<int N>
1914    void moveTo(const TinyVector<int,N>& i)
1915    {
1916        iter1_.moveTo(i);
1917        iter2_.moveTo(i);
1918        iter3_.moveTo(i);
1919        iter4_.moveTo(i);
1920    }
1921
1922    void prettyPrint(BZ_STD_SCOPE(string) &str, 
1923        prettyPrintFormat& format) const
1924    {
1925        T_op::prettyPrint(str, format, iter1_, iter2_, iter3_, iter4_);
1926    }
1927
1928    template<typename T_shape>
1929    bool shapeCheck(const T_shape& shape) const
1930    {
1931        return iter1_.shapeCheck(shape)
1932            && iter2_.shapeCheck(shape)
1933            && iter3_.shapeCheck(shape)
1934            && iter4_.shapeCheck(shape);
1935    }
1936
1937
1938  // sliceinfo for expressions
1939  template<typename T1, typename T2 = nilArraySection, 
1940           class T3 = nilArraySection, typename T4 = nilArraySection, 
1941           class T5 = nilArraySection, typename T6 = nilArraySection, 
1942           class T7 = nilArraySection, typename T8 = nilArraySection, 
1943           class T9 = nilArraySection, typename T10 = nilArraySection, 
1944           class T11 = nilArraySection>
1945  class SliceInfo {
1946  public:
1947    typedef typename T_expr1::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice1;
1948    typedef typename T_expr2::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice2;
1949    typedef typename T_expr3::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice3;
1950    typedef typename T_expr4::template SliceInfo<T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11>::T_slice T_slice4;
1951    typedef _bz_ArrayExprQuaternaryOp<T_slice1, T_slice2, T_slice3, T_slice4, T_op> T_slice;
1952};
1953
1954    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
1955        typename T7, typename T8, typename T9, typename T10, typename T11>
1956    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
1957    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9, T10 r10, T11 r11) const
1958    {
1959      return typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
1960        (iter1_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11),
1961         iter2_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11),
1962         iter3_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11),
1963         iter4_(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11));
1964    }
1965
1966protected:
1967    _bz_ArrayExprQuaternaryOp() { }
1968
1969    T_expr1 iter1_;
1970    T_expr2 iter2_; 
1971    T_expr3 iter3_; 
1972    T_expr4 iter4_; 
1973};
1974
1975template<typename P_numtype>
1976class _bz_ArrayExprConstant {
1977public:
1978    typedef P_numtype T_numtype;
1979  typedef typename opType<T_numtype>::T_optype T_optype;
1980  typedef typename asET<T_numtype>::T_wrapped T_typeprop;
1981  typedef typename unwrapET<T_typeprop>::T_unwrapped T_result;
1982
1983    typedef T_numtype T_ctorArg1;
1984    typedef int       T_ctorArg2;    // dummy
1985  typedef _bz_ArrayExprConstant<P_numtype> T_range_result;
1986    static const int 
1987        numArrayOperands = 0, 
1988        numTVOperands = 0, 
1989        numTMOperands = 0, 
1990        numIndexPlaceholders = 0, 
1991      minWidth = simdTypes<T_numtype>::vecWidth,
1992      maxWidth = simdTypes<T_numtype>::vecWidth,
1993        rank_ = 0;
1994
1995  /** For the purpose of vectorizing across the container (as opposed
1996      to for operating on multicomponent types), a constant is always
1997      a constant. */
1998  template<int N> struct tvresult {
1999    typedef _bz_ArrayExprConstant<T_numtype> Type;
2000  };
2001
2002    _bz_ArrayExprConstant(const _bz_ArrayExprConstant<T_numtype>& a)
2003        : value_(a.value_)
2004    { }
2005
2006    _bz_ArrayExprConstant(T_numtype value)
2007        : value_(BZ_NO_PROPAGATE(value))
2008    { 
2009    }
2010
2011    // tiny() and huge() return the smallest and largest representable
2012    // integer values.  See <blitz/numinquire.h>
2013    // NEEDS_WORK: use tiny(int()) once numeric_limits<T> available on
2014    // all platforms
2015
2016    int ascending(const int) const { return INT_MIN; }
2017    int ordering(const int)  const { return INT_MIN; }
2018    int lbound(const int)    const { return INT_MIN; }
2019    int ubound(const int)    const { return INT_MAX; }
2020
2021  // there is no rank so we use highest possible
2022  RectDomain<12> domain() const;
2023
2024    // NEEDS_WORK: use huge(int()) once numeric_limits<T> available on
2025    // all platforms
2026
2027    T_result operator*()   const { return value_; }
2028    T_result first_value() const { return value_; }
2029
2030#ifdef BZ_ARRAY_EXPR_PASS_INDEX_BY_VALUE
2031    template<int N_rank>
2032    T_result operator()(const TinyVector<int,N_rank>) const
2033    { return value_; }
2034#else
2035    template<int N_rank>
2036    T_result operator()(const TinyVector<int,N_rank>&) const
2037    { return value_; }
2038#endif
2039
2040  template<int N_rank>
2041  const _bz_ArrayExprConstant& operator()(const RectDomain<N_rank>& d) const
2042  {
2043    return *this;
2044  }
2045
2046    void push(int) { }
2047    void pop(int) { }
2048    void advance() { }
2049    void advance(int) { }
2050    void loadStride(int) { }
2051
2052    bool isUnitStride(int) const
2053    { return true; }
2054
2055    bool isUnitStride() const
2056    { return true; }
2057
2058    void advanceUnitStride()
2059    { }
2060
2061    bool canCollapse(int,int) const 
2062    { return true; }
2063
2064    T_numtype operator[](int) const
2065    { return value_; }
2066
2067  // this must return by reference, because for multicomponent arrays
2068  // it gets turned into an iterator by the containing expression.
2069  const T_numtype& fastRead(diffType) const
2070    { return value_; }
2071
2072  template<int N>
2073  typename tvresult<N>::Type fastRead_tv(diffType i) const
2074  { return value_; }
2075
2076  bool isVectorAligned(diffType offset) const 
2077  { return true; }
2078
2079  // this is needed for the stencil expression fastRead to work
2080  void _bz_offsetData(sizeType i) const{};
2081
2082    // and these are needed for stencil expression shift to work
2083  void _bz_offsetData(sizeType offset, int dim) const {};
2084 
2085  void _bz_offsetData(sizeType offset1, int dim1, sizeType offset2, int dim2) const {};
2086
2087    diffType suggestStride(int) const
2088    { return 1; }
2089
2090    bool isStride(int,diffType) const
2091    { return true; }
2092
2093    void moveTo(int) const { }
2094
2095    T_result shift(int offset, int dim) const {return value_;}
2096
2097    T_result shift(int offset1, int dim1,int offset2, int dim2) const 
2098    { return value_;}
2099
2100    template<int N_rank>
2101    void moveTo(const TinyVector<int,N_rank>&) const { }
2102
2103    void prettyPrint(BZ_STD_SCOPE(string) &str, 
2104        prettyPrintFormat& format) const
2105    {
2106        if (format.tersePrintingSelected())
2107            str += format.nextScalarOperandSymbol();
2108        else
2109            str += BZ_DEBUG_TEMPLATE_AS_STRING_LITERAL(T_numtype);
2110    }
2111
2112    template<typename T_shape>
2113    bool shapeCheck(const T_shape&) const
2114    { return true; }
2115
2116
2117  // sliceinfo for expressions
2118  template<typename T1, typename T2 = nilArraySection, 
2119           class T3 = nilArraySection, typename T4 = nilArraySection, 
2120           class T5 = nilArraySection, typename T6 = nilArraySection, 
2121           class T7 = nilArraySection, typename T8 = nilArraySection, 
2122           class T9 = nilArraySection, typename T10 = nilArraySection, 
2123           class T11 = nilArraySection>
2124  class SliceInfo {
2125  public:
2126    typedef _bz_ArrayExprConstant<T_numtype> T_slice;
2127};
2128
2129    template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6,
2130        typename T7, typename T8, typename T9, typename T10, typename T11>
2131    typename SliceInfo<T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11>::T_slice
2132    operator()(T1 r1, T2 r2, T3 r3, T4 r4, T5 r5, T6 r6, T7 r7, T8 r8, T9 r9, T10 r10, T11 r11) const
2133    {
2134      return *this;
2135    }
2136
2137protected:
2138    _bz_ArrayExprConstant() { }
2139
2140    T_numtype value_;
2141};
2142   
2143
2144BZ_NAMESPACE_END
2145
2146#include <blitz/array/asexpr.h>
2147
2148#endif // BZ_ARRAYEXPR_H
2149
Note: See TracBrowser for help on using the repository browser.