1 /***************************************************************************
2 * blitz/array/reduce.cc Array reductions.
3 *
4 * $Id$
5 *
6 * Copyright (C) 1997-2011 Todd Veldhuizen <tveldhui@acm.org>
7 *
8 * This file is a part of Blitz.
9 *
10 * Blitz is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation, either version 3
13 * of the License, or (at your option) any later version.
14 *
15 * Blitz is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU Lesser General Public License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with Blitz. If not, see <http://www.gnu.org/licenses/>.
22 *
23 * Suggestions: blitz-devel@lists.sourceforge.net
24 * Bugs: blitz-support@lists.sourceforge.net
25 *
26 * For more information, please see the Blitz++ Home Page:
27 * https://sourceforge.net/projects/blitz/
28 *
29 ****************************************************************************/
30 #ifndef BZ_ARRAYREDUCE_H
31 #error <blitz/array/reduce.cc> must be included via <blitz/array/reduce.h>
32 #endif
33
34 namespace blitz {
35
36 template<typename T_expr, typename T_reduction>
37 _bz_typename T_reduction::T_resulttype
_bz_ArrayExprFullReduce(T_expr expr,T_reduction reduction)38 _bz_ArrayExprFullReduce(T_expr expr, T_reduction reduction)
39 {
40 #ifdef BZ_TAU_PROFILING
41 // Tau profiling code. Provide Tau with a pretty-printed version of
42 // the expression.
43 static std::string exprDescription;
44 if (!exprDescription.length()) // faked static initializer
45 {
46 exprDescription = T_reduction::name();
47 exprDescription += "(";
48 prettyPrintFormat format(true); // Terse mode on
49 expr.prettyPrint(exprDescription, format);
50 exprDescription += ")";
51 }
52 TAU_PROFILE(" ", exprDescription, TAU_BLITZ);
53 #endif // BZ_TAU_PROFILING
54
55 return _bz_reduceWithIndexTraversal(expr, reduction);
56
57 #ifdef BZ_NOT_IMPLEMENTED_FLAG
58 if ((T_expr::numIndexPlaceholders > 0) || (T_reduction::needIndex))
59 {
60 // The expression involves index placeholders, so have to
61 // use index traversal rather than stack traversal.
62 return reduceWithIndexTraversal(expr, reduction);
63 }
64 else {
65 // Use a stack traversal
66 return reduceWithStackTraversal(expr, reduction);
67 }
68 #endif
69 }
70
71 template <typename T_index> struct _bz_IndexingVariant;
72
73 template <>
74 struct _bz_IndexingVariant<int> {
75 template <int N>
indexblitz::_bz_IndexingVariant76 static int index(const TinyVector<int,N>& ind,const int i) {
77 return ind[i];
78 }
79 };
80
81 template <int N>
82 struct _bz_IndexingVariant<TinyVector<int,N> > {
indexblitz::_bz_IndexingVariant83 static const TinyVector<int,N>& index(const TinyVector<int,N>& ind,const int) {
84 return ind;
85 }
86 };
87
88 template<typename T_index, typename T_expr, typename T_reduction>
89 _bz_typename T_reduction::T_resulttype
_bz_reduceWithIndexTraversalGeneric(T_expr expr,T_reduction reduction)90 _bz_reduceWithIndexTraversalGeneric(T_expr expr, T_reduction reduction)
91 {
92 // This is optimized assuming C-style arrays.
93
94 const int rank = T_expr::rank_;
95
96 TinyVector<int,T_expr::rank_> index, first, last;
97
98 unsigned long count = 1;
99
100 for (int i=0; i < rank; ++i) {
101 first(i) = expr.lbound(i);
102 last(i) = expr.ubound(i) + 1;
103 index(i) = first(i);
104 count *= last(i) - first(i);
105 }
106
107 const int maxRank = rank - 1;
108 const int lastlbound = expr.lbound(maxRank);
109 const int lastubound = expr.ubound(maxRank);
110
111 const int lastIndex = lastubound + 1;
112
113 typedef _bz_IndexingVariant<T_index> adapter;
114
115 _bz_ReduceReset<T_reduction::needIndex,T_reduction::needInit> reset;
116 reset(reduction,first,expr);
117
118 while(true) {
119 for (index[maxRank]=lastlbound;index[maxRank]<lastIndex;++index[maxRank])
120 if (!reduction(expr(index),adapter::index(index,maxRank)))
121 return reduction.result(count);
122
123 int j = rank-2;
124 for (;j>=0;--j) {
125 index(j+1) = first(j+1);
126 ++index(j);
127 if (index(j) < last(j))
128 break;
129 }
130
131 if (j<0)
132 return reduction.result(count);
133 }
134 }
135
136 template<typename T_expr, typename T_reduction>
137 _bz_typename T_reduction::T_resulttype
_bz_reduceWithIndexTraversal(T_expr expr,T_reduction reduction)138 _bz_reduceWithIndexTraversal(T_expr expr, T_reduction reduction)
139 {
140 return _bz_reduceWithIndexTraversalGeneric<int>(expr,reduction);
141 }
142
143 // This version is for reductions that require a vector of index positions.
144
145 template<typename T_expr, typename T_reduction>
146 _bz_typename T_reduction::T_resulttype
_bz_reduceWithIndexVectorTraversal(T_expr expr,T_reduction reduction)147 _bz_reduceWithIndexVectorTraversal(T_expr expr, T_reduction reduction)
148 {
149 // We are doing minIndex/maxIndex, so initialize with lower bound
150 return _bz_reduceWithIndexTraversalGeneric<TinyVector<int,T_expr::rank_> >(expr,reduction);
151 }
152
153 }
154
155