GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
mx-inlines.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1993-2025 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if ! defined (octave_mx_inlines_h)
27#define octave_mx_inlines_h 1
28
29// This file should *not* include config.h. It is only included in other C++
30// source files that should have included config.h before including this file.
31
32#include <cstddef>
33#include <cmath>
34
35#include <algorithm>
36
37#include "Array-util.h"
38#include "Array.h"
39#include "bsxfun.h"
40#include "oct-cmplx.h"
41#include "oct-inttypes-fwd.h"
42#include "oct-locbuf.h"
43
44// Provides some commonly repeated, basic loop templates.
45
46template <typename R, typename S>
47inline void
48mx_inline_fill (std::size_t n, R *r, S s)
49{
50 for (std::size_t i = 0; i < n; i++)
51 r[i] = s;
52}
53
54template <typename R, typename X>
55inline void
56mx_inline_uminus (std::size_t n, R *r, const X *x)
57{
58 for (std::size_t i = 0; i < n; i++)
59 r[i] = -x[i];
60}
61
62template <typename R>
63inline void
64mx_inline_uminus2 (std::size_t n, R *r)
65{
66 for (std::size_t i = 0; i < n; i++)
67 r[i] = -r[i];
68}
69
70template <typename X>
71inline void
72mx_inline_iszero (std::size_t n, bool *r, const X *x)
73{
74 const X zero = X ();
75 for (std::size_t i = 0; i < n; i++)
76 r[i] = x[i] == zero;
77}
78
79template <typename X>
80inline void
81mx_inline_notzero (std::size_t n, bool *r, const X *x)
82{
83 const X zero = X ();
84 for (std::size_t i = 0; i < n; i++)
85 r[i] = x[i] != zero;
86}
87
88#define DEFMXBINOP(F, OP) \
89 template <typename R, typename X, typename Y> \
90 inline void F (std::size_t n, R *r, const X *x, const Y *y) \
91 { \
92 for (std::size_t i = 0; i < n; i++) \
93 r[i] = x[i] OP y[i]; \
94 } \
95 template <typename R, typename X, typename Y> \
96 inline void F (std::size_t n, R *r, const X *x, Y y) \
97 { \
98 for (std::size_t i = 0; i < n; i++) \
99 r[i] = x[i] OP y; \
100 } \
101 template <typename R, typename X, typename Y> \
102 inline void F (std::size_t n, R *r, X x, const Y *y) \
103 { \
104 for (std::size_t i = 0; i < n; i++) \
105 r[i] = x OP y[i]; \
106 }
107
108DEFMXBINOP (mx_inline_add, +)
109DEFMXBINOP (mx_inline_sub, -)
110DEFMXBINOP (mx_inline_mul, *)
111DEFMXBINOP (mx_inline_div, /)
112
113#define DEFMXBINOPEQ(F, OP) \
114 template <typename R, typename X> \
115 inline void F (std::size_t n, R *r, const X *x) \
116 { \
117 for (std::size_t i = 0; i < n; i++) \
118 r[i] OP x[i]; \
119 } \
120 template <typename R, typename X> \
121 inline void F (std::size_t n, R *r, X x) \
122 { \
123 for (std::size_t i = 0; i < n; i++) \
124 r[i] OP x; \
125 }
126
127DEFMXBINOPEQ (mx_inline_add2, +=)
128DEFMXBINOPEQ (mx_inline_sub2, -=)
129DEFMXBINOPEQ (mx_inline_mul2, *=)
130DEFMXBINOPEQ (mx_inline_div2, /=)
131
132#define DEFMXCMPOP(F, OP) \
133 template <typename X, typename Y> \
134 inline void F (std::size_t n, bool *r, const X *x, const Y *y) \
135 { \
136 for (std::size_t i = 0; i < n; i++) \
137 r[i] = x[i] OP y[i]; \
138 } \
139 template <typename X, typename Y> \
140 inline void F (std::size_t n, bool *r, const X *x, Y y) \
141 { \
142 for (std::size_t i = 0; i < n; i++) \
143 r[i] = x[i] OP y; \
144 } \
145 template <typename X, typename Y> \
146 inline void F (std::size_t n, bool *r, X x, const Y *y) \
147 { \
148 for (std::size_t i = 0; i < n; i++) \
149 r[i] = x OP y[i]; \
150 }
151
152DEFMXCMPOP (mx_inline_lt, <)
153DEFMXCMPOP (mx_inline_le, <=)
154DEFMXCMPOP (mx_inline_gt, >)
155DEFMXCMPOP (mx_inline_ge, >=)
156DEFMXCMPOP (mx_inline_eq, ==)
157DEFMXCMPOP (mx_inline_ne, !=)
158
159// Convert to logical value, for logical op purposes.
160template <typename T>
161inline bool
162logical_value (T x)
163{
164 return x;
165}
166
167template <typename T>
168inline bool
169logical_value (const std::complex<T>& x)
170{
171 return x.real () != 0 || x.imag () != 0;
172}
173
174template <typename T>
175inline bool
176logical_value (const octave_int<T>& x)
177{
178 return x.value ();
179}
180
181template <typename X>
182void
183mx_inline_not (std::size_t n, bool *r, const X *x)
184{
185 for (std::size_t i = 0; i < n; i++)
186 r[i] = ! logical_value (x[i]);
187}
188
189inline void
190mx_inline_not2 (std::size_t n, bool *r)
191{
192 for (std::size_t i = 0; i < n; i++)
193 r[i] = ! r[i];
194}
195
196#define DEFMXBOOLOP(F, NOT1, OP, NOT2) \
197 template <typename X, typename Y> \
198 inline void F (std::size_t n, bool *r, const X *x, const Y *y) \
199 { \
200 for (std::size_t i = 0; i < n; i++) \
201 r[i] = ((NOT1 logical_value (x[i])) \
202 OP (NOT2 logical_value (y[i]))); \
203 } \
204 template <typename X, typename Y> \
205 inline void F (std::size_t n, bool *r, const X *x, Y y) \
206 { \
207 const bool yy = (NOT2 logical_value (y)); \
208 for (std::size_t i = 0; i < n; i++) \
209 r[i] = (NOT1 logical_value (x[i])) OP yy; \
210 } \
211 template <typename X, typename Y> \
212 inline void F (std::size_t n, bool *r, X x, const Y *y) \
213 { \
214 const bool xx = (NOT1 logical_value (x)); \
215 for (std::size_t i = 0; i < n; i++) \
216 r[i] = xx OP (NOT2 logical_value (y[i])); \
217 }
218
219DEFMXBOOLOP (mx_inline_and, , &&, )
220DEFMXBOOLOP (mx_inline_or, , ||, )
221DEFMXBOOLOP (mx_inline_not_and, !, &&, )
222DEFMXBOOLOP (mx_inline_not_or, !, ||, )
223DEFMXBOOLOP (mx_inline_and_not, , &&, !)
224DEFMXBOOLOP (mx_inline_or_not, , ||, !)
225
226template <typename X>
227inline void
228mx_inline_and2 (std::size_t n, bool *r, const X *x)
229{
230 for (std::size_t i = 0; i < n; i++)
231 r[i] &= logical_value (x[i]);
232}
233
234template <typename X>
235inline void
236mx_inline_and2 (std::size_t n, bool *r, X x)
237{
238 for (std::size_t i = 0; i < n; i++)
239 r[i] &= x;
240}
241
242template <typename X>
243inline void
244mx_inline_or2 (std::size_t n, bool *r, const X *x)
245{
246 for (std::size_t i = 0; i < n; i++)
247 r[i] |= logical_value (x[i]);
248}
249
250template <typename X>
251inline void
252mx_inline_or2 (std::size_t n, bool *r, X x)
253{
254 for (std::size_t i = 0; i < n; i++)
255 r[i] |= x;
256}
257
258template <typename T>
259inline bool
260mx_inline_any_nan (std::size_t n, const T *x)
261{
262 for (std::size_t i = 0; i < n; i++)
263 {
264 if (octave::math::isnan (x[i]))
265 return true;
266 }
267
268 return false;
269}
270
271template <typename T>
272inline bool
273mx_inline_all_finite (std::size_t n, const T *x)
274{
275 for (std::size_t i = 0; i < n; i++)
276 {
277 if (! octave::math::isfinite (x[i]))
278 return false;
279 }
280
281 return true;
282}
283
284template <typename T>
285inline bool
286mx_inline_any_negative (std::size_t n, const T *x)
287{
288 for (std::size_t i = 0; i < n; i++)
289 {
290 if (x[i] < 0)
291 return true;
292 }
293
294 return false;
295}
296
297template <typename T>
298inline bool
299mx_inline_any_positive (std::size_t n, const T *x)
300{
301 for (std::size_t i = 0; i < n; i++)
302 {
303 if (x[i] > 0)
304 return true;
305 }
306
307 return false;
308}
309
310template <typename T>
311inline bool
312mx_inline_all_real (std::size_t n, const std::complex<T> *x)
313{
314 for (std::size_t i = 0; i < n; i++)
315 {
316 if (x[i].imag () != 0)
317 return false;
318 }
319
320 return true;
321}
322
323template <typename T>
324inline void
325mx_inline_real (std::size_t n, T *r, const std::complex<T> *x)
326{
327 for (std::size_t i = 0; i < n; i++)
328 r[i] = x[i].real ();
329}
330
331template <typename T>
332inline void
333mx_inline_imag (std::size_t n, T *r, const std::complex<T> *x)
334{
335 for (std::size_t i = 0; i < n; i++)
336 r[i] = x[i].imag ();
337}
338
339template <typename T>
340inline void
341mx_inline_xmin (std::size_t n, T *r, const T *x, const T *y)
342{
343 for (std::size_t i = 0; i < n; i++)
344 r[i] = octave::math::min (x[i], y[i]);
345}
346
347template <typename T>
348inline void
349mx_inline_xmin (std::size_t n, T *r, const T *x, T y)
350{
351 for (std::size_t i = 0; i < n; i++)
352 r[i] = octave::math::min (x[i], y);
353}
354
355template <typename T>
356inline void
357mx_inline_xmin (std::size_t n, T *r, T x, const T *y)
358{
359 for (std::size_t i = 0; i < n; i++)
360 r[i] = octave::math::min (x, y[i]);
361}
362
363template <typename T>
364inline void
365mx_inline_xmax (std::size_t n, T *r, const T *x, const T *y)
366{
367 for (std::size_t i = 0; i < n; i++)
368 r[i] = octave::math::max (x[i], y[i]);
369}
370
371template <typename T>
372inline void
373mx_inline_xmax (std::size_t n, T *r, const T *x, T y)
374{
375 for (std::size_t i = 0; i < n; i++)
376 r[i] = octave::math::max (x[i], y);
377}
378
379template <typename T>
380inline void
381mx_inline_xmax (std::size_t n, T *r, T x, const T *y)
382{
383 for (std::size_t i = 0; i < n; i++)
384 r[i] = octave::math::max (x, y[i]);
385}
386
387// Specialize array-scalar max/min
388#define DEFMINMAXSPEC(T, F, OP) \
389 template <> \
390 inline void F<T> (std::size_t n, T *r, const T *x, T y) \
391 { \
392 if (octave::math::isnan (y)) \
393 std::memcpy (r, x, n * sizeof (T)); \
394 else \
395 for (std::size_t i = 0; i < n; i++) \
396 r[i] = (x[i] OP y ? x[i] : y); \
397 } \
398 template <> \
399 inline void F<T> (std::size_t n, T *r, T x, const T *y) \
400 { \
401 if (octave::math::isnan (x)) \
402 std::memcpy (r, y, n * sizeof (T)); \
403 else \
404 for (std::size_t i = 0; i < n; i++) \
405 r[i] = (y[i] OP x ? y[i] : x); \
406 }
407
408DEFMINMAXSPEC (double, mx_inline_xmin, <=)
409DEFMINMAXSPEC (double, mx_inline_xmax, >=)
410DEFMINMAXSPEC (float, mx_inline_xmin, <=)
411DEFMINMAXSPEC (float, mx_inline_xmax, >=)
412
413// FIXME: Is this comment correct anymore? It seems like std::pow is chosen.
414// Let the compiler decide which pow to use, whichever best matches the
415// arguments provided.
416
417template <typename R, typename X, typename Y>
418inline void
419mx_inline_pow (std::size_t n, R *r, const X *x, const Y *y)
420{
421 using std::pow;
422
423 for (std::size_t i = 0; i < n; i++)
424 r[i] = pow (x[i], y[i]);
425}
426
427template <typename R, typename X, typename Y>
428inline void
429mx_inline_pow (std::size_t n, R *r, const X *x, Y y)
430{
431 using std::pow;
432
433 for (std::size_t i = 0; i < n; i++)
434 r[i] = pow (x[i], y);
435}
436
437template <typename R, typename X, typename Y>
438inline void
439mx_inline_pow (std::size_t n, R *r, X x, const Y *y)
440{
441 using std::pow;
442
443 for (std::size_t i = 0; i < n; i++)
444 r[i] = pow (x, y[i]);
445}
446
447// Arbitrary function appliers.
448// The function is a template parameter to enable inlining.
449template <typename R, typename X, R fcn (X x)>
450inline void
451mx_inline_map (std::size_t n, R *r, const X *x)
452{
453 for (std::size_t i = 0; i < n; i++)
454 r[i] = fcn (x[i]);
455}
456
457template <typename R, typename X, R fcn (const X& x)>
458inline void
459mx_inline_map (std::size_t n, R *r, const X *x)
460{
461 for (std::size_t i = 0; i < n; i++)
462 r[i] = fcn (x[i]);
463}
464
465// Appliers. Since these call the operation just once, we pass it as
466// a pointer, to allow the compiler reduce number of instances.
467
468template <typename R, typename X>
469inline Array<R>
470do_mx_unary_op (const Array<X>& x,
471 void (*op) (std::size_t, R *, const X *))
472{
473 Array<R> r (x.dims ());
474 op (r.numel (), r.rwdata (), x.data ());
475 return r;
476}
477
478// Shortcuts for applying mx_inline_map.
479
480template <typename R, typename X, R fcn (X)>
481inline Array<R>
482do_mx_unary_map (const Array<X>& x)
483{
484 return do_mx_unary_op<R, X> (x, mx_inline_map<R, X, fcn>);
485}
486
487template <typename R, typename X, R fcn (const X&)>
488inline Array<R>
489do_mx_unary_map (const Array<X>& x)
490{
491 return do_mx_unary_op<R, X> (x, mx_inline_map<R, X, fcn>);
492}
493
494template <typename R>
495inline Array<R>&
496do_mx_inplace_op (Array<R>& r,
497 void (*op) (std::size_t, R *))
498{
499 op (r.numel (), r.rwdata ());
500 return r;
501}
502
503template <typename R, typename X, typename Y>
504inline Array<R>
505do_mm_binary_op (const Array<X>& x, const Array<Y>& y,
506 void (*op) (std::size_t, R *, const X *, const Y *),
507 void (*op1) (std::size_t, R *, X, const Y *),
508 void (*op2) (std::size_t, R *, const X *, Y),
509 const char *opname)
510{
511 const dim_vector& dx = x.dims ();
512 const dim_vector& dy = y.dims ();
513 if (dx == dy)
514 {
515 Array<R> r (dx);
516 op (r.numel (), r.rwdata (), x.data (), y.data ());
517 return r;
518 }
519 else if (is_valid_bsxfun (opname, dx, dy))
520 {
521 return do_bsxfun_op (x, y, op, op1, op2);
522 }
523 else
524 octave::err_nonconformant (opname, dx, dy);
525}
526
527template <typename R, typename X, typename Y>
528inline Array<R>
529do_ms_binary_op (const Array<X>& x, const Y& y,
530 void (*op) (std::size_t, R *, const X *, Y))
531{
532 Array<R> r (x.dims ());
533 op (r.numel (), r.rwdata (), x.data (), y);
534 return r;
535}
536
537template <typename R, typename X, typename Y>
538inline Array<R>
539do_sm_binary_op (const X& x, const Array<Y>& y,
540 void (*op) (std::size_t, R *, X, const Y *))
541{
542 Array<R> r (y.dims ());
543 op (r.numel (), r.rwdata (), x, y.data ());
544 return r;
545}
546
547template <typename R, typename X>
548inline Array<R>&
549do_mm_inplace_op (Array<R>& r, const Array<X>& x,
550 void (*op) (std::size_t, R *, const X *),
551 void (*op1) (std::size_t, R *, X),
552 const char *opname)
553{
554 const dim_vector& dr = r.dims ();
555 const dim_vector& dx = x.dims ();
556 if (dr == dx)
557 op (r.numel (), r.rwdata (), x.data ());
558 else if (is_valid_inplace_bsxfun (opname, dr, dx))
559 do_inplace_bsxfun_op (r, x, op, op1);
560 else
561 octave::err_nonconformant (opname, dr, dx);
562
563 return r;
564}
565
566template <typename R, typename X>
567inline Array<R>&
568do_ms_inplace_op (Array<R>& r, const X& x,
569 void (*op) (std::size_t, R *, X))
570{
571 op (r.numel (), r.rwdata (), x);
572 return r;
573}
574
575template <typename T1, typename T2>
576inline bool
577mx_inline_equal (std::size_t n, const T1 *x, const T2 *y)
578{
579 for (std::size_t i = 0; i < n; i++)
580 if (x[i] != y[i])
581 return false;
582 return true;
583}
584
585template <typename T>
586inline bool
587do_mx_check (const Array<T>& a,
588 bool (*op) (std::size_t, const T *))
589{
590 return op (a.numel (), a.data ());
591}
592
593// NOTE: we don't use std::norm because it typically does some heavyweight
594// magic to avoid underflows, which we don't need here.
595template <typename T>
596inline T
597cabsq (const std::complex<T>& c)
598{
599 return c.real () * c.real () + c.imag () * c.imag ();
600}
601
602// default. works for integers and bool.
603template <typename T>
604inline bool
605xis_true (T x)
606{
607 return x;
608}
609
610template <typename T>
611inline bool
612xis_false (T x)
613{
614 return ! x;
615}
616
617// for octave_ints
618template <typename T>
619inline bool
620xis_true (const octave_int<T>& x)
621{
622 return x.value ();
623}
624
625template <typename T>
626inline bool
627xis_false (const octave_int<T>& x)
628{
629 return ! x.value ();
630}
631
632// for reals, we want to ignore NaNs.
633inline bool
634xis_true (double x)
635{
636 return ! octave::math::isnan (x) && x != 0.0;
637}
638
639inline bool
640xis_false (double x)
641{
642 return x == 0.0;
643}
644
645inline bool
646xis_true (float x)
647{
648 return ! octave::math::isnan (x) && x != 0.0f;
649}
650
651inline bool
652xis_false (float x)
653{
654 return x == 0.0f;
655}
656
657// Ditto for complex.
658inline bool
659xis_true (const Complex& x)
660{
661 return ! octave::math::isnan (x) && x != 0.0;
662}
663
664inline bool
665xis_false (const Complex& x)
666{
667 return x == 0.0;
668}
669
670inline bool
671xis_true (const FloatComplex& x)
672{
673 return ! octave::math::isnan (x) && x != 0.0f;
674}
675
676inline bool
677xis_false (const FloatComplex& x)
678{
679 return x == 0.0f;
680}
681
682#define OP_RED_SUM(ac, el) ac += el
683#define OP_RED_PROD(ac, el) ac *= el
684#define OP_RED_SUMSQ(ac, el) ac += ((el)*(el))
685#define OP_RED_SUMSQC(ac, el) ac += cabsq (el)
686
687inline void
688op_dble_prod (double& ac, float el)
689{
690 ac *= el;
691}
692
693// FIXME: guaranteed?
694inline void
695op_dble_prod (Complex& ac, const FloatComplex& el)
696{
697 ac *= el;
698}
699
700template <typename T>
701inline void
702op_dble_prod (double& ac, const octave_int<T>& el)
703{
704 ac *= el.double_value ();
705}
706
707inline void
708op_dble_sum (double& ac, float el)
709{
710 ac += el;
711}
712
713// FIXME: guaranteed?
714inline void
715op_dble_sum (Complex& ac, const FloatComplex& el)
716{
717 ac += el;
718}
719
720template <typename T>
721inline void
722op_dble_sum (double& ac, const octave_int<T>& el)
723{
724 ac += el.double_value ();
725}
726
727// The following two implement a simple short-circuiting.
728#define OP_RED_ANYC(ac, el) \
729 if (xis_true (el)) \
730 { \
731 ac = true; \
732 break; \
733 } \
734 else \
735 continue
736
737#define OP_RED_ALLC(ac, el) \
738 if (xis_false (el)) \
739 { \
740 ac = false; \
741 break; \
742 } \
743 else \
744 continue
745
746#define OP_RED_FCN(F, TSRC, TRES, OP, ZERO) \
747 template <typename T> \
748 inline TRES \
749 F (const TSRC *v, octave_idx_type n) \
750 { \
751 TRES ac = ZERO; \
752 for (octave_idx_type i = 0; i < n; i++) \
753 OP(ac, v[i]); \
754 return ac; \
755 }
756
757#define PROMOTE_DOUBLE(T) \
758 typename subst_template_param<std::complex, T, double>::type
759
760OP_RED_FCN (mx_inline_sum, T, T, OP_RED_SUM, 0)
761OP_RED_FCN (mx_inline_dsum, T, PROMOTE_DOUBLE(T), op_dble_sum, 0.0)
762OP_RED_FCN (mx_inline_count, bool, T, OP_RED_SUM, 0)
763OP_RED_FCN (mx_inline_prod, T, T, OP_RED_PROD, 1)
764OP_RED_FCN (mx_inline_dprod, T, PROMOTE_DOUBLE(T), op_dble_prod, 1)
765OP_RED_FCN (mx_inline_sumsq, T, T, OP_RED_SUMSQ, 0)
766OP_RED_FCN (mx_inline_sumsq, std::complex<T>, T, OP_RED_SUMSQC, 0)
767OP_RED_FCN (mx_inline_any, T, bool, OP_RED_ANYC, false)
768OP_RED_FCN (mx_inline_all, T, bool, OP_RED_ALLC, true)
769
770#define OP_RED_FCN2(F, TSRC, TRES, OP, ZERO) \
771 template <typename T> \
772 inline void \
773 F (const TSRC *v, TRES *r, octave_idx_type m, octave_idx_type n) \
774 { \
775 for (octave_idx_type i = 0; i < m; i++) \
776 r[i] = ZERO; \
777 for (octave_idx_type j = 0; j < n; j++) \
778 { \
779 for (octave_idx_type i = 0; i < m; i++) \
780 OP(r[i], v[i]); \
781 v += m; \
782 } \
783 }
784
785OP_RED_FCN2 (mx_inline_sum, T, T, OP_RED_SUM, 0)
786OP_RED_FCN2 (mx_inline_dsum, T, PROMOTE_DOUBLE(T), op_dble_sum, 0.0)
787OP_RED_FCN2 (mx_inline_count, bool, T, OP_RED_SUM, 0)
788OP_RED_FCN2 (mx_inline_prod, T, T, OP_RED_PROD, 1)
789OP_RED_FCN2 (mx_inline_dprod, T, PROMOTE_DOUBLE(T), op_dble_prod, 0.0)
790OP_RED_FCN2 (mx_inline_sumsq, T, T, OP_RED_SUMSQ, 0)
791OP_RED_FCN2 (mx_inline_sumsq, std::complex<T>, T, OP_RED_SUMSQC, 0)
792
793#define OP_RED_ANYR(ac, el) ac |= xis_true (el)
794#define OP_RED_ALLR(ac, el) ac &= xis_true (el)
795
796OP_RED_FCN2 (mx_inline_any_r, T, bool, OP_RED_ANYR, false)
797OP_RED_FCN2 (mx_inline_all_r, T, bool, OP_RED_ALLR, true)
798
799// Using the general code for any/all would sacrifice short-circuiting.
800// OTOH, going by rows would sacrifice cache-coherence. The following
801// algorithm will achieve both, at the cost of a temporary octave_idx_type
802// array.
803
804#define OP_ROW_SHORT_CIRCUIT(F, PRED, ZERO) \
805 template <typename T> \
806 inline void \
807 F (const T *v, bool *r, octave_idx_type m, octave_idx_type n) \
808 { \
809 if (n <= 8) \
810 return F ## _r (v, r, m, n); \
811 \
812 /* FIXME: it may be sub-optimal to allocate the buffer here. */ \
813 OCTAVE_LOCAL_BUFFER (octave_idx_type, iact, m); \
814 for (octave_idx_type i = 0; i < m; i++) iact[i] = i; \
815 octave_idx_type nact = m; \
816 for (octave_idx_type j = 0; j < n; j++) \
817 { \
818 octave_idx_type k = 0; \
819 for (octave_idx_type i = 0; i < nact; i++) \
820 { \
821 octave_idx_type ia = iact[i]; \
822 if (! PRED (v[ia])) \
823 iact[k++] = ia; \
824 } \
825 nact = k; \
826 v += m; \
827 } \
828 for (octave_idx_type i = 0; i < m; i++) r[i] = ! ZERO; \
829 for (octave_idx_type i = 0; i < nact; i++) r[iact[i]] = ZERO; \
830 }
831
832OP_ROW_SHORT_CIRCUIT (mx_inline_any, xis_true, false)
833OP_ROW_SHORT_CIRCUIT (mx_inline_all, xis_false, true)
834
835#define OP_RED_FCNN(F, TSRC, TRES) \
836 template <typename T> \
837 inline void \
838 F (const TSRC *v, TRES *r, octave_idx_type l, \
839 octave_idx_type n, octave_idx_type u) \
840 { \
841 if (l == 1) \
842 { \
843 for (octave_idx_type i = 0; i < u; i++) \
844 { \
845 r[i] = F<T> (v, n); \
846 v += n; \
847 } \
848 } \
849 else \
850 { \
851 for (octave_idx_type i = 0; i < u; i++) \
852 { \
853 F (v, r, l, n); \
854 v += l*n; \
855 r += l; \
856 } \
857 } \
858 }
859
860OP_RED_FCNN (mx_inline_sum, T, T)
861OP_RED_FCNN (mx_inline_dsum, T, PROMOTE_DOUBLE(T))
862OP_RED_FCNN (mx_inline_count, bool, T)
863OP_RED_FCNN (mx_inline_prod, T, T)
864OP_RED_FCNN (mx_inline_dprod, T, PROMOTE_DOUBLE(T))
865OP_RED_FCNN (mx_inline_sumsq, T, T)
866OP_RED_FCNN (mx_inline_sumsq, std::complex<T>, T)
867OP_RED_FCNN (mx_inline_any, T, bool)
868OP_RED_FCNN (mx_inline_all, T, bool)
869
870#define OP_CUM_FCN(F, TSRC, TRES, OP) \
871 template <typename T> \
872 inline void \
873 F (const TSRC *v, TRES *r, octave_idx_type n) \
874 { \
875 if (n) \
876 { \
877 TRES t = r[0] = v[0]; \
878 for (octave_idx_type i = 1; i < n; i++) \
879 r[i] = t = t OP v[i]; \
880 } \
881 }
882
883OP_CUM_FCN (mx_inline_cumsum, T, T, +)
884OP_CUM_FCN (mx_inline_cumprod, T, T, *)
885OP_CUM_FCN (mx_inline_cumcount, bool, T, +)
886
887#define OP_CUM_FCN2(F, TSRC, TRES, OP) \
888 template <typename T> \
889 inline void \
890 F (const TSRC *v, TRES *r, octave_idx_type m, octave_idx_type n) \
891 { \
892 if (n) \
893 { \
894 for (octave_idx_type i = 0; i < m; i++) \
895 r[i] = v[i]; \
896 const T *r0 = r; \
897 for (octave_idx_type j = 1; j < n; j++) \
898 { \
899 r += m; v += m; \
900 for (octave_idx_type i = 0; i < m; i++) \
901 r[i] = r0[i] OP v[i]; \
902 r0 += m; \
903 } \
904 } \
905 }
906
907OP_CUM_FCN2 (mx_inline_cumsum, T, T, +)
908OP_CUM_FCN2 (mx_inline_cumprod, T, T, *)
909OP_CUM_FCN2 (mx_inline_cumcount, bool, T, +)
910
911#define OP_CUM_FCNN(F, TSRC, TRES) \
912 template <typename T> \
913 inline void \
914 F (const TSRC *v, TRES *r, octave_idx_type l, \
915 octave_idx_type n, octave_idx_type u) \
916 { \
917 if (l == 1) \
918 { \
919 for (octave_idx_type i = 0; i < u; i++) \
920 { \
921 F (v, r, n); \
922 v += n; \
923 r += n; \
924 } \
925 } \
926 else \
927 { \
928 for (octave_idx_type i = 0; i < u; i++) \
929 { \
930 F (v, r, l, n); \
931 v += l*n; \
932 r += l*n; \
933 } \
934 } \
935 }
936
937OP_CUM_FCNN (mx_inline_cumsum, T, T)
938OP_CUM_FCNN (mx_inline_cumprod, T, T)
939OP_CUM_FCNN (mx_inline_cumcount, bool, T)
940
941#define OP_MINMAX_FCN(F, OP) \
942 template <typename T> \
943 void F (const T *v, T *r, octave_idx_type n) \
944 { \
945 if (! n) \
946 return; \
947 T tmp = v[0]; \
948 octave_idx_type i = 1; \
949 if (octave::math::isnan (tmp)) \
950 { \
951 for (; i < n && octave::math::isnan (v[i]); i++) ; \
952 if (i < n) \
953 tmp = v[i]; \
954 } \
955 for (; i < n; i++) \
956 if (v[i] OP tmp) \
957 tmp = v[i]; \
958 *r = tmp; \
959 } \
960 template <typename T> \
961 void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n) \
962 { \
963 if (! n) \
964 return; \
965 T tmp = v[0]; \
966 octave_idx_type tmpi = 0; \
967 octave_idx_type i = 1; \
968 if (octave::math::isnan (tmp)) \
969 { \
970 for (; i < n && octave::math::isnan (v[i]); i++) ; \
971 if (i < n) \
972 { \
973 tmp = v[i]; \
974 tmpi = i; \
975 } \
976 } \
977 for (; i < n; i++) \
978 if (v[i] OP tmp) \
979 { \
980 tmp = v[i]; \
981 tmpi = i; \
982 } \
983 *r = tmp; \
984 *ri = tmpi; \
985 }
986
987OP_MINMAX_FCN (mx_inline_min, <)
988OP_MINMAX_FCN (mx_inline_max, >)
989
990// Row reductions will be slightly complicated. We will proceed with checks
991// for NaNs until we detect that no row will yield a NaN, in which case we
992// proceed to a faster code.
993
994#define OP_MINMAX_FCN2(F, OP) \
995 template <typename T> \
996 inline void \
997 F (const T *v, T *r, octave_idx_type m, octave_idx_type n) \
998 { \
999 if (! n) \
1000 return; \
1001 bool nan = false; \
1002 octave_idx_type j = 0; \
1003 for (octave_idx_type i = 0; i < m; i++) \
1004 { \
1005 r[i] = v[i]; \
1006 if (octave::math::isnan (v[i])) \
1007 nan = true; \
1008 } \
1009 j++; \
1010 v += m; \
1011 while (nan && j < n) \
1012 { \
1013 nan = false; \
1014 for (octave_idx_type i = 0; i < m; i++) \
1015 { \
1016 if (octave::math::isnan (v[i])) \
1017 nan = true; \
1018 else if (octave::math::isnan (r[i]) || v[i] OP r[i]) \
1019 r[i] = v[i]; \
1020 } \
1021 j++; \
1022 v += m; \
1023 } \
1024 while (j < n) \
1025 { \
1026 for (octave_idx_type i = 0; i < m; i++) \
1027 if (v[i] OP r[i]) \
1028 r[i] = v[i]; \
1029 j++; \
1030 v += m; \
1031 } \
1032 } \
1033 template <typename T> \
1034 inline void \
1035 F (const T *v, T *r, octave_idx_type *ri, \
1036 octave_idx_type m, octave_idx_type n) \
1037 { \
1038 if (! n) \
1039 return; \
1040 bool nan = false; \
1041 octave_idx_type j = 0; \
1042 for (octave_idx_type i = 0; i < m; i++) \
1043 { \
1044 r[i] = v[i]; \
1045 ri[i] = j; \
1046 if (octave::math::isnan (v[i])) \
1047 nan = true; \
1048 } \
1049 j++; \
1050 v += m; \
1051 while (nan && j < n) \
1052 { \
1053 nan = false; \
1054 for (octave_idx_type i = 0; i < m; i++) \
1055 { \
1056 if (octave::math::isnan (v[i])) \
1057 nan = true; \
1058 else if (octave::math::isnan (r[i]) || v[i] OP r[i]) \
1059 { \
1060 r[i] = v[i]; \
1061 ri[i] = j; \
1062 } \
1063 } \
1064 j++; \
1065 v += m; \
1066 } \
1067 while (j < n) \
1068 { \
1069 for (octave_idx_type i = 0; i < m; i++) \
1070 if (v[i] OP r[i]) \
1071 { \
1072 r[i] = v[i]; \
1073 ri[i] = j; \
1074 } \
1075 j++; \
1076 v += m; \
1077 } \
1078 }
1079
1080OP_MINMAX_FCN2 (mx_inline_min, <)
1081OP_MINMAX_FCN2 (mx_inline_max, >)
1082
1083#define OP_MINMAX_FCNN(F) \
1084 template <typename T> \
1085 inline void \
1086 F (const T *v, T *r, octave_idx_type l, \
1087 octave_idx_type n, octave_idx_type u) \
1088 { \
1089 if (! n) \
1090 return; \
1091 if (l == 1) \
1092 { \
1093 for (octave_idx_type i = 0; i < u; i++) \
1094 { \
1095 F (v, r, n); \
1096 v += n; \
1097 r++; \
1098 } \
1099 } \
1100 else \
1101 { \
1102 for (octave_idx_type i = 0; i < u; i++) \
1103 { \
1104 F (v, r, l, n); \
1105 v += l*n; \
1106 r += l; \
1107 } \
1108 } \
1109 } \
1110 template <typename T> \
1111 inline void \
1112 F (const T *v, T *r, octave_idx_type *ri, \
1113 octave_idx_type l, octave_idx_type n, octave_idx_type u) \
1114 { \
1115 if (! n) return; \
1116 if (l == 1) \
1117 { \
1118 for (octave_idx_type i = 0; i < u; i++) \
1119 { \
1120 F (v, r, ri, n); \
1121 v += n; \
1122 r++; \
1123 ri++; \
1124 } \
1125 } \
1126 else \
1127 { \
1128 for (octave_idx_type i = 0; i < u; i++) \
1129 { \
1130 F (v, r, ri, l, n); \
1131 v += l*n; \
1132 r += l; \
1133 ri += l; \
1134 } \
1135 } \
1136 }
1137
1138OP_MINMAX_FCNN (mx_inline_min)
1139OP_MINMAX_FCNN (mx_inline_max)
1140
1141#define OP_CUMMINMAX_FCN(F, OP) \
1142 template <typename T> \
1143 void F (const T *v, T *r, octave_idx_type n) \
1144 { \
1145 if (! n) \
1146 return; \
1147 T tmp = v[0]; \
1148 octave_idx_type i = 1; \
1149 octave_idx_type j = 0; \
1150 if (octave::math::isnan (tmp)) \
1151 { \
1152 for (; i < n && octave::math::isnan (v[i]); i++) ; \
1153 for (; j < i; j++) \
1154 r[j] = tmp; \
1155 if (i < n) \
1156 tmp = v[i]; \
1157 } \
1158 for (; i < n; i++) \
1159 if (v[i] OP tmp) \
1160 { \
1161 for (; j < i; j++) \
1162 r[j] = tmp; \
1163 tmp = v[i]; \
1164 } \
1165 for (; j < i; j++) \
1166 r[j] = tmp; \
1167 } \
1168 template <typename T> \
1169 void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n) \
1170 { \
1171 if (! n) \
1172 return; \
1173 T tmp = v[0]; \
1174 octave_idx_type tmpi = 0; \
1175 octave_idx_type i = 1; \
1176 octave_idx_type j = 0; \
1177 if (octave::math::isnan (tmp)) \
1178 { \
1179 for (; i < n && octave::math::isnan (v[i]); i++) ; \
1180 for (; j < i; j++) \
1181 { \
1182 r[j] = tmp; \
1183 ri[j] = tmpi; \
1184 } \
1185 if (i < n) \
1186 { \
1187 tmp = v[i]; \
1188 tmpi = i; \
1189 } \
1190 } \
1191 for (; i < n; i++) \
1192 if (v[i] OP tmp) \
1193 { \
1194 for (; j < i; j++) \
1195 { \
1196 r[j] = tmp; \
1197 ri[j] = tmpi; \
1198 } \
1199 tmp = v[i]; \
1200 tmpi = i; \
1201 } \
1202 for (; j < i; j++) \
1203 { \
1204 r[j] = tmp; \
1205 ri[j] = tmpi; \
1206 } \
1207 }
1208
1209OP_CUMMINMAX_FCN (mx_inline_cummin, <)
1210OP_CUMMINMAX_FCN (mx_inline_cummax, >)
1211
1212// Row reductions will be slightly complicated. We will proceed with checks
1213// for NaNs until we detect that no row will yield a NaN, in which case we
1214// proceed to a faster code.
1215
1216#define OP_CUMMINMAX_FCN2(F, OP) \
1217 template <typename T> \
1218 inline void \
1219 F (const T *v, T *r, octave_idx_type m, octave_idx_type n) \
1220 { \
1221 if (! n) \
1222 return; \
1223 bool nan = false; \
1224 const T *r0; \
1225 octave_idx_type j = 0; \
1226 for (octave_idx_type i = 0; i < m; i++) \
1227 { \
1228 r[i] = v[i]; \
1229 if (octave::math::isnan (v[i])) \
1230 nan = true; \
1231 } \
1232 j++; \
1233 v += m; \
1234 r0 = r; \
1235 r += m; \
1236 while (nan && j < n) \
1237 { \
1238 nan = false; \
1239 for (octave_idx_type i = 0; i < m; i++) \
1240 { \
1241 if (octave::math::isnan (v[i])) \
1242 { \
1243 r[i] = r0[i]; \
1244 nan = true; \
1245 } \
1246 else if (octave::math::isnan (r0[i]) || v[i] OP r0[i]) \
1247 r[i] = v[i]; \
1248 else \
1249 r[i] = r0[i]; \
1250 } \
1251 j++; \
1252 v += m; \
1253 r0 = r; \
1254 r += m; \
1255 } \
1256 while (j < n) \
1257 { \
1258 for (octave_idx_type i = 0; i < m; i++) \
1259 if (v[i] OP r0[i]) \
1260 r[i] = v[i]; \
1261 else \
1262 r[i] = r0[i]; \
1263 j++; \
1264 v += m; \
1265 r0 = r; \
1266 r += m; \
1267 } \
1268 } \
1269 template <typename T> \
1270 inline void \
1271 F (const T *v, T *r, octave_idx_type *ri, \
1272 octave_idx_type m, octave_idx_type n) \
1273 { \
1274 if (! n) \
1275 return; \
1276 bool nan = false; \
1277 const T *r0; \
1278 const octave_idx_type *r0i; \
1279 octave_idx_type j = 0; \
1280 for (octave_idx_type i = 0; i < m; i++) \
1281 { \
1282 r[i] = v[i]; ri[i] = 0; \
1283 if (octave::math::isnan (v[i])) \
1284 nan = true; \
1285 } \
1286 j++; \
1287 v += m; \
1288 r0 = r; \
1289 r += m; \
1290 r0i = ri; \
1291 ri += m; \
1292 while (nan && j < n) \
1293 { \
1294 nan = false; \
1295 for (octave_idx_type i = 0; i < m; i++) \
1296 { \
1297 if (octave::math::isnan (v[i])) \
1298 { \
1299 r[i] = r0[i]; \
1300 ri[i] = r0i[i]; \
1301 nan = true; \
1302 } \
1303 else if (octave::math::isnan (r0[i]) || v[i] OP r0[i]) \
1304 { \
1305 r[i] = v[i]; \
1306 ri[i] = j; \
1307 } \
1308 else \
1309 { \
1310 r[i] = r0[i]; \
1311 ri[i] = r0i[i]; \
1312 } \
1313 } \
1314 j++; \
1315 v += m; \
1316 r0 = r; \
1317 r += m; \
1318 r0i = ri; \
1319 ri += m; \
1320 } \
1321 while (j < n) \
1322 { \
1323 for (octave_idx_type i = 0; i < m; i++) \
1324 if (v[i] OP r0[i]) \
1325 { \
1326 r[i] = v[i]; \
1327 ri[i] = j; \
1328 } \
1329 else \
1330 { \
1331 r[i] = r0[i]; \
1332 ri[i] = r0i[i]; \
1333 } \
1334 j++; \
1335 v += m; \
1336 r0 = r; \
1337 r += m; \
1338 r0i = ri; \
1339 ri += m; \
1340 } \
1341 }
1342
1343OP_CUMMINMAX_FCN2 (mx_inline_cummin, <)
1344OP_CUMMINMAX_FCN2 (mx_inline_cummax, >)
1345
1346#define OP_CUMMINMAX_FCNN(F) \
1347 template <typename T> \
1348 inline void \
1349 F (const T *v, T *r, octave_idx_type l, \
1350 octave_idx_type n, octave_idx_type u) \
1351 { \
1352 if (! n) \
1353 return; \
1354 if (l == 1) \
1355 { \
1356 for (octave_idx_type i = 0; i < u; i++) \
1357 { \
1358 F (v, r, n); \
1359 v += n; \
1360 r += n; \
1361 } \
1362 } \
1363 else \
1364 { \
1365 for (octave_idx_type i = 0; i < u; i++) \
1366 { \
1367 F (v, r, l, n); \
1368 v += l*n; \
1369 r += l*n; \
1370 } \
1371 } \
1372 } \
1373 template <typename T> \
1374 inline void \
1375 F (const T *v, T *r, octave_idx_type *ri, \
1376 octave_idx_type l, octave_idx_type n, octave_idx_type u) \
1377 { \
1378 if (! n) \
1379 return; \
1380 if (l == 1) \
1381 { \
1382 for (octave_idx_type i = 0; i < u; i++) \
1383 { \
1384 F (v, r, ri, n); \
1385 v += n; \
1386 r += n; \
1387 ri += n; \
1388 } \
1389 } \
1390 else \
1391 { \
1392 for (octave_idx_type i = 0; i < u; i++) \
1393 { \
1394 F (v, r, ri, l, n); \
1395 v += l*n; \
1396 r += l*n; \
1397 ri += l*n; \
1398 } \
1399 } \
1400 }
1401
1402OP_CUMMINMAX_FCNN (mx_inline_cummin)
1403OP_CUMMINMAX_FCNN (mx_inline_cummax)
1404
1405template <typename T>
1406void mx_inline_diff (const T *v, T *r, octave_idx_type n,
1407 octave_idx_type order)
1408{
1409 switch (order)
1410 {
1411 case 1:
1412 for (octave_idx_type i = 0; i < n-1; i++)
1413 r[i] = v[i+1] - v[i];
1414 break;
1415 case 2:
1416 if (n > 1)
1417 {
1418 T lst = v[1] - v[0];
1419 for (octave_idx_type i = 0; i < n-2; i++)
1420 {
1421 T dif = v[i+2] - v[i+1];
1422 r[i] = dif - lst;
1423 lst = dif;
1424 }
1425 }
1426 break;
1427 default:
1428 {
1429 OCTAVE_LOCAL_BUFFER (T, buf, n-1);
1430
1431 for (octave_idx_type i = 0; i < n-1; i++)
1432 buf[i] = v[i+1] - v[i];
1433
1434 for (octave_idx_type o = 2; o <= order; o++)
1435 {
1436 for (octave_idx_type i = 0; i < n-o; i++)
1437 buf[i] = buf[i+1] - buf[i];
1438 }
1439
1440 for (octave_idx_type i = 0; i < n-order; i++)
1441 r[i] = buf[i];
1442 }
1443 }
1444}
1445
1446template <typename T>
1447void
1448mx_inline_diff (const T *v, T *r,
1449 octave_idx_type m, octave_idx_type n,
1450 octave_idx_type order)
1451{
1452 switch (order)
1453 {
1454 case 1:
1455 for (octave_idx_type i = 0; i < m*(n-1); i++)
1456 r[i] = v[i+m] - v[i];
1457 break;
1458 case 2:
1459 for (octave_idx_type i = 0; i < n-2; i++)
1460 {
1461 for (octave_idx_type j = i*m; j < i*m+m; j++)
1462 r[j] = (v[j+m+m] - v[j+m]) - (v[j+m] - v[j]);
1463 }
1464 break;
1465 default:
1466 {
1467 OCTAVE_LOCAL_BUFFER (T, buf, n-1);
1468
1469 for (octave_idx_type j = 0; j < m; j++)
1470 {
1471 for (octave_idx_type i = 0; i < n-1; i++)
1472 buf[i] = v[i*m+j+m] - v[i*m+j];
1473
1474 for (octave_idx_type o = 2; o <= order; o++)
1475 {
1476 for (octave_idx_type i = 0; i < n-o; i++)
1477 buf[i] = buf[i+1] - buf[i];
1478 }
1479
1480 for (octave_idx_type i = 0; i < n-order; i++)
1481 r[i*m+j] = buf[i];
1482 }
1483 }
1484 }
1485}
1486
1487template <typename T>
1488inline void
1489mx_inline_diff (const T *v, T *r,
1490 octave_idx_type l, octave_idx_type n, octave_idx_type u,
1491 octave_idx_type order)
1492{
1493 if (! n) return;
1494 if (l == 1)
1495 {
1496 for (octave_idx_type i = 0; i < u; i++)
1497 {
1498 mx_inline_diff (v, r, n, order);
1499 v += n; r += n-order;
1500 }
1501 }
1502 else
1503 {
1504 for (octave_idx_type i = 0; i < u; i++)
1505 {
1506 mx_inline_diff (v, r, l, n, order);
1507 v += l*n;
1508 r += l*(n-order);
1509 }
1510 }
1511}
1512
1513// Assistant function
1514
1515inline void
1516get_extent_triplet (const dim_vector& dims, int& dim,
1517 octave_idx_type& l, octave_idx_type& n,
1518 octave_idx_type& u)
1519{
1520 octave_idx_type ndims = dims.ndims ();
1521 if (dim >= ndims)
1522 {
1523 l = dims.numel ();
1524 n = 1;
1525 u = 1;
1526 }
1527 else
1528 {
1529 if (dim < 0)
1530 dim = dims.first_non_singleton ();
1531
1532 // calculate extent triplet.
1533 l = 1, n = dims(dim), u = 1;
1534 for (octave_idx_type i = 0; i < dim; i++)
1535 l *= dims(i);
1536 for (octave_idx_type i = dim + 1; i < ndims; i++)
1537 u *= dims(i);
1538 }
1539}
1540
1541// Appliers.
1542// FIXME: is this the best design? C++ gives a lot of options here...
1543// maybe it can be done without an explicit parameter?
1544
1545template <typename R, typename T>
1546inline Array<R>
1547do_mx_red_op (const Array<T>& src, int dim,
1548 void (*mx_red_op) (const T *, R *, octave_idx_type,
1549 octave_idx_type, octave_idx_type))
1550{
1551 octave_idx_type l, n, u;
1552 dim_vector dims = src.dims ();
1553 // M*b inconsistency: sum ([]) = 0 etc.
1554 if (dims.ndims () == 2 && dims(0) == 0 && dims(1) == 0)
1555 dims(1) = 1;
1556
1557 get_extent_triplet (dims, dim, l, n, u);
1558
1559 // Reduction operation reduces the array size.
1560 if (dim < dims.ndims ()) dims(dim) = 1;
1561 dims.chop_trailing_singletons ();
1562
1563 Array<R> ret (dims);
1564 mx_red_op (src.data (), ret.rwdata (), l, n, u);
1565
1566 return ret;
1567}
1568
1569template <typename R, typename T>
1570inline Array<R>
1571do_mx_cum_op (const Array<T>& src, int dim,
1572 void (*mx_cum_op) (const T *, R *, octave_idx_type,
1573 octave_idx_type, octave_idx_type))
1574{
1575 octave_idx_type l, n, u;
1576 const dim_vector& dims = src.dims ();
1577 get_extent_triplet (dims, dim, l, n, u);
1578
1579 // Cumulative operation doesn't reduce the array size.
1580 Array<R> ret (dims);
1581 mx_cum_op (src.data (), ret.rwdata (), l, n, u);
1582
1583 return ret;
1584}
1585
1586template <typename R>
1587inline Array<R>
1588do_mx_minmax_op (const Array<R>& src, int dim,
1589 void (*mx_minmax_op) (const R *, R *, octave_idx_type,
1590 octave_idx_type, octave_idx_type))
1591{
1592 octave_idx_type l, n, u;
1593 dim_vector dims = src.dims ();
1594 get_extent_triplet (dims, dim, l, n, u);
1595
1596 // If the dimension is zero, we don't do anything.
1597 if (dim < dims.ndims () && dims(dim) != 0) dims(dim) = 1;
1598 dims.chop_trailing_singletons ();
1599
1600 Array<R> ret (dims);
1601 mx_minmax_op (src.data (), ret.rwdata (), l, n, u);
1602
1603 return ret;
1604}
1605
1606template <typename R>
1607inline Array<R>
1608do_mx_minmax_op (const Array<R>& src, Array<octave_idx_type>& idx, int dim,
1609 void (*mx_minmax_op) (const R *, R *, octave_idx_type *,
1610 octave_idx_type, octave_idx_type, octave_idx_type))
1611{
1612 octave_idx_type l, n, u;
1613 dim_vector dims = src.dims ();
1614 get_extent_triplet (dims, dim, l, n, u);
1615
1616 // If the dimension is zero, we don't do anything.
1617 if (dim < dims.ndims () && dims(dim) != 0) dims(dim) = 1;
1618 dims.chop_trailing_singletons ();
1619
1620 Array<R> ret (dims);
1621 if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
1622
1623 mx_minmax_op (src.data (), ret.rwdata (), idx.rwdata (),
1624 l, n, u);
1625
1626 return ret;
1627}
1628
1629template <typename R>
1630inline Array<R>
1631do_mx_cumminmax_op (const Array<R>& src, int dim,
1632 void (*mx_cumminmax_op) (const R *, R *, octave_idx_type,
1633 octave_idx_type, octave_idx_type))
1634{
1635 octave_idx_type l, n, u;
1636 const dim_vector& dims = src.dims ();
1637 get_extent_triplet (dims, dim, l, n, u);
1638
1639 Array<R> ret (dims);
1640 mx_cumminmax_op (src.data (), ret.rwdata (), l, n, u);
1641
1642 return ret;
1643}
1644
1645template <typename R>
1646inline Array<R>
1647do_mx_cumminmax_op (const Array<R>& src, Array<octave_idx_type>& idx, int dim,
1648 void (*mx_cumminmax_op) (const R *, R *, octave_idx_type *,
1649 octave_idx_type, octave_idx_type, octave_idx_type))
1650{
1651 octave_idx_type l, n, u;
1652 const dim_vector& dims = src.dims ();
1653 get_extent_triplet (dims, dim, l, n, u);
1654
1655 Array<R> ret (dims);
1656 if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
1657
1658 mx_cumminmax_op (src.data (), ret.rwdata (), idx.rwdata (),
1659 l, n, u);
1660
1661 return ret;
1662}
1663
1664template <typename R>
1665inline Array<R>
1666do_mx_diff_op (const Array<R>& src, int dim, octave_idx_type order,
1667 void (*mx_diff_op) (const R *, R *,
1668 octave_idx_type, octave_idx_type,
1669 octave_idx_type, octave_idx_type))
1670{
1671 octave_idx_type l, n, u;
1672 if (order <= 0)
1673 return src;
1674
1675 dim_vector dims = src.dims ();
1676
1677 get_extent_triplet (dims, dim, l, n, u);
1678 if (dim >= dims.ndims ())
1679 dims.resize (dim+1, 1);
1680
1681 if (dims(dim) <= order)
1682 {
1683 dims(dim) = 0;
1684 return Array<R> (dims);
1685 }
1686 else
1687 {
1688 dims(dim) -= order;
1689 }
1690
1691 Array<R> ret (dims);
1692 mx_diff_op (src.data (), ret.rwdata (), l, n, u, order);
1693
1694 return ret;
1695}
1696
1697// Fast extra-precise summation. According to
1698// T. Ogita, S. M. Rump, S. Oishi:
1699// Accurate Sum And Dot Product,
1700// SIAM J. Sci. Computing, Vol. 26, 2005
1701
1702template <typename T>
1703inline void
1704twosum_accum (T& s, T& e,
1705 const T& x)
1706{
1707 T s1 = s + x;
1708 T t = s1 - s;
1709 T e1 = (s - (s1 - t)) + (x - t);
1710 s = s1;
1711 e += e1;
1712}
1713
1714template <typename T>
1715inline T
1716mx_inline_xsum (const T *v, octave_idx_type n)
1717{
1718 T s, e;
1719 s = e = 0;
1720 for (octave_idx_type i = 0; i < n; i++)
1721 twosum_accum (s, e, v[i]);
1722
1723 return s + e;
1724}
1725
1726template <typename T>
1727inline void
1728mx_inline_xsum (const T *v, T *r,
1729 octave_idx_type m, octave_idx_type n)
1730{
1731 OCTAVE_LOCAL_BUFFER (T, e, m);
1732 for (octave_idx_type i = 0; i < m; i++)
1733 e[i] = r[i] = T ();
1734
1735 for (octave_idx_type j = 0; j < n; j++)
1736 {
1737 for (octave_idx_type i = 0; i < m; i++)
1738 twosum_accum (r[i], e[i], v[i]);
1739
1740 v += m;
1741 }
1742
1743 for (octave_idx_type i = 0; i < m; i++)
1744 r[i] += e[i];
1745}
1746
1747OP_RED_FCNN (mx_inline_xsum, T, T)
1748
1749#endif
do_bsxfun_op
Array< R > do_bsxfun_op(const Array< X > &x, const Array< Y > &y, void(*op_vv)(std::size_t, R *, const X *, const Y *), void(*op_sv)(std::size_t, R *, X, const Y *), void(*op_vs)(std::size_t, R *, const X *, Y))
Definition bsxfun-defs.cc:41
do_inplace_bsxfun_op
void do_inplace_bsxfun_op(Array< R > &r, const Array< X > &x, void(*op_vv)(std::size_t, R *, const X *), void(*op_vs)(std::size_t, R *, X))
Definition bsxfun-defs.cc:142
is_valid_inplace_bsxfun
bool is_valid_inplace_bsxfun(const std::string &name, const dim_vector &rdv, const dim_vector &xdv)
Definition bsxfun.h:63
is_valid_bsxfun
bool is_valid_bsxfun(const std::string &name, const dim_vector &xdv, const dim_vector &ydv)
Definition bsxfun.h:39
Array
N Dimensional Array with copy-on-write semantics.
Definition Array.h:130
Array::dims
const dim_vector & dims() const
Return a const-reference so that dims ()(i) works efficiently.
Definition Array.h:507
Array::data
const T * data() const
Size of the specified dimension.
Definition Array.h:665
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array.h:418
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:90
dim_vector::numel
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
Definition dim-vector.h:331
dim_vector::chop_trailing_singletons
void chop_trailing_singletons()
Definition dim-vector.h:160
dim_vector::resize
void resize(int n, int fill_value=0)
Definition dim-vector.h:268
dim_vector::ndims
octave_idx_type ndims() const
Number of dimensions.
Definition dim-vector.h:253
dim_vector::first_non_singleton
int first_non_singleton(int def=0) const
Definition dim-vector.h:440
octave_idx_type
octave_int
Definition oct-inttypes.h:783
octave_int::double_value
double double_value() const
Definition oct-inttypes.h:843
real
ColumnVector real(const ComplexColumnVector &a)
Definition dColVector.cc:137
imag
ColumnVector imag(const ComplexColumnVector &a)
Definition dColVector.cc:143
x
F77_RET_T const F77_DBLE * x
Definition lo-slatec-proto.h:38
mx_inline_map
void mx_inline_map(std::size_t n, R *r, const X *x)
Definition mx-inlines.cc:451
do_mx_unary_op
Array< R > do_mx_unary_op(const Array< X > &x, void(*op)(std::size_t, R *, const X *))
Definition mx-inlines.cc:470
mx_inline_any
bool mx_inline_any(const T *v, octave_idx_type n)
Definition mx-inlines.cc:767
mx_inline_div2
void mx_inline_div2(std::size_t n, R *r, const X *x)
Definition mx-inlines.cc:130
do_mx_check
bool do_mx_check(const Array< T > &a, bool(*op)(std::size_t, const T *))
Definition mx-inlines.cc:587
mx_inline_sub2
void mx_inline_sub2(std::size_t n, R *r, const X *x)
Definition mx-inlines.cc:128
mx_inline_notzero
void mx_inline_notzero(std::size_t n, bool *r, const X *x)
Definition mx-inlines.cc:81
DEFMXCMPOP
#define DEFMXCMPOP(F, OP)
Definition mx-inlines.cc:132
mx_inline_xmin
void mx_inline_xmin(std::size_t n, T *r, const T *x, const T *y)
Definition mx-inlines.cc:341
mx_inline_not_and
void mx_inline_not_and(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:221
OP_RED_ANYR
#define OP_RED_ANYR(ac, el)
Definition mx-inlines.cc:793
xis_false
bool xis_false(T x)
Definition mx-inlines.cc:612
do_mm_inplace_op
Array< R > & do_mm_inplace_op(Array< R > &r, const Array< X > &x, void(*op)(std::size_t, R *, const X *), void(*op1)(std::size_t, R *, X), const char *opname)
Definition mx-inlines.cc:549
OP_CUM_FCN2
#define OP_CUM_FCN2(F, TSRC, TRES, OP)
Definition mx-inlines.cc:887
do_mx_cumminmax_op
Array< R > do_mx_cumminmax_op(const Array< R > &src, int dim, void(*mx_cumminmax_op)(const R *, R *, octave_idx_type, octave_idx_type, octave_idx_type))
Definition mx-inlines.cc:1631
mx_inline_le
void mx_inline_le(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:153
OP_ROW_SHORT_CIRCUIT
#define OP_ROW_SHORT_CIRCUIT(F, PRED, ZERO)
Definition mx-inlines.cc:804
mx_inline_xsum
T mx_inline_xsum(const T *v, octave_idx_type n)
Definition mx-inlines.cc:1716
OP_RED_SUMSQC
#define OP_RED_SUMSQC(ac, el)
Definition mx-inlines.cc:685
mx_inline_add2
void mx_inline_add2(std::size_t n, R *r, const X *x)
Definition mx-inlines.cc:127
do_ms_binary_op
Array< R > do_ms_binary_op(const Array< X > &x, const Y &y, void(*op)(std::size_t, R *, const X *, Y))
Definition mx-inlines.cc:529
mx_inline_and_not
void mx_inline_and_not(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:223
OP_RED_FCNN
#define OP_RED_FCNN(F, TSRC, TRES)
Definition mx-inlines.cc:835
OP_CUMMINMAX_FCN
#define OP_CUMMINMAX_FCN(F, OP)
Definition mx-inlines.cc:1141
mx_inline_sub
void mx_inline_sub(std::size_t n, R *r, const X *x, const Y *y)
Definition mx-inlines.cc:109
mx_inline_uminus
void mx_inline_uminus(std::size_t n, R *r, const X *x)
Definition mx-inlines.cc:56
do_mx_inplace_op
Array< R > & do_mx_inplace_op(Array< R > &r, void(*op)(std::size_t, R *))
Definition mx-inlines.cc:496
mx_inline_dprod
subst_template_param< std::complex, T, double >::type mx_inline_dprod(const T *v, octave_idx_type n)
Definition mx-inlines.cc:764
mx_inline_uminus2
void mx_inline_uminus2(std::size_t n, R *r)
Definition mx-inlines.cc:64
mx_inline_eq
void mx_inline_eq(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:156
OP_RED_ALLC
#define OP_RED_ALLC(ac, el)
Definition mx-inlines.cc:737
mx_inline_cummin
void mx_inline_cummin(const T *v, T *r, octave_idx_type n)
Definition mx-inlines.cc:1209
DEFMXBOOLOP
#define DEFMXBOOLOP(F, NOT1, OP, NOT2)
Definition mx-inlines.cc:196
mx_inline_ge
void mx_inline_ge(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:155
mx_inline_all_real
bool mx_inline_all_real(std::size_t n, const std::complex< T > *x)
Definition mx-inlines.cc:312
PROMOTE_DOUBLE
#define PROMOTE_DOUBLE(T)
Definition mx-inlines.cc:757
mx_inline_or
void mx_inline_or(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:220
mx_inline_div
void mx_inline_div(std::size_t n, R *r, const X *x, const Y *y)
Definition mx-inlines.cc:111
OP_MINMAX_FCN2
#define OP_MINMAX_FCN2(F, OP)
Definition mx-inlines.cc:994
mx_inline_cumprod
void mx_inline_cumprod(const T *v, T *r, octave_idx_type n)
Definition mx-inlines.cc:884
do_mx_cum_op
Array< R > do_mx_cum_op(const Array< T > &src, int dim, void(*mx_cum_op)(const T *, R *, octave_idx_type, octave_idx_type, octave_idx_type))
Definition mx-inlines.cc:1571
mx_inline_cumsum
void mx_inline_cumsum(const T *v, T *r, octave_idx_type n)
Definition mx-inlines.cc:883
twosum_accum
void twosum_accum(T &s, T &e, const T &x)
Definition mx-inlines.cc:1704
mx_inline_any_nan
bool mx_inline_any_nan(std::size_t n, const T *x)
Definition mx-inlines.cc:260
mx_inline_and2
void mx_inline_and2(std::size_t n, bool *r, const X *x)
Definition mx-inlines.cc:228
do_mx_red_op
Array< R > do_mx_red_op(const Array< T > &src, int dim, void(*mx_red_op)(const T *, R *, octave_idx_type, octave_idx_type, octave_idx_type))
Definition mx-inlines.cc:1547
DEFMXBINOPEQ
#define DEFMXBINOPEQ(F, OP)
Definition mx-inlines.cc:113
mx_inline_max
void mx_inline_max(const T *v, T *r, octave_idx_type n)
Definition mx-inlines.cc:988
OP_CUM_FCN
#define OP_CUM_FCN(F, TSRC, TRES, OP)
Definition mx-inlines.cc:870
mx_inline_iszero
void mx_inline_iszero(std::size_t n, bool *r, const X *x)
Definition mx-inlines.cc:72
mx_inline_xmax
void mx_inline_xmax(std::size_t n, T *r, const T *x, const T *y)
Definition mx-inlines.cc:365
mx_inline_any_r
void mx_inline_any_r(const T *v, bool *r, octave_idx_type m, octave_idx_type n)
Definition mx-inlines.cc:796
xis_true
bool xis_true(T x)
Definition mx-inlines.cc:605
mx_inline_count
T mx_inline_count(const bool *v, octave_idx_type n)
Definition mx-inlines.cc:762
mx_inline_and
void mx_inline_and(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:219
OP_MINMAX_FCNN
#define OP_MINMAX_FCNN(F)
Definition mx-inlines.cc:1083
mx_inline_add
void mx_inline_add(std::size_t n, R *r, const X *x, const Y *y)
Definition mx-inlines.cc:108
mx_inline_not
void mx_inline_not(std::size_t n, bool *r, const X *x)
Definition mx-inlines.cc:183
OP_CUM_FCNN
#define OP_CUM_FCNN(F, TSRC, TRES)
Definition mx-inlines.cc:911
do_ms_inplace_op
Array< R > & do_ms_inplace_op(Array< R > &r, const X &x, void(*op)(std::size_t, R *, X))
Definition mx-inlines.cc:568
mx_inline_not_or
void mx_inline_not_or(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:222
do_mm_binary_op
Array< R > do_mm_binary_op(const Array< X > &x, const Array< Y > &y, void(*op)(std::size_t, R *, const X *, const Y *), void(*op1)(std::size_t, R *, X, const Y *), void(*op2)(std::size_t, R *, const X *, Y), const char *opname)
Definition mx-inlines.cc:505
mx_inline_cummax
void mx_inline_cummax(const T *v, T *r, octave_idx_type n)
Definition mx-inlines.cc:1210
OP_CUMMINMAX_FCN2
#define OP_CUMMINMAX_FCN2(F, OP)
Definition mx-inlines.cc:1216
do_sm_binary_op
Array< R > do_sm_binary_op(const X &x, const Array< Y > &y, void(*op)(std::size_t, R *, X, const Y *))
Definition mx-inlines.cc:539
op_dble_prod
void op_dble_prod(double &ac, float el)
Definition mx-inlines.cc:688
cabsq
T cabsq(const std::complex< T > &c)
Definition mx-inlines.cc:597
OP_RED_FCN2
#define OP_RED_FCN2(F, TSRC, TRES, OP, ZERO)
Definition mx-inlines.cc:770
mx_inline_gt
void mx_inline_gt(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:154
op_dble_sum
void op_dble_sum(double &ac, float el)
Definition mx-inlines.cc:708
mx_inline_all_r
void mx_inline_all_r(const T *v, bool *r, octave_idx_type m, octave_idx_type n)
Definition mx-inlines.cc:797
DEFMINMAXSPEC
#define DEFMINMAXSPEC(T, F, OP)
Definition mx-inlines.cc:388
mx_inline_cumcount
void mx_inline_cumcount(const bool *v, T *r, octave_idx_type n)
Definition mx-inlines.cc:885
mx_inline_mul2
void mx_inline_mul2(std::size_t n, R *r, const X *x)
Definition mx-inlines.cc:129
OP_RED_ANYC
#define OP_RED_ANYC(ac, el)
Definition mx-inlines.cc:728
mx_inline_lt
void mx_inline_lt(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:152
mx_inline_all
bool mx_inline_all(const T *v, octave_idx_type n)
Definition mx-inlines.cc:768
mx_inline_equal
bool mx_inline_equal(std::size_t n, const T1 *x, const T2 *y)
Definition mx-inlines.cc:577
mx_inline_real
void mx_inline_real(std::size_t n, T *r, const std::complex< T > *x)
Definition mx-inlines.cc:325
mx_inline_sumsq
T mx_inline_sumsq(const T *v, octave_idx_type n)
Definition mx-inlines.cc:765
do_mx_unary_map
Array< R > do_mx_unary_map(const Array< X > &x)
Definition mx-inlines.cc:482
OP_CUMMINMAX_FCNN
#define OP_CUMMINMAX_FCNN(F)
Definition mx-inlines.cc:1346
mx_inline_any_positive
bool mx_inline_any_positive(std::size_t n, const T *x)
Definition mx-inlines.cc:299
mx_inline_imag
void mx_inline_imag(std::size_t n, T *r, const std::complex< T > *x)
Definition mx-inlines.cc:333
mx_inline_sum
T mx_inline_sum(const T *v, octave_idx_type n)
Definition mx-inlines.cc:760
mx_inline_min
void mx_inline_min(const T *v, T *r, octave_idx_type n)
Definition mx-inlines.cc:987
OP_RED_SUM
#define OP_RED_SUM(ac, el)
Definition mx-inlines.cc:682
mx_inline_prod
T mx_inline_prod(const T *v, octave_idx_type n)
Definition mx-inlines.cc:763
OP_MINMAX_FCN
#define OP_MINMAX_FCN(F, OP)
Definition mx-inlines.cc:941
OP_RED_FCN
#define OP_RED_FCN(F, TSRC, TRES, OP, ZERO)
Definition mx-inlines.cc:746
mx_inline_diff
void mx_inline_diff(const T *v, T *r, octave_idx_type n, octave_idx_type order)
Definition mx-inlines.cc:1406
mx_inline_or2
void mx_inline_or2(std::size_t n, bool *r, const X *x)
Definition mx-inlines.cc:244
mx_inline_all_finite
bool mx_inline_all_finite(std::size_t n, const T *x)
Definition mx-inlines.cc:273
mx_inline_mul
void mx_inline_mul(std::size_t n, R *r, const X *x, const Y *y)
Definition mx-inlines.cc:110
OP_RED_PROD
#define OP_RED_PROD(ac, el)
Definition mx-inlines.cc:683
mx_inline_any_negative
bool mx_inline_any_negative(std::size_t n, const T *x)
Definition mx-inlines.cc:286
logical_value
bool logical_value(T x)
Definition mx-inlines.cc:162
get_extent_triplet
void get_extent_triplet(const dim_vector &dims, int &dim, octave_idx_type &l, octave_idx_type &n, octave_idx_type &u)
Definition mx-inlines.cc:1516
mx_inline_pow
void mx_inline_pow(std::size_t n, R *r, const X *x, const Y *y)
Definition mx-inlines.cc:419
do_mx_diff_op
Array< R > do_mx_diff_op(const Array< R > &src, int dim, octave_idx_type order, void(*mx_diff_op)(const R *, R *, octave_idx_type, octave_idx_type, octave_idx_type, octave_idx_type))
Definition mx-inlines.cc:1666
mx_inline_fill
void mx_inline_fill(std::size_t n, R *r, S s)
Definition mx-inlines.cc:48
mx_inline_not2
void mx_inline_not2(std::size_t n, bool *r)
Definition mx-inlines.cc:190
mx_inline_ne
void mx_inline_ne(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:157
OP_RED_ALLR
#define OP_RED_ALLR(ac, el)
Definition mx-inlines.cc:794
mx_inline_dsum
subst_template_param< std::complex, T, double >::type mx_inline_dsum(const T *v, octave_idx_type n)
Definition mx-inlines.cc:761
OP_RED_SUMSQ
#define OP_RED_SUMSQ(ac, el)
Definition mx-inlines.cc:684
do_mx_minmax_op
Array< R > do_mx_minmax_op(const Array< R > &src, int dim, void(*mx_minmax_op)(const R *, R *, octave_idx_type, octave_idx_type, octave_idx_type))
Definition mx-inlines.cc:1588
mx_inline_or_not
void mx_inline_or_not(std::size_t n, bool *r, const X *x, const Y *y)
Definition mx-inlines.cc:224
DEFMXBINOP
#define DEFMXBINOP(F, OP)
Definition mx-inlines.cc:88
Complex
std::complex< double > Complex
Definition oct-cmplx.h:33
FloatComplex
std::complex< float > FloatComplex
Definition oct-cmplx.h:34
pow
octave_int< T > pow(const octave_int< T > &a, const octave_int< T > &b)
Definition oct-inttypes.cc:740
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition oct-locbuf.h:44