Mercurial > octave-nkf
annotate liboctave/mx-inlines.cc @ 8743:1bd918cfb6e2
reimplement any & all using the new reduction code
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Sat, 14 Feb 2009 19:50:43 +0100 |
| parents | 53b4fdeacc2e |
| children | 9f7ce4bf7650 |
| rev | line source |
|---|---|
| 3 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2001, 2002, |
| 4 2003, 2004, 2005, 2006, 2007 John W. Eaton | |
| 3 | 5 |
| 6 This file is part of Octave. | |
| 7 | |
| 8 Octave is free software; you can redistribute it and/or modify it | |
| 9 under the terms of the GNU General Public License as published by the | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 3 | 12 |
| 13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 16 for more details. | |
| 17 | |
| 18 You should have received a copy of the GNU General Public License | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 3 | 21 |
| 22 */ | |
| 23 | |
| 2828 | 24 #if !defined (octave_mx_inlines_h) |
| 25 #define octave_mx_inlines_h 1 | |
| 2804 | 26 |
| 27 #include <cstddef> | |
|
8650
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8380
diff
changeset
|
28 #include <cmath> |
| 2804 | 29 |
| 5525 | 30 #include "quit.h" |
| 31 | |
| 1650 | 32 #include "oct-cmplx.h" |
| 461 | 33 |
|
8380
dbe67764e628
fix & improve speed of diagonal matrix multiplication
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
34 template <class R, class S> |
|
dbe67764e628
fix & improve speed of diagonal matrix multiplication
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
35 inline void |
|
dbe67764e628
fix & improve speed of diagonal matrix multiplication
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
36 mx_inline_fill_vs (R *r, size_t n, S s) |
|
dbe67764e628
fix & improve speed of diagonal matrix multiplication
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
37 { |
|
dbe67764e628
fix & improve speed of diagonal matrix multiplication
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
38 for (size_t i = 0; i < n; i++) |
|
dbe67764e628
fix & improve speed of diagonal matrix multiplication
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
39 r[i] = s; |
|
dbe67764e628
fix & improve speed of diagonal matrix multiplication
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
40 } |
|
dbe67764e628
fix & improve speed of diagonal matrix multiplication
Jaroslav Hajek <highegg@gmail.com>
parents:
7789
diff
changeset
|
41 |
| 2811 | 42 #define VS_OP_FCN(F, OP) \ |
| 43 template <class R, class V, class S> \ | |
| 3262 | 44 inline void \ |
| 2811 | 45 F ## _vs (R *r, const V *v, size_t n, S s) \ |
| 46 { \ | |
| 47 for (size_t i = 0; i < n; i++) \ | |
| 48 r[i] = v[i] OP s; \ | |
| 49 } | |
| 50 | |
| 3769 | 51 VS_OP_FCN (mx_inline_add, +) |
| 52 VS_OP_FCN (mx_inline_subtract, -) | |
| 53 VS_OP_FCN (mx_inline_multiply, *) | |
| 54 VS_OP_FCN (mx_inline_divide, /) | |
| 3 | 55 |
| 2804 | 56 #define VS_OP(F, OP, R, V, S) \ |
| 57 static inline R * \ | |
| 58 F (const V *v, size_t n, S s) \ | |
| 59 { \ | |
| 60 R *r = 0; \ | |
| 61 if (n > 0) \ | |
| 62 { \ | |
| 63 r = new R [n]; \ | |
| 2811 | 64 F ## _vs (r, v, n, s); \ |
| 2804 | 65 } \ |
| 66 return r; \ | |
| 67 } | |
| 3 | 68 |
| 2804 | 69 #define VS_OPS(R, V, S) \ |
| 3769 | 70 VS_OP (mx_inline_add, +, R, V, S) \ |
| 71 VS_OP (mx_inline_subtract, -, R, V, S) \ | |
| 72 VS_OP (mx_inline_multiply, *, R, V, S) \ | |
| 73 VS_OP (mx_inline_divide, /, R, V, S) | |
| 3 | 74 |
| 2804 | 75 VS_OPS (double, double, double) |
| 76 VS_OPS (Complex, double, Complex) | |
| 77 VS_OPS (Complex, Complex, double) | |
| 78 VS_OPS (Complex, Complex, Complex) | |
| 3 | 79 |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
80 VS_OPS (float, float, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
81 VS_OPS (FloatComplex, float, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
82 VS_OPS (FloatComplex, FloatComplex, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
83 VS_OPS (FloatComplex, FloatComplex, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
84 |
| 2811 | 85 #define SV_OP_FCN(F, OP) \ |
| 86 template <class R, class S, class V> \ | |
| 3262 | 87 inline void \ |
| 2811 | 88 F ## _sv (R *r, S s, const V *v, size_t n) \ |
| 89 { \ | |
| 90 for (size_t i = 0; i < n; i++) \ | |
| 91 r[i] = s OP v[i]; \ | |
| 92 } \ | |
| 93 | |
| 3769 | 94 SV_OP_FCN (mx_inline_add, +) |
| 95 SV_OP_FCN (mx_inline_subtract, -) | |
| 96 SV_OP_FCN (mx_inline_multiply, *) | |
| 97 SV_OP_FCN (mx_inline_divide, /) | |
| 2811 | 98 |
| 2804 | 99 #define SV_OP(F, OP, R, S, V) \ |
| 100 static inline R * \ | |
| 101 F (S s, const V *v, size_t n) \ | |
| 102 { \ | |
| 103 R *r = 0; \ | |
| 104 if (n > 0) \ | |
| 105 { \ | |
| 106 r = new R [n]; \ | |
| 2811 | 107 F ## _sv (r, s, v, n); \ |
| 2804 | 108 } \ |
| 109 return r; \ | |
| 110 } | |
| 3 | 111 |
| 2804 | 112 #define SV_OPS(R, S, V) \ |
| 3769 | 113 SV_OP (mx_inline_add, +, R, S, V) \ |
| 114 SV_OP (mx_inline_subtract, -, R, S, V) \ | |
| 115 SV_OP (mx_inline_multiply, *, R, S, V) \ | |
| 116 SV_OP (mx_inline_divide, /, R, S, V) | |
| 3 | 117 |
| 2804 | 118 SV_OPS (double, double, double) |
| 119 SV_OPS (Complex, double, Complex) | |
| 120 SV_OPS (Complex, Complex, double) | |
| 121 SV_OPS (Complex, Complex, Complex) | |
| 3 | 122 |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
123 SV_OPS (float, float, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
124 SV_OPS (FloatComplex, float, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
125 SV_OPS (FloatComplex, FloatComplex, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
126 SV_OPS (FloatComplex, FloatComplex, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
127 |
| 2811 | 128 #define VV_OP_FCN(F, OP) \ |
| 129 template <class R, class T1, class T2> \ | |
| 3262 | 130 inline void \ |
| 2811 | 131 F ## _vv (R *r, const T1 *v1, const T2 *v2, size_t n) \ |
| 132 { \ | |
| 133 for (size_t i = 0; i < n; i++) \ | |
| 134 r[i] = v1[i] OP v2[i]; \ | |
| 135 } \ | |
| 136 | |
| 3769 | 137 VV_OP_FCN (mx_inline_add, +) |
| 138 VV_OP_FCN (mx_inline_subtract, -) | |
| 139 VV_OP_FCN (mx_inline_multiply, *) | |
| 140 VV_OP_FCN (mx_inline_divide, /) | |
| 2811 | 141 |
| 2804 | 142 #define VV_OP(F, OP, R, T1, T2) \ |
| 143 static inline R * \ | |
| 144 F (const T1 *v1, const T2 *v2, size_t n) \ | |
| 145 { \ | |
| 146 R *r = 0; \ | |
| 147 if (n > 0) \ | |
| 148 { \ | |
| 149 r = new R [n]; \ | |
| 2811 | 150 F ## _vv (r, v1, v2, n); \ |
| 2804 | 151 } \ |
| 152 return r; \ | |
| 153 } | |
| 3 | 154 |
| 2804 | 155 #define VV_OPS(R, T1, T2) \ |
| 3769 | 156 VV_OP (mx_inline_add, +, R, T1, T2) \ |
| 157 VV_OP (mx_inline_subtract, -, R, T1, T2) \ | |
| 158 VV_OP (mx_inline_multiply, *, R, T1, T2) \ | |
| 159 VV_OP (mx_inline_divide, /, R, T1, T2) | |
| 3 | 160 |
| 2804 | 161 VV_OPS (double, double, double) |
| 162 VV_OPS (Complex, double, Complex) | |
| 163 VV_OPS (Complex, Complex, double) | |
| 164 VV_OPS (Complex, Complex, Complex) | |
| 3 | 165 |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
166 VV_OPS (float, float, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
167 VV_OPS (FloatComplex, float, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
168 VV_OPS (FloatComplex, FloatComplex, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
169 VV_OPS (FloatComplex, FloatComplex, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
170 |
| 2804 | 171 #define VS_OP2(F, OP, V, S) \ |
| 172 static inline V * \ | |
| 173 F (V *v, size_t n, S s) \ | |
| 174 { \ | |
| 175 for (size_t i = 0; i < n; i++) \ | |
| 176 v[i] OP s; \ | |
| 177 return v; \ | |
| 178 } | |
| 3 | 179 |
| 2804 | 180 #define VS_OP2S(V, S) \ |
| 3769 | 181 VS_OP2 (mx_inline_add2, +=, V, S) \ |
| 182 VS_OP2 (mx_inline_subtract2, -=, V, S) \ | |
| 183 VS_OP2 (mx_inline_multiply2, *=, V, S) \ | |
| 184 VS_OP2 (mx_inline_divide2, /=, V, S) \ | |
| 185 VS_OP2 (mx_inline_copy, =, V, S) | |
| 3 | 186 |
| 2804 | 187 VS_OP2S (double, double) |
| 188 VS_OP2S (Complex, double) | |
| 189 VS_OP2S (Complex, Complex) | |
| 3 | 190 |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
191 VS_OP2S (float, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
192 VS_OP2S (FloatComplex, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
193 VS_OP2S (FloatComplex, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
194 |
| 2804 | 195 #define VV_OP2(F, OP, T1, T2) \ |
| 196 static inline T1 * \ | |
| 197 F (T1 *v1, const T2 *v2, size_t n) \ | |
| 198 { \ | |
| 199 for (size_t i = 0; i < n; i++) \ | |
| 200 v1[i] OP v2[i]; \ | |
| 201 return v1; \ | |
| 202 } | |
| 3 | 203 |
| 2804 | 204 #define VV_OP2S(T1, T2) \ |
| 3769 | 205 VV_OP2 (mx_inline_add2, +=, T1, T2) \ |
| 206 VV_OP2 (mx_inline_subtract2, -=, T1, T2) \ | |
| 207 VV_OP2 (mx_inline_multiply2, *=, T1, T2) \ | |
| 208 VV_OP2 (mx_inline_divide2, /=, T1, T2) \ | |
| 209 VV_OP2 (mx_inline_copy, =, T1, T2) | |
| 3 | 210 |
| 2804 | 211 VV_OP2S (double, double) |
| 212 VV_OP2S (Complex, double) | |
| 213 VV_OP2S (Complex, Complex) | |
| 3 | 214 |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
215 VV_OP2S (float, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
216 VV_OP2S (FloatComplex, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
217 VV_OP2S (FloatComplex, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
218 |
| 2804 | 219 #define OP_EQ_FCN(T1, T2) \ |
| 220 static inline bool \ | |
| 3769 | 221 mx_inline_equal (const T1 *x, const T2 *y, size_t n) \ |
| 2804 | 222 { \ |
| 223 for (size_t i = 0; i < n; i++) \ | |
| 224 if (x[i] != y[i]) \ | |
| 225 return false; \ | |
| 226 return true; \ | |
| 227 } | |
| 3 | 228 |
| 2828 | 229 OP_EQ_FCN (bool, bool) |
| 2804 | 230 OP_EQ_FCN (char, char) |
| 231 OP_EQ_FCN (double, double) | |
| 232 OP_EQ_FCN (Complex, Complex) | |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
233 OP_EQ_FCN (float, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
234 OP_EQ_FCN (FloatComplex, FloatComplex) |
| 3 | 235 |
| 2804 | 236 #define OP_DUP_FCN(OP, F, R, T) \ |
| 237 static inline R * \ | |
| 238 F (const T *x, size_t n) \ | |
| 239 { \ | |
| 240 R *r = 0; \ | |
| 241 if (n > 0) \ | |
| 242 { \ | |
| 243 r = new R [n]; \ | |
| 244 for (size_t i = 0; i < n; i++) \ | |
| 245 r[i] = OP (x[i]); \ | |
| 246 } \ | |
| 247 return r; \ | |
| 248 } | |
| 3 | 249 |
| 3769 | 250 OP_DUP_FCN (, mx_inline_dup, double, double) |
| 251 OP_DUP_FCN (, mx_inline_dup, Complex, Complex) | |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
252 OP_DUP_FCN (, mx_inline_dup, float, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
253 OP_DUP_FCN (, mx_inline_dup, FloatComplex, FloatComplex) |
| 3 | 254 |
| 2804 | 255 // These should really return a bool *. Also, they should probably be |
| 256 // in with a collection of other element-by-element boolean ops. | |
| 3769 | 257 OP_DUP_FCN (0.0 ==, mx_inline_not, double, double) |
| 258 OP_DUP_FCN (0.0 ==, mx_inline_not, double, Complex) | |
| 3 | 259 |
| 3769 | 260 OP_DUP_FCN (, mx_inline_make_complex, Complex, double) |
| 2804 | 261 |
| 3769 | 262 OP_DUP_FCN (-, mx_inline_change_sign, double, double) |
| 263 OP_DUP_FCN (-, mx_inline_change_sign, Complex, Complex) | |
| 3 | 264 |
|
8650
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8380
diff
changeset
|
265 OP_DUP_FCN (std::abs, mx_inline_fabs_dup, double, double) |
|
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8380
diff
changeset
|
266 OP_DUP_FCN (std::abs, mx_inline_cabs_dup, double, Complex) |
| 3769 | 267 OP_DUP_FCN (real, mx_inline_real_dup, double, Complex) |
| 268 OP_DUP_FCN (imag, mx_inline_imag_dup, double, Complex) | |
| 269 OP_DUP_FCN (conj, mx_inline_conj_dup, Complex, Complex) | |
| 3 | 270 |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
271 OP_DUP_FCN (0.0 ==, mx_inline_not, float, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
272 OP_DUP_FCN (static_cast<float>(0.0) ==, mx_inline_not, float, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
273 |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
274 OP_DUP_FCN (, mx_inline_make_complex, FloatComplex, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
275 |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
276 OP_DUP_FCN (-, mx_inline_change_sign, float, float) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
277 OP_DUP_FCN (-, mx_inline_change_sign, FloatComplex, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
278 |
|
8650
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8380
diff
changeset
|
279 OP_DUP_FCN (std::abs, mx_inline_fabs_dup, float, float) |
|
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8380
diff
changeset
|
280 OP_DUP_FCN (std::abs, mx_inline_cabs_dup, float, FloatComplex) |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
281 OP_DUP_FCN (real, mx_inline_real_dup, float, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
282 OP_DUP_FCN (imag, mx_inline_imag_dup, float, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
283 OP_DUP_FCN (conj, mx_inline_conj_dup, FloatComplex, FloatComplex) |
|
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
284 |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
285 // NOTE: std::norm is NOT equivalent |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
286 template <class T> |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
287 inline T cabsq (const std::complex<T>& c) |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
288 { return c.real () * c.real () + c.imag () * c.imag (); } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
289 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
290 #define OP_RED_SUM(ac, el) ac += el |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
291 #define OP_RED_PROD(ac, el) ac *= el |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
292 #define OP_RED_SUMSQ(ac, el) ac += el*el |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
293 #define OP_RED_SUMSQC(ac, el) ac += cabsq (el) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
294 |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
295 // default. works for integers and bool. |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
296 template <class T> |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
297 inline bool xis_true (T x) { return x; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
298 template <class T> |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
299 inline bool xis_false (T x) { return ! x; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
300 // for octave_ints |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
301 template <class T> |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
302 inline bool xis_true (const octave_int<T>& x) { return x.value (); } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
303 template <class T> |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
304 inline bool xis_false (const octave_int<T>& x) { return ! x.value (); } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
305 // for reals, we want to ignore NaNs. |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
306 inline bool xis_true (double x) { return ! xisnan (x) && x != 0.0; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
307 inline bool xis_false (double x) { return x == 0.0; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
308 inline bool xis_true (float x) { return ! xisnan (x) && x != 0.0f; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
309 inline bool xis_false (float x) { return x == 0.0f; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
310 // Ditto for complex. |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
311 inline bool xis_true (const Complex& x) { return ! xisnan (x) && x != 0.0; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
312 inline bool xis_false (const Complex& x) { return x == 0.0; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
313 inline bool xis_true (const FloatComplex& x) { return ! xisnan (x) && x != 0.0f; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
314 inline bool xis_false (const FloatComplex& x) { return x == 0.0f; } |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
315 |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
316 // The following two implement a simple short-circuiting. |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
317 #define OP_RED_ANYC(ac, el) if (xis_true (el)) { ac = true; break; } else continue |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
318 #define OP_RED_ALLC(ac, el) if (xis_false (el)) { ac = false; break; } else continue |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
319 |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
320 // Row any/all reductions are a tradeoff - we traverse the array by |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
321 // columns to gain cache coherence, but sacrifice short-circuiting for that. |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
322 // For certain logical arrays, this could mean a significant loss. |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
323 // A more sophisticated implementation could introduce a buffer of active |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
324 // row indices to achieve both. Right now, I don't see the operation as |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
325 // important enough. |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
326 |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
327 #define OP_RED_ANYR(ac, el) if (xis_true (el)) ac = true |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
328 #define OP_RED_ALLR(ac, el) if (xis_false (el)) ac = false |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
329 |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
330 #define OP_RED_FCN(F, TSRC, TRES, OP, ZERO) \ |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
331 template <class T> \ |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
332 inline TRES \ |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
333 F (const TSRC* v, octave_idx_type n) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
334 { \ |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
335 TRES ac = ZERO; \ |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
336 for (octave_idx_type i = 0; i < n; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
337 OP(ac, v[i]); \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
338 return ac; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
339 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
340 |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
341 OP_RED_FCN (mx_inline_sum, T, T, OP_RED_SUM, 0) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
342 OP_RED_FCN (mx_inline_prod, T, T, OP_RED_PROD, 1) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
343 OP_RED_FCN (mx_inline_sumsq, T, T, OP_RED_SUMSQ, 0) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
344 OP_RED_FCN (mx_inline_sumsq, std::complex<T>, T, OP_RED_SUMSQC, 0) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
345 OP_RED_FCN (mx_inline_any, T, bool, OP_RED_ANYC, false) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
346 OP_RED_FCN (mx_inline_all, T, bool, OP_RED_ALLC, true) |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
347 |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
348 |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
349 #define OP_RED_FCN2(F, TSRC, TRES, OP, ZERO) \ |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
350 template <class T> \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
351 inline void \ |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
352 F (const TSRC* v, TRES *r, octave_idx_type m, octave_idx_type n) \ |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
353 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
354 for (octave_idx_type i = 0; i < m; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
355 r[i] = ZERO; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
356 for (octave_idx_type j = 0; j < n; j++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
357 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
358 for (octave_idx_type i = 0; i < m; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
359 OP(r[i], v[i]); \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
360 v += m; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
361 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
362 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
363 |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
364 OP_RED_FCN2 (mx_inline_sum, T, T, OP_RED_SUM, 0) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
365 OP_RED_FCN2 (mx_inline_prod, T, T, OP_RED_PROD, 1) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
366 OP_RED_FCN2 (mx_inline_sumsq, T, T, OP_RED_SUMSQ, 0) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
367 OP_RED_FCN2 (mx_inline_sumsq, std::complex<T>, T, OP_RED_SUMSQC, 0) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
368 OP_RED_FCN2 (mx_inline_any, T, bool, OP_RED_ANYR, false) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
369 OP_RED_FCN2 (mx_inline_all, T, bool, OP_RED_ALLR, true) |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
370 |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
371 #define OP_RED_FCNN(F, TSRC, TRES) \ |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
372 template <class T> \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
373 inline void \ |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
374 F (const TSRC *v, TRES *r, octave_idx_type l, \ |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
375 octave_idx_type n, octave_idx_type u) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
376 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
377 if (l == 1) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
378 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
379 for (octave_idx_type i = 0; i < u; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
380 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
381 r[i] = F (v, n); \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
382 v += n; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
383 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
384 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
385 else \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
386 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
387 for (octave_idx_type i = 0; i < u; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
388 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
389 F (v, r, l, n); \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
390 v += l*n; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
391 r += l; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
392 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
393 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
394 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
395 |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
396 OP_RED_FCNN (mx_inline_sum, T, T) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
397 OP_RED_FCNN (mx_inline_prod, T, T) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
398 OP_RED_FCNN (mx_inline_sumsq, T, T) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
399 OP_RED_FCNN (mx_inline_sumsq, std::complex<T>, T) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
400 OP_RED_FCNN (mx_inline_any, T, bool) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
401 OP_RED_FCNN (mx_inline_all, T, bool) |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
402 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
403 #define OP_CUM_FCN(F, OP) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
404 template <class T> \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
405 inline void \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
406 F (const T *v, T *r, octave_idx_type n) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
407 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
408 if (n) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
409 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
410 T t = r[0] = v[0]; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
411 for (octave_idx_type i = 1; i < n; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
412 r[i] = t = t OP v[i]; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
413 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
414 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
415 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
416 OP_CUM_FCN (mx_inline_cumsum, +) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
417 OP_CUM_FCN (mx_inline_cumprod, *) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
418 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
419 #define OP_CUM_FCN2(F, OP) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
420 template <class T> \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
421 inline void \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
422 F (const T *v, T *r, octave_idx_type m, octave_idx_type n) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
423 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
424 if (n) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
425 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
426 for (octave_idx_type i = 0; i < m; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
427 r[i] = v[i]; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
428 const T *r0 = r; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
429 for (octave_idx_type j = 1; j < n; j++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
430 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
431 r += m; v += m; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
432 for (octave_idx_type i = 0; i < m; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
433 r[i] = v[i] OP r0[i]; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
434 r0 += m; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
435 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
436 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
437 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
438 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
439 OP_CUM_FCN2 (mx_inline_cumsum, +) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
440 OP_CUM_FCN2 (mx_inline_cumprod, *) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
441 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
442 #define OP_CUM_FCNN(F) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
443 template <class T> \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
444 inline void \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
445 F (const T *v, T *r, octave_idx_type l, \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
446 octave_idx_type n, octave_idx_type u) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
447 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
448 if (l == 1) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
449 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
450 for (octave_idx_type i = 0; i < u; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
451 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
452 F (v, r, n); \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
453 v += n; r += n; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
454 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
455 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
456 else \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
457 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
458 for (octave_idx_type i = 0; i < u; i++) \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
459 { \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
460 F (v, r, l, n); \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
461 v += l*n; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
462 r += l*n; \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
463 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
464 } \ |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
465 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
466 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
467 OP_CUM_FCNN (mx_inline_cumsum) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
468 OP_CUM_FCNN (mx_inline_cumprod) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
469 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
470 // Assistant function |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
471 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
472 inline void |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
473 get_extent_triplet (const dim_vector& dims, int& dim, |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
474 octave_idx_type& l, octave_idx_type& n, |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
475 octave_idx_type& u) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
476 { |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
477 octave_idx_type ndims = dims.length (); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
478 if (dim >= ndims) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
479 { |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
480 l = dims.numel (); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
481 n = 1; |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
482 u = 1; |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
483 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
484 else |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
485 { |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
486 if (dim < 0) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
487 { |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
488 // find first non-singleton dim |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
489 for (dim = 0; dims(dim) == 1 && dim < ndims - 1; dim++) ; |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
490 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
491 // calculate extent triplet. |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
492 l = 1, n = dims(dim), u = 1; |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
493 for (octave_idx_type i = 0; i < dim; i++) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
494 l *= dims (i); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
495 for (octave_idx_type i = dim + 1; i < ndims; i++) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
496 u *= dims (i); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
497 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
498 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
499 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
500 // Appliers. |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
501 // FIXME: is this the best design? C++ gives a lot of options here... |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
502 // maybe it can be done without an explicit parameter? |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
503 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
504 template <class ArrayType, class T> |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
505 inline ArrayType |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
506 do_mx_red_op (const Array<T>& src, int dim, |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
507 void (*mx_red_op) (const T *, typename ArrayType::element_type *, |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
508 octave_idx_type, octave_idx_type, octave_idx_type)) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
509 { |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
510 octave_idx_type l, n, u; |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
511 dim_vector dims = src.dims (); |
|
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
512 // M*b inconsistency: sum([]) = 0 etc. |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
513 if (dims.length () == 2 && dims(0) == 0 && dims(1) == 0) |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
514 dims (1) = 1; |
|
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
515 |
|
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
516 get_extent_triplet (dims, dim, l, n, u); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
517 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
518 // Reduction operation reduces the array size. |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
519 if (dim < dims.length ()) dims(dim) = 1; |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
520 dims.chop_trailing_singletons (); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
521 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
522 ArrayType ret (dims); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
523 mx_red_op (src.data (), ret.fortran_vec (), l, n, u); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
524 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
525 return ret; |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
526 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
527 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
528 template <class ArrayType, class T> |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
529 inline ArrayType |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
530 do_mx_cum_op (const Array<T>& src, int dim, |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
531 void (*mx_cum_op) (const T *, typename ArrayType::element_type *, |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
532 octave_idx_type, octave_idx_type, octave_idx_type)) |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
533 { |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
534 octave_idx_type l, n, u; |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
535 dim_vector dims = src.dims (); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
536 get_extent_triplet (dims, dim, l, n, u); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
537 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
538 // Cumulative operation doesn't reduce the array size. |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
539 ArrayType ret (dims); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
540 mx_cum_op (src.data (), ret.fortran_vec (), l, n, u); |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
541 |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
542 return ret; |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
543 } |
|
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
544 |
| 3864 | 545 // Avoid some code duplication. Maybe we should use templates. |
| 546 | |
| 4015 | 547 #define MX_CUMULATIVE_OP(RET_TYPE, ELT_TYPE, OP) \ |
| 3864 | 548 \ |
| 5275 | 549 octave_idx_type nr = rows (); \ |
| 550 octave_idx_type nc = cols (); \ | |
| 3864 | 551 \ |
| 552 RET_TYPE retval (nr, nc); \ | |
| 553 \ | |
| 554 if (nr > 0 && nc > 0) \ | |
| 555 { \ | |
| 556 if ((nr == 1 && dim == -1) || dim == 1) \ | |
| 557 { \ | |
| 5275 | 558 for (octave_idx_type i = 0; i < nr; i++) \ |
| 3864 | 559 { \ |
| 560 ELT_TYPE t = elem (i, 0); \ | |
| 5275 | 561 for (octave_idx_type j = 0; j < nc; j++) \ |
| 3864 | 562 { \ |
| 563 retval.elem (i, j) = t; \ | |
| 564 if (j < nc - 1) \ | |
| 565 t OP elem (i, j+1); \ | |
| 566 } \ | |
| 567 } \ | |
| 568 } \ | |
| 569 else \ | |
| 570 { \ | |
| 5275 | 571 for (octave_idx_type j = 0; j < nc; j++) \ |
| 3864 | 572 { \ |
| 573 ELT_TYPE t = elem (0, j); \ | |
| 5275 | 574 for (octave_idx_type i = 0; i < nr; i++) \ |
| 3864 | 575 { \ |
| 576 retval.elem (i, j) = t; \ | |
| 577 if (i < nr - 1) \ | |
| 578 t OP elem (i+1, j); \ | |
| 579 } \ | |
| 580 } \ | |
| 581 } \ | |
| 582 } \ | |
| 583 \ | |
| 584 return retval | |
| 585 | |
| 586 #define MX_BASE_REDUCTION_OP(RET_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, \ | |
| 587 MT_RESULT) \ | |
| 588 \ | |
| 5275 | 589 octave_idx_type nr = rows (); \ |
| 590 octave_idx_type nc = cols (); \ | |
| 3864 | 591 \ |
| 592 RET_TYPE retval; \ | |
| 593 \ | |
| 594 if (nr > 0 && nc > 0) \ | |
| 595 { \ | |
| 596 if ((nr == 1 && dim == -1) || dim == 1) \ | |
| 597 { \ | |
| 598 retval.resize (nr, 1); \ | |
| 5275 | 599 for (octave_idx_type i = 0; i < nr; i++) \ |
| 3864 | 600 { \ |
| 601 retval.elem (i, 0) = INIT_VAL; \ | |
| 5275 | 602 for (octave_idx_type j = 0; j < nc; j++) \ |
| 3864 | 603 { \ |
| 604 ROW_EXPR; \ | |
| 605 } \ | |
| 606 } \ | |
| 607 } \ | |
| 608 else \ | |
| 609 { \ | |
| 610 retval.resize (1, nc); \ | |
| 5275 | 611 for (octave_idx_type j = 0; j < nc; j++) \ |
| 3864 | 612 { \ |
| 613 retval.elem (0, j) = INIT_VAL; \ | |
| 5275 | 614 for (octave_idx_type i = 0; i < nr; i++) \ |
| 3864 | 615 { \ |
| 616 COL_EXPR; \ | |
| 617 } \ | |
| 618 } \ | |
| 619 } \ | |
| 620 } \ | |
| 4139 | 621 else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \ |
| 622 retval.resize (1, 1, MT_RESULT); \ | |
| 4015 | 623 else if (nr == 0 && (dim == 0 || dim == -1)) \ |
| 624 retval.resize (1, nc, MT_RESULT); \ | |
| 625 else if (nc == 0 && dim == 1) \ | |
| 626 retval.resize (nr, 1, MT_RESULT); \ | |
| 627 else \ | |
| 4139 | 628 retval.resize (nr > 0, nc > 0); \ |
| 3864 | 629 \ |
| 630 return retval | |
| 631 | |
| 632 #define MX_REDUCTION_OP_ROW_EXPR(OP) \ | |
| 633 retval.elem (i, 0) OP elem (i, j) | |
| 634 | |
| 635 #define MX_REDUCTION_OP_COL_EXPR(OP) \ | |
| 636 retval.elem (0, j) OP elem (i, j) | |
| 637 | |
| 638 #define MX_REDUCTION_OP(RET_TYPE, OP, INIT_VAL, MT_RESULT) \ | |
| 639 MX_BASE_REDUCTION_OP (RET_TYPE, \ | |
| 640 MX_REDUCTION_OP_ROW_EXPR (OP), \ | |
| 641 MX_REDUCTION_OP_COL_EXPR (OP), \ | |
| 642 INIT_VAL, MT_RESULT) | |
| 4015 | 643 |
| 644 #define MX_ANY_ALL_OP_ROW_CODE(TEST_OP, TEST_TRUE_VAL) \ | |
| 645 if (elem (i, j) TEST_OP 0.0) \ | |
| 646 { \ | |
| 647 retval.elem (i, 0) = TEST_TRUE_VAL; \ | |
| 648 break; \ | |
| 649 } | |
| 650 | |
| 651 #define MX_ANY_ALL_OP_COL_CODE(TEST_OP, TEST_TRUE_VAL) \ | |
| 652 if (elem (i, j) TEST_OP 0.0) \ | |
| 653 { \ | |
| 654 retval.elem (0, j) = TEST_TRUE_VAL; \ | |
| 655 break; \ | |
| 656 } | |
| 657 | |
| 658 #define MX_ANY_ALL_OP(DIM, INIT_VAL, TEST_OP, TEST_TRUE_VAL) \ | |
| 659 MX_BASE_REDUCTION_OP (boolMatrix, \ | |
| 660 MX_ANY_ALL_OP_ROW_CODE (TEST_OP, TEST_TRUE_VAL), \ | |
| 661 MX_ANY_ALL_OP_COL_CODE (TEST_OP, TEST_TRUE_VAL), \ | |
| 662 INIT_VAL, INIT_VAL) | |
| 663 | |
| 664 #define MX_ALL_OP(DIM) MX_ANY_ALL_OP (DIM, true, ==, false) | |
| 665 | |
| 666 #define MX_ANY_OP(DIM) MX_ANY_ALL_OP (DIM, false, !=, true) | |
| 667 | |
| 4563 | 668 #define MX_ND_ALL_EXPR elem (iter_idx) == 0 |
| 4556 | 669 |
| 4563 | 670 #define MX_ND_ANY_EXPR elem (iter_idx) != 0 |
| 4556 | 671 |
| 672 #define MX_ND_ALL_EVAL(TEST_EXPR) \ | |
| 5520 | 673 if (retval(result_idx) && (TEST_EXPR)) \ |
| 674 retval(result_idx) = 0; | |
| 4556 | 675 |
| 676 #define MX_ND_ANY_EVAL(TEST_EXPR) \ | |
| 5520 | 677 if (retval(result_idx) || (TEST_EXPR)) \ |
| 678 retval(result_idx) = 1; | |
| 4569 | 679 |
| 5520 | 680 #define MX_ND_REDUCTION(EVAL_EXPR, INIT_VAL, RET_TYPE) \ |
| 4556 | 681 \ |
| 4569 | 682 RET_TYPE retval; \ |
| 4556 | 683 \ |
| 4932 | 684 dim_vector dv = this->dims (); \ |
| 5520 | 685 int nd = this->ndims (); \ |
| 4563 | 686 \ |
| 5955 | 687 int empty = false; \ |
| 4563 | 688 \ |
| 5520 | 689 for (int i = 0; i < nd; i++) \ |
| 4563 | 690 { \ |
| 5955 | 691 if (dv(i) == 0) \ |
| 4563 | 692 { \ |
| 5955 | 693 empty = true; \ |
| 4563 | 694 break; \ |
| 695 } \ | |
| 696 } \ | |
| 697 \ | |
| 5972 | 698 if (nd == 2 && dv(0) == 0 && dv(1) == 0) \ |
| 699 { \ | |
| 700 retval.resize (dim_vector (1, 1), INIT_VAL); \ | |
| 701 return retval; \ | |
| 702 } \ | |
| 703 \ | |
| 5520 | 704 /* We need to find first non-singleton dim. */ \ |
| 705 \ | |
| 706 if (dim == -1) \ | |
| 4556 | 707 { \ |
| 5520 | 708 dim = 0; \ |
| 709 \ | |
| 710 for (int i = 0; i < nd; i++) \ | |
| 4563 | 711 { \ |
| 5520 | 712 if (dv(i) != 1) \ |
| 4563 | 713 { \ |
| 714 dim = i; \ | |
| 715 break; \ | |
| 716 } \ | |
| 717 } \ | |
| 718 } \ | |
| 5520 | 719 else if (dim >= nd) \ |
| 4563 | 720 { \ |
| 5520 | 721 dim = nd++; \ |
| 722 dv.resize (nd, 1); \ | |
| 4563 | 723 } \ |
| 724 \ | |
| 5523 | 725 /* R = op (A, DIM) \ |
| 4563 | 726 \ |
| 5523 | 727 The strategy here is to access the elements of A along the \ |
| 728 dimension specified by DIM. This means that we loop over each \ | |
| 5615 | 729 element of R and adjust the index into A as needed. Store the \ |
| 730 cummulative product of all dimensions of A in CP_SZ. The last \ | |
| 731 element of CP_SZ is the total number of elements of A. */ \ | |
| 4563 | 732 \ |
| 5615 | 733 Array<octave_idx_type> cp_sz (nd+1, 1); \ |
| 734 for (int i = 1; i <= nd; i++) \ | |
| 5523 | 735 cp_sz(i) = cp_sz(i-1)*dv(i-1); \ |
| 5520 | 736 \ |
| 5523 | 737 octave_idx_type reset_at = cp_sz(dim); \ |
| 738 octave_idx_type base_incr = cp_sz(dim+1); \ | |
| 739 octave_idx_type incr = cp_sz(dim); \ | |
| 740 octave_idx_type base = 0; \ | |
| 741 octave_idx_type next_base = base + base_incr; \ | |
| 742 octave_idx_type iter_idx = base; \ | |
| 743 octave_idx_type n_elts = dv(dim); \ | |
| 4556 | 744 \ |
| 5520 | 745 dv(dim) = 1; \ |
| 4556 | 746 \ |
| 5520 | 747 retval.resize (dv, INIT_VAL); \ |
| 4556 | 748 \ |
| 5955 | 749 if (! empty) \ |
| 4556 | 750 { \ |
| 5955 | 751 octave_idx_type nel = dv.numel (); \ |
| 5523 | 752 \ |
| 5955 | 753 octave_idx_type k = 1; \ |
| 754 \ | |
| 755 for (octave_idx_type result_idx = 0; result_idx < nel; result_idx++) \ | |
| 5523 | 756 { \ |
| 5955 | 757 OCTAVE_QUIT; \ |
| 5523 | 758 \ |
| 5955 | 759 for (octave_idx_type j = 0; j < n_elts; j++) \ |
| 760 { \ | |
| 761 OCTAVE_QUIT; \ | |
| 762 \ | |
| 763 EVAL_EXPR; \ | |
| 5520 | 764 \ |
| 5955 | 765 iter_idx += incr; \ |
| 766 } \ | |
| 5523 | 767 \ |
| 5955 | 768 if (k == reset_at) \ |
| 769 { \ | |
| 770 base = next_base; \ | |
| 771 next_base += base_incr; \ | |
| 772 iter_idx = base; \ | |
| 773 k = 1; \ | |
| 774 } \ | |
| 775 else \ | |
| 776 { \ | |
| 777 base++; \ | |
| 778 iter_idx = base; \ | |
| 779 k++; \ | |
| 780 } \ | |
| 781 } \ | |
| 4556 | 782 } \ |
| 783 \ | |
| 4871 | 784 retval.chop_trailing_singletons (); \ |
| 785 \ | |
| 4556 | 786 return retval |
| 4569 | 787 |
| 788 #define MX_ND_REAL_OP_REDUCTION(ASN_EXPR, INIT_VAL) \ | |
| 5520 | 789 MX_ND_REDUCTION (retval(result_idx) ASN_EXPR, INIT_VAL, NDArray) |
| 4569 | 790 |
| 791 #define MX_ND_COMPLEX_OP_REDUCTION(ASN_EXPR, INIT_VAL) \ | |
| 5520 | 792 MX_ND_REDUCTION (retval(result_idx) ASN_EXPR, INIT_VAL, ComplexNDArray) |
| 4569 | 793 |
| 794 #define MX_ND_ANY_ALL_REDUCTION(EVAL_EXPR, VAL) \ | |
| 5520 | 795 MX_ND_REDUCTION (EVAL_EXPR, VAL, boolNDArray) |
| 4556 | 796 |
| 5523 | 797 #define MX_ND_CUMULATIVE_OP(RET_TYPE, ACC_TYPE, INIT_VAL, OP) \ |
| 798 \ | |
| 4584 | 799 RET_TYPE retval; \ |
| 800 \ | |
| 4932 | 801 dim_vector dv = this->dims (); \ |
| 5523 | 802 int nd = this->ndims (); \ |
| 4584 | 803 \ |
| 5955 | 804 bool empty = false; \ |
| 4584 | 805 \ |
| 5523 | 806 for (int i = 0; i < nd; i++) \ |
| 4584 | 807 { \ |
| 5955 | 808 if (dv(i) == 0) \ |
| 4584 | 809 { \ |
| 5955 | 810 empty = true; \ |
| 4584 | 811 break; \ |
| 812 } \ | |
| 813 } \ | |
| 814 \ | |
| 5523 | 815 /* We need to find first non-singleton dim. */ \ |
| 816 \ | |
| 4584 | 817 if (dim == -1) \ |
| 818 { \ | |
| 5523 | 819 dim = 0; \ |
| 820 \ | |
| 821 for (int i = 0; i < nd; i++) \ | |
| 4584 | 822 { \ |
| 5523 | 823 if (dv(i) != 1) \ |
| 4584 | 824 { \ |
| 825 dim = i; \ | |
| 826 break; \ | |
| 827 } \ | |
| 828 } \ | |
| 829 } \ | |
| 5523 | 830 else if (dim >= nd) \ |
| 4584 | 831 { \ |
| 5523 | 832 dim = nd++; \ |
| 833 dv.resize (nd, 1); \ | |
| 4584 | 834 } \ |
| 835 \ | |
| 5523 | 836 /* R = op (A, DIM) \ |
| 4584 | 837 \ |
| 5523 | 838 The strategy here is to access the elements of A along the \ |
| 839 dimension specified by DIM. This means that we loop over each \ | |
| 5611 | 840 element of R and adjust the index into A as needed. Store the \ |
| 5614 | 841 cummulative product of all dimensions of A in CP_SZ. The last \ |
| 842 element of CP_SZ is the total number of elements of A. */ \ | |
| 4584 | 843 \ |
| 5611 | 844 Array<octave_idx_type> cp_sz (nd+1, 1); \ |
| 845 for (int i = 1; i <= nd; i++) \ | |
| 5523 | 846 cp_sz(i) = cp_sz(i-1)*dv(i-1); \ |
| 4584 | 847 \ |
| 5523 | 848 octave_idx_type reset_at = cp_sz(dim); \ |
| 849 octave_idx_type base_incr = cp_sz(dim+1); \ | |
| 850 octave_idx_type incr = cp_sz(dim); \ | |
| 851 octave_idx_type base = 0; \ | |
| 852 octave_idx_type next_base = base + base_incr; \ | |
| 853 octave_idx_type iter_idx = base; \ | |
| 854 octave_idx_type n_elts = dv(dim); \ | |
| 4584 | 855 \ |
| 5523 | 856 retval.resize (dv, INIT_VAL); \ |
| 857 \ | |
| 5955 | 858 if (! empty) \ |
| 859 { \ | |
| 860 octave_idx_type nel = dv.numel () / n_elts; \ | |
| 5523 | 861 \ |
| 5955 | 862 octave_idx_type k = 1; \ |
| 4584 | 863 \ |
| 5955 | 864 for (octave_idx_type i = 0; i < nel; i++) \ |
| 865 { \ | |
| 866 OCTAVE_QUIT; \ | |
| 867 \ | |
| 868 ACC_TYPE prev_val = INIT_VAL; \ | |
| 869 \ | |
| 870 for (octave_idx_type j = 0; j < n_elts; j++) \ | |
| 871 { \ | |
| 872 OCTAVE_QUIT; \ | |
| 5523 | 873 \ |
| 5955 | 874 if (j == 0) \ |
| 875 { \ | |
| 876 retval(iter_idx) = elem (iter_idx); \ | |
| 877 prev_val = retval(iter_idx); \ | |
| 878 } \ | |
| 879 else \ | |
| 880 { \ | |
| 881 prev_val = prev_val OP elem (iter_idx); \ | |
| 882 retval(iter_idx) = prev_val; \ | |
| 883 } \ | |
| 5523 | 884 \ |
| 5955 | 885 iter_idx += incr; \ |
| 886 } \ | |
| 4584 | 887 \ |
| 5955 | 888 if (k == reset_at) \ |
| 5523 | 889 { \ |
| 5955 | 890 base = next_base; \ |
| 891 next_base += base_incr; \ | |
| 892 iter_idx = base; \ | |
| 893 k = 1; \ | |
| 4584 | 894 } \ |
| 895 else \ | |
| 5523 | 896 { \ |
| 5955 | 897 base++; \ |
| 898 iter_idx = base; \ | |
| 899 k++; \ | |
| 900 } \ | |
| 4584 | 901 } \ |
| 5523 | 902 } \ |
| 4584 | 903 \ |
| 5523 | 904 retval.chop_trailing_singletons (); \ |
| 905 \ | |
| 4584 | 906 return retval |
| 907 | |
| 2804 | 908 #endif |
| 3 | 909 |
| 910 /* | |
| 911 ;;; Local Variables: *** | |
| 912 ;;; mode: C++ *** | |
| 913 ;;; End: *** | |
| 914 */ |