GNU Octave 11.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-2026 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#include <limits>
37
38#include "Array-util.h"
39#include "Array-oct.h"
40#include "bsxfun.h"
41#include "oct-cmplx.h"
42#include "oct-inttypes-fwd.h"
43#include "oct-locbuf.h"
44
45// Provides some commonly repeated, basic loop templates.
46
47template <typename R, typename S>
48inline void
49mx_inline_fill (std::size_t n, R *r, S s)
50{
51 for (std::size_t i = 0; i < n; i++)
52 r[i] = s;
53}
54
55template <typename R, typename X>
56inline void
57mx_inline_uminus (std::size_t n, R *r, const X *x)
58{
59 for (std::size_t i = 0; i < n; i++)
60 r[i] = -x[i];
61}
62
63template <typename R>
64inline void
65mx_inline_uminus2 (std::size_t n, R *r)
66{
67 for (std::size_t i = 0; i < n; i++)
68 r[i] = -r[i];
69}
70
71template <typename X>
72inline void
73mx_inline_iszero (std::size_t n, bool *r, const X *x)
74{
75 const X zero = X ();
76 for (std::size_t i = 0; i < n; i++)
77 r[i] = x[i] == zero;
78}
79
80template <typename X>
81inline void
82mx_inline_notzero (std::size_t n, bool *r, const X *x)
83{
84 const X zero = X ();
85 for (std::size_t i = 0; i < n; i++)
86 r[i] = x[i] != zero;
87}
88
89#define DEFMXBINOP(F, OP) \
90 template <typename R, typename X, typename Y> \
91 inline void F (std::size_t n, R *r, const X *x, const Y *y) \
92 { \
93 for (std::size_t i = 0; i < n; i++) \
94 r[i] = x[i] OP y[i]; \
95 } \
96 template <typename R, typename X, typename Y> \
97 inline void F (std::size_t n, R *r, const X *x, Y y) \
98 { \
99 for (std::size_t i = 0; i < n; i++) \
100 r[i] = x[i] OP y; \
101 } \
102 template <typename R, typename X, typename Y> \
103 inline void F (std::size_t n, R *r, X x, const Y *y) \
104 { \
105 for (std::size_t i = 0; i < n; i++) \
106 r[i] = x OP y[i]; \
107 }
108
113
114#define DEFMXBINOPEQ(F, OP) \
115 template <typename R, typename X> \
116 inline void F (std::size_t n, R *r, const X *x) \
117 { \
118 for (std::size_t i = 0; i < n; i++) \
119 r[i] OP x[i]; \
120 } \
121 template <typename R, typename X> \
122 inline void F (std::size_t n, R *r, X x) \
123 { \
124 for (std::size_t i = 0; i < n; i++) \
125 r[i] OP x; \
126 }
127
132
133#define DEFMXCMPOP(F, OP) \
134 template <typename X, typename Y> \
135 inline void F (std::size_t n, bool *r, const X *x, const Y *y) \
136 { \
137 for (std::size_t i = 0; i < n; i++) \
138 r[i] = x[i] OP y[i]; \
139 } \
140 template <typename X, typename Y> \
141 inline void F (std::size_t n, bool *r, const X *x, Y y) \
142 { \
143 for (std::size_t i = 0; i < n; i++) \
144 r[i] = x[i] OP y; \
145 } \
146 template <typename X, typename Y> \
147 inline void F (std::size_t n, bool *r, X x, const Y *y) \
148 { \
149 for (std::size_t i = 0; i < n; i++) \
150 r[i] = x OP y[i]; \
151 }
152
159
160// Convert to logical value, for logical op purposes.
161template <typename T>
162inline bool
164{
165 return x;
166}
167
168template <typename T>
169inline bool
170logical_value (const std::complex<T>& x)
171{
172 return x.real () != 0 || x.imag () != 0;
173}
174
175template <typename T>
176inline bool
178{
179 return x.value ();
180}
181
182template <typename X>
183void
184mx_inline_not (std::size_t n, bool *r, const X *x)
185{
186 for (std::size_t i = 0; i < n; i++)
187 r[i] = ! logical_value (x[i]);
188}
189
190inline void
191mx_inline_not2 (std::size_t n, bool *r)
192{
193 for (std::size_t i = 0; i < n; i++)
194 r[i] = ! r[i];
195}
196
197#define DEFMXBOOLOP(F, NOT1, OP, NOT2) \
198 template <typename X, typename Y> \
199 inline void F (std::size_t n, bool *r, const X *x, const Y *y) \
200 { \
201 for (std::size_t i = 0; i < n; i++) \
202 r[i] = ((NOT1 logical_value (x[i])) \
203 OP (NOT2 logical_value (y[i]))); \
204 } \
205 template <typename X, typename Y> \
206 inline void F (std::size_t n, bool *r, const X *x, Y y) \
207 { \
208 const bool yy = (NOT2 logical_value (y)); \
209 for (std::size_t i = 0; i < n; i++) \
210 r[i] = (NOT1 logical_value (x[i])) OP yy; \
211 } \
212 template <typename X, typename Y> \
213 inline void F (std::size_t n, bool *r, X x, const Y *y) \
214 { \
215 const bool xx = (NOT1 logical_value (x)); \
216 for (std::size_t i = 0; i < n; i++) \
217 r[i] = xx OP (NOT2 logical_value (y[i])); \
218 }
219
226
227template <typename X>
228inline void
229mx_inline_and2 (std::size_t n, bool *r, const X *x)
230{
231 for (std::size_t i = 0; i < n; i++)
232 r[i] &= logical_value (x[i]);
233}
234
235template <typename X>
236inline void
237mx_inline_and2 (std::size_t n, bool *r, X x)
238{
239 for (std::size_t i = 0; i < n; i++)
240 r[i] &= x;
241}
242
243template <typename X>
244inline void
245mx_inline_or2 (std::size_t n, bool *r, const X *x)
246{
247 for (std::size_t i = 0; i < n; i++)
248 r[i] |= logical_value (x[i]);
249}
250
251template <typename X>
252inline void
253mx_inline_or2 (std::size_t n, bool *r, X x)
254{
255 for (std::size_t i = 0; i < n; i++)
256 r[i] |= x;
257}
258
259template <typename T>
260inline bool
261mx_inline_any_nan (std::size_t n, const T *x)
262{
263 for (std::size_t i = 0; i < n; i++)
264 {
265 if (octave::math::isnan (x[i]))
266 return true;
267 }
268
269 return false;
270}
271
272template <typename T>
273inline bool
274mx_inline_all_finite (std::size_t n, const T *x)
275{
276 for (std::size_t i = 0; i < n; i++)
277 {
278 if (! octave::math::isfinite (x[i]))
279 return false;
280 }
281
282 return true;
283}
284
285template <typename T>
286inline bool
287mx_inline_any_negative (std::size_t n, const T *x)
288{
289 for (std::size_t i = 0; i < n; i++)
290 {
291 if (x[i] < 0)
292 return true;
293 }
294
295 return false;
296}
297
298template <typename T>
299inline bool
300mx_inline_any_positive (std::size_t n, const T *x)
301{
302 for (std::size_t i = 0; i < n; i++)
303 {
304 if (x[i] > 0)
305 return true;
306 }
307
308 return false;
309}
310
311template <typename T>
312inline bool
313mx_inline_all_real (std::size_t n, const std::complex<T> *x)
314{
315 for (std::size_t i = 0; i < n; i++)
316 {
317 if (x[i].imag () != 0)
318 return false;
319 }
320
321 return true;
322}
323
324template <typename T>
325inline void
326mx_inline_real (std::size_t n, T *r, const std::complex<T> *x)
327{
328 for (std::size_t i = 0; i < n; i++)
329 r[i] = x[i].real ();
330}
331
332template <typename T>
333inline void
334mx_inline_imag (std::size_t n, T *r, const std::complex<T> *x)
335{
336 for (std::size_t i = 0; i < n; i++)
337 r[i] = x[i].imag ();
338}
339
340
341template <typename T>
342inline void
343mx_inline_xmin (std::size_t n, T *r, const T *x, const T *y)
344{
345 for (std::size_t i = 0; i < n; i++)
346 r[i] = octave::math::min (x[i], y[i]);
347}
348
349template <typename T>
350inline void
351mx_inline_xmin (std::size_t n, T *r, const T *x, const T *y, const bool nanflag)
352{
353 for (std::size_t i = 0; i < n; i++)
354 r[i] = octave::math::min (x[i], y[i], nanflag);
355}
356
357template <typename T>
358inline void
359mx_inline_xmin (std::size_t n, T *r, const T *x, const T *y,
360 const bool nanflag, const bool realabs)
361{
362 for (std::size_t i = 0; i < n; i++)
363 r[i] = octave::math::min (x[i], y[i], nanflag, realabs);
364}
365
366template <typename T>
367inline void
368mx_inline_xmin (std::size_t n, T *r, const T *x, T y)
369{
370 for (std::size_t i = 0; i < n; i++)
371 r[i] = octave::math::min (x[i], y);
372}
373
374template <typename T>
375inline void
376mx_inline_xmin (std::size_t n, T *r, const T *x, T y, const bool nanflag)
377{
378 for (std::size_t i = 0; i < n; i++)
379 r[i] = octave::math::min (x[i], y, nanflag);
380}
381
382template <typename T>
383inline void
384mx_inline_xmin (std::size_t n, T *r, const T *x, T y,
385 const bool nanflag, const bool realabs)
386{
387 for (std::size_t i = 0; i < n; i++)
388 r[i] = octave::math::min (x[i], y, nanflag, realabs);
389}
390
391template <typename T>
392inline void
393mx_inline_xmin (std::size_t n, T *r, T x, const T *y)
394{
395 for (std::size_t i = 0; i < n; i++)
396 r[i] = octave::math::min (x, y[i]);
397}
398
399template <typename T>
400inline void
401mx_inline_xmin (std::size_t n, T *r, T x, const T *y, const bool nanflag)
402{
403 for (std::size_t i = 0; i < n; i++)
404 r[i] = octave::math::min (x, y[i], nanflag);
405}
406
407template <typename T>
408inline void
409mx_inline_xmin (std::size_t n, T *r, T x, const T *y,
410 const bool nanflag, const bool realabs)
411{
412 for (std::size_t i = 0; i < n; i++)
413 r[i] = octave::math::min (x, y[i], nanflag, realabs);
414}
415
416template <typename T>
417inline void
418mx_inline_xmax (std::size_t n, T *r, const T *x, const T *y)
419{
420 for (std::size_t i = 0; i < n; i++)
421 r[i] = octave::math::max (x[i], y[i]);
422}
423
424template <typename T>
425inline void
426mx_inline_xmax (std::size_t n, T *r, const T *x, const T *y, const bool nanflag)
427{
428 for (std::size_t i = 0; i < n; i++)
429 r[i] = octave::math::max (x[i], y[i], nanflag);
430}
431
432template <typename T>
433inline void
434mx_inline_xmax (std::size_t n, T *r, const T *x, const T *y,
435 const bool nanflag, const bool realabs)
436{
437 for (std::size_t i = 0; i < n; i++)
438 r[i] = octave::math::max (x[i], y[i], nanflag, realabs);
439}
440
441template <typename T>
442inline void
443mx_inline_xmax (std::size_t n, T *r, const T *x, T y)
444{
445 for (std::size_t i = 0; i < n; i++)
446 r[i] = octave::math::max (x[i], y);
447}
448
449template <typename T>
450inline void
451mx_inline_xmax (std::size_t n, T *r, const T *x, T y, const bool nanflag)
452{
453 for (std::size_t i = 0; i < n; i++)
454 r[i] = octave::math::max (x[i], y, nanflag);
455}
456
457template <typename T>
458inline void
459mx_inline_xmax (std::size_t n, T *r, const T *x, T y,
460 const bool nanflag, const bool realabs)
461{
462 for (std::size_t i = 0; i < n; i++)
463 r[i] = octave::math::max (x[i], y, nanflag, realabs);
464}
465
466template <typename T>
467inline void
468mx_inline_xmax (std::size_t n, T *r, T x, const T *y)
469{
470 for (std::size_t i = 0; i < n; i++)
471 r[i] = octave::math::max (x, y[i]);
472}
473
474template <typename T>
475inline void
476mx_inline_xmax (std::size_t n, T *r, T x, const T *y, const bool nanflag)
477{
478 for (std::size_t i = 0; i < n; i++)
479 r[i] = octave::math::max (x, y[i], nanflag);
480}
481
482template <typename T>
483inline void
484mx_inline_xmax (std::size_t n, T *r, T x, const T *y,
485 const bool nanflag, const bool realabs)
486{
487 for (std::size_t i = 0; i < n; i++)
488 r[i] = octave::math::max (x, y[i], nanflag, realabs);
489}
490
491// FIXME: Is this comment correct anymore? It seems like std::pow is chosen.
492// Let the compiler decide which pow to use, whichever best matches the
493// arguments provided.
494
495template <typename R, typename X, typename Y>
496inline void
497mx_inline_pow (std::size_t n, R *r, const X *x, const Y *y)
498{
499 using std::pow;
500
501 for (std::size_t i = 0; i < n; i++)
502 r[i] = pow (x[i], y[i]);
503}
504
505template <typename R, typename X, typename Y>
506inline void
507mx_inline_pow (std::size_t n, R *r, const X *x, Y y)
508{
509 using std::pow;
510
511 for (std::size_t i = 0; i < n; i++)
512 r[i] = pow (x[i], y);
513}
514
515template <typename R, typename X, typename Y>
516inline void
517mx_inline_pow (std::size_t n, R *r, X x, const Y *y)
518{
519 using std::pow;
520
521 for (std::size_t i = 0; i < n; i++)
522 r[i] = pow (x, y[i]);
523}
524
525// Arbitrary function appliers.
526// The function is a template parameter to enable inlining.
527template <typename R, typename X, R fcn (X x)>
528inline void
529mx_inline_map (std::size_t n, R *r, const X *x)
530{
531 for (std::size_t i = 0; i < n; i++)
532 r[i] = fcn (x[i]);
533}
534
535template <typename R, typename X, R fcn (const X& x)>
536inline void
537mx_inline_map (std::size_t n, R *r, const X *x)
538{
539 for (std::size_t i = 0; i < n; i++)
540 r[i] = fcn (x[i]);
541}
542
543// Appliers. Since these call the operation just once, we pass it as
544// a pointer, to allow the compiler reduce number of instances.
545
546template <typename R, typename X>
547inline Array<R>
549 void (*op) (std::size_t, R *, const X *))
550{
551 Array<R> r (x.dims ());
552 op (r.numel (), r.rwdata (), x.data ());
553 return r;
554}
555
556// Shortcuts for applying mx_inline_map.
557
558template <typename R, typename X, R fcn (X)>
559inline Array<R>
561{
562 return do_mx_unary_op<R, X> (x, mx_inline_map<R, X, fcn>);
563}
564
565template <typename R, typename X, R fcn (const X&)>
566inline Array<R>
568{
569 return do_mx_unary_op<R, X> (x, mx_inline_map<R, X, fcn>);
570}
571
572template <typename R>
573inline Array<R>&
575 void (*op) (std::size_t, R *))
576{
577 op (r.numel (), r.rwdata ());
578 return r;
579}
580
581template <typename R, typename X, typename Y>
582inline Array<R>
584 void (*op) (std::size_t, R *, const X *, const Y *),
585 void (*op1) (std::size_t, R *, X, const Y *),
586 void (*op2) (std::size_t, R *, const X *, Y),
587 const char *opname)
588{
589 const dim_vector& dx = x.dims ();
590 const dim_vector& dy = y.dims ();
591 if (dx == dy)
592 {
593 Array<R> r (dx);
594 op (r.numel (), r.rwdata (), x.data (), y.data ());
595 return r;
596 }
597 else if (is_valid_bsxfun (dx, dy))
598 {
599 return do_bsxfun_op (x, y, op, op1, op2);
600 }
601 else
602 octave::err_nonconformant (opname, dx, dy);
603}
604
605template <typename R, typename X, typename Y>
606inline Array<R>
607do_ms_binary_op (const Array<X>& x, const Y& y,
608 void (*op) (std::size_t, R *, const X *, Y))
609{
610 Array<R> r (x.dims ());
611 op (r.numel (), r.rwdata (), x.data (), y);
612 return r;
613}
614
615template <typename R, typename X, typename Y>
616inline Array<R>
617do_sm_binary_op (const X& x, const Array<Y>& y,
618 void (*op) (std::size_t, R *, X, const Y *))
619{
620 Array<R> r (y.dims ());
621 op (r.numel (), r.rwdata (), x, y.data ());
622 return r;
623}
624
625// These are required for min/max binary operations with extra
626// nanflag and realabs arguments
627
628template <typename R, typename X, typename Y>
629inline Array<R>
630do_mm_binary_op (const Array<X>& x, const Array<Y>& y, bool nanflag,
631 void (*op) (std::size_t, R *, const X *, const Y *, bool),
632 void (*op1) (std::size_t, R *, X, const Y *, bool),
633 void (*op2) (std::size_t, R *, const X *, Y, bool),
634 const char *opname)
635{
636 const dim_vector& dx = x.dims ();
637 const dim_vector& dy = y.dims ();
638 if (dx == dy)
639 {
640 Array<R> r (dx);
641 op (r.numel (), r.rwdata (), x.data (), y.data (), nanflag);
642 return r;
643 }
644 else if (is_valid_bsxfun (dx, dy))
645 {
646 return do_bsxfun1_op (x, y, nanflag, op, op1, op2);
647 }
648 else
649 octave::err_nonconformant (opname, dx, dy);
650}
651
652template <typename R, typename X, typename Y>
653inline Array<R>
654do_ms_binary_op (const Array<X>& x, const Y& y, bool nanflag,
655 void (*op) (std::size_t, R *, const X *, Y, bool))
656{
657 Array<R> r (x.dims ());
658 op (r.numel (), r.rwdata (), x.data (), y, nanflag);
659 return r;
660}
661
662template <typename R, typename X, typename Y>
663inline Array<R>
664do_sm_binary_op (const X& x, const Array<Y>& y, bool nanflag,
665 void (*op) (std::size_t, R *, X, const Y *, bool))
666{
667 Array<R> r (y.dims ());
668 op (r.numel (), r.rwdata (), x, y.data (), nanflag);
669 return r;
670}
671
672template <typename R, typename X, typename Y>
673inline Array<R>
675 bool nanflag, bool realabs,
676 void (*op) (std::size_t, R *, const X *, const Y *,
677 bool, bool),
678 void (*op1) (std::size_t, R *, X, const Y *, bool, bool),
679 void (*op2) (std::size_t, R *, const X *, Y, bool, bool),
680 const char *opname)
681{
682 const dim_vector& dx = x.dims ();
683 const dim_vector& dy = y.dims ();
684 if (dx == dy)
685 {
686 Array<R> r (dx);
687 op (r.numel (), r.rwdata (), x.data (), y.data (), nanflag, realabs);
688 return r;
689 }
690 else if (is_valid_bsxfun (dx, dy))
691 {
692 return do_bsxfun2_op (x, y, nanflag, realabs, op, op1, op2);
693 }
694 else
695 octave::err_nonconformant (opname, dx, dy);
696}
697
698template <typename R, typename X, typename Y>
699inline Array<R>
700do_ms_binary_op (const Array<X>& x, const Y& y, bool nanflag, bool realabs,
701 void (*op) (std::size_t, R *, const X *, Y, bool, bool))
702{
703 Array<R> r (x.dims ());
704 op (r.numel (), r.rwdata (), x.data (), y, nanflag, realabs);
705 return r;
706}
707
708template <typename R, typename X, typename Y>
709inline Array<R>
710do_sm_binary_op (const X& x, const Array<Y>& y, bool nanflag, bool realabs,
711 void (*op) (std::size_t, R *, X, const Y *, bool, bool))
712{
713 Array<R> r (y.dims ());
714 op (r.numel (), r.rwdata (), x, y.data (), nanflag, realabs);
715 return r;
716}
717
718template <typename R, typename X>
719inline Array<R>&
721 void (*op) (std::size_t, R *, const X *),
722 void (*op1) (std::size_t, R *, X),
723 const char *opname)
724{
725 const dim_vector& dr = r.dims ();
726 const dim_vector& dx = x.dims ();
727 if (dr == dx)
728 op (r.numel (), r.rwdata (), x.data ());
729 else if (is_valid_inplace_bsxfun (dr, dx))
730 do_inplace_bsxfun_op (r, x, op, op1);
731 else
732 octave::err_nonconformant (opname, dr, dx);
733
734 return r;
735}
736
737template <typename R, typename X>
738inline Array<R>&
740 void (*op) (std::size_t, R *, X))
741{
742 op (r.numel (), r.rwdata (), x);
743 return r;
744}
745
746template <typename T1, typename T2>
747inline bool
748mx_inline_equal (std::size_t n, const T1 *x, const T2 *y)
749{
750 for (std::size_t i = 0; i < n; i++)
751 if (x[i] != y[i])
752 return false;
753 return true;
754}
755
756template <typename T>
757inline bool
759 bool (*op) (std::size_t, const T *))
760{
761 return op (a.numel (), a.data ());
762}
763
764// NOTE: we don't use std::norm because it typically does some heavyweight
765// magic to avoid underflows, which we don't need here.
766template <typename T>
767inline T
768cabsq (const std::complex<T>& c)
769{
770 return c.real () * c.real () + c.imag () * c.imag ();
771}
772
773// default. works for integers and bool.
774template <typename T>
775inline bool
777{
778 return x;
779}
780
781template <typename T>
782inline bool
784{
785 return ! x;
786}
787
788// for octave_ints
789template <typename T>
790inline bool
792{
793 return x.value ();
794}
795
796template <typename T>
797inline bool
799{
800 return ! x.value ();
801}
802
803// for reals, we want to ignore NaNs.
804inline bool
805xis_true (double x)
806{
807 return ! octave::math::isnan (x) && x != 0.0;
808}
809
810inline bool
811xis_false (double x)
812{
813 return x == 0.0;
814}
815
816inline bool
817xis_true (float x)
818{
819 return ! octave::math::isnan (x) && x != 0.0f;
820}
821
822inline bool
824{
825 return x == 0.0f;
826}
827
828// Ditto for complex.
829inline bool
831{
832 return ! octave::math::isnan (x) && x != 0.0;
833}
834
835inline bool
837{
838 return x == 0.0;
839}
840
841inline bool
843{
844 return ! octave::math::isnan (x) && x != 0.0f;
845}
846
847inline bool
849{
850 return x == 0.0f;
851}
852
853#define OP_RED_SUM(ac, el) ac += el
854#define OP_RED_PROD(ac, el) ac *= el
855#define OP_RED_SUMSQ(ac, el) ac += ((el)*(el))
856#define OP_RED_SUMSQC(ac, el) ac += cabsq (el)
857
858inline void
859op_dble_prod (double& ac, float el)
860{
861 ac *= el;
862}
863
864// FIXME: guaranteed?
865inline void
867{
868 ac *= el;
869}
870
871template <typename T>
872inline void
873op_dble_prod (double& ac, const octave_int<T>& el)
874{
875 ac *= el.double_value ();
876}
877
878inline void
879op_dble_sum (double& ac, float el)
880{
881 ac += el;
882}
883
884// FIXME: guaranteed?
885inline void
887{
888 ac += el;
889}
890
891template <typename T>
892inline void
893op_dble_sum (double& ac, const octave_int<T>& el)
894{
895 ac += el.double_value ();
896}
897
898// The following two implement a simple short-circuiting.
899#define OP_RED_ANYC(ac, el) \
900 if (xis_true (el)) \
901 { \
902 ac = true; \
903 break; \
904 } \
905 else \
906 continue
907
908#define OP_RED_ALLC(ac, el) \
909 if (xis_false (el)) \
910 { \
911 ac = false; \
912 break; \
913 } \
914 else \
915 continue
916
917#define OP_RED_FCN(F, TSRC, TRES, OP, ZERO) \
918 template <typename T> \
919 inline TRES \
920 F (const TSRC *v, octave_idx_type n) \
921 { \
922 TRES ac = ZERO; \
923 for (octave_idx_type i = 0; i < n; i++) \
924 OP(ac, v[i]); \
925 return ac; \
926 }
927
928#define PROMOTE_DOUBLE(T) \
929 typename subst_template_param<std::complex, T, double>::type
930
937OP_RED_FCN (mx_inline_sumsq, std::complex<T>, T, OP_RED_SUMSQC, 0)
942
943#define OP_RED_NAN_FCN(F, TSRC, TRES, OP, ZERO) \
944 template <typename T> \
945 inline TRES \
946 F (const TSRC *v, octave_idx_type n, bool nanflag) \
947 { \
948 TRES ac = ZERO; \
949 if (nanflag) \
950 { \
951 for (octave_idx_type i = 0; i < n; i++) \
952 if (! octave::math::isnan (v[i])) \
953 OP(ac, v[i]); \
954 } \
955 else \
956 { \
957 for (octave_idx_type i = 0; i < n; i++) \
958 OP(ac, v[i]); \
959 } \
960 return ac; \
961 }
962
971
972#define OP_RED_FCN2(F, TSRC, TRES, OP, ZERO) \
973 template <typename T> \
974 inline void \
975 F (const TSRC *v, TRES *r, octave_idx_type m, octave_idx_type n) \
976 { \
977 for (octave_idx_type i = 0; i < m; i++) \
978 r[i] = ZERO; \
979 for (octave_idx_type j = 0; j < n; j++) \
980 { \
981 for (octave_idx_type i = 0; i < m; i++) \
982 OP(r[i], v[i]); \
983 v += m; \
984 } \
985 }
986
996
997#define OP_RED_NAN_FCN2(F, TSRC, TRES, OP, ZERO) \
998 template <typename T> \
999 inline void \
1000 F (const TSRC *v, TRES *r, octave_idx_type m, octave_idx_type n, \
1001 bool nanflag) \
1002 { \
1003 for (octave_idx_type i = 0; i < m; i++) \
1004 r[i] = ZERO; \
1005 \
1006 if (nanflag) \
1007 { \
1008 /* omitnan: skip NaN, result stays at identity if all NaN */ \
1009 for (octave_idx_type j = 0; j < n; j++) \
1010 { \
1011 for (octave_idx_type i = 0; i < m; i++) \
1012 if (! octave::math::isnan (v[i])) \
1013 OP(r[i], v[i]); \
1014 v += m; \
1015 } \
1016 } \
1017 else \
1018 { \
1019 /* includenan: NaN propagates naturally */ \
1020 for (octave_idx_type j = 0; j < n; j++) \
1021 { \
1022 for (octave_idx_type i = 0; i < m; i++) \
1023 OP(r[i], v[i]); \
1024 v += m; \
1025 } \
1026 } \
1027 }
1028
1037
1038#define OP_RED_ANYR(ac, el) ac |= xis_true (el)
1039#define OP_RED_ALLR(ac, el) ac &= xis_true (el)
1040
1043
1044// Using the general code for any/all would sacrifice short-circuiting.
1045// OTOH, going by rows would sacrifice cache-coherence. The following
1046// algorithm will achieve both, at the cost of a temporary octave_idx_type
1047// array.
1048
1049#define OP_ROW_SHORT_CIRCUIT(F, PRED, ZERO) \
1050 template <typename T> \
1051 inline void \
1052 F (const T *v, bool *r, octave_idx_type m, octave_idx_type n) \
1053 { \
1054 if (n <= 8) \
1055 return F ## _r (v, r, m, n); \
1056 \
1057 /* FIXME: it may be sub-optimal to allocate the buffer here. */ \
1058 OCTAVE_LOCAL_BUFFER (octave_idx_type, iact, m); \
1059 for (octave_idx_type i = 0; i < m; i++) iact[i] = i; \
1060 octave_idx_type nact = m; \
1061 for (octave_idx_type j = 0; j < n; j++) \
1062 { \
1063 octave_idx_type k = 0; \
1064 for (octave_idx_type i = 0; i < nact; i++) \
1065 { \
1066 octave_idx_type ia = iact[i]; \
1067 if (! PRED (v[ia])) \
1068 iact[k++] = ia; \
1069 } \
1070 nact = k; \
1071 v += m; \
1072 } \
1073 for (octave_idx_type i = 0; i < m; i++) r[i] = ! ZERO; \
1074 for (octave_idx_type i = 0; i < nact; i++) r[iact[i]] = ZERO; \
1075 }
1076
1079
1080#define OP_RED_FCNN(F, TSRC, TRES) \
1081 template <typename T> \
1082 inline void \
1083 F (const TSRC *v, TRES *r, octave_idx_type l, \
1084 octave_idx_type n, octave_idx_type u) \
1085 { \
1086 if (l == 1) \
1087 { \
1088 for (octave_idx_type i = 0; i < u; i++) \
1089 { \
1090 r[i] = F<T> (v, n); \
1091 v += n; \
1092 } \
1093 } \
1094 else \
1095 { \
1096 for (octave_idx_type i = 0; i < u; i++) \
1097 { \
1098 F (v, r, l, n); \
1099 v += l*n; \
1100 r += l; \
1101 } \
1102 } \
1103 }
1104
1111OP_RED_FCNN (mx_inline_sumsq, std::complex<T>, T)
1116
1117#define OP_RED_NAN_FCNN(F, TSRC, TRES) \
1118 template <typename T> \
1119 inline void \
1120 F (const TSRC *v, TRES *r, octave_idx_type l, \
1121 octave_idx_type n, octave_idx_type u, bool nanflag) \
1122 { \
1123 if (l == 1) \
1124 { \
1125 for (octave_idx_type i = 0; i < u; i++) \
1126 { \
1127 r[i] = F<T> (v, n, nanflag); \
1128 v += n; \
1129 } \
1130 } \
1131 else \
1132 { \
1133 for (octave_idx_type i = 0; i < u; i++) \
1134 { \
1135 F (v, r, l, n, nanflag); \
1136 v += l*n; \
1137 r += l; \
1138 } \
1139 } \
1140 }
1141
1150
1151#define OP_CUM_FCN(F, TSRC, TRES, OP) \
1152 template <typename T> \
1153 inline void \
1154 F (const TSRC *v, TRES *r, octave_idx_type n) \
1155 { \
1156 if (n) \
1157 { \
1158 TRES t = r[0] = v[0]; \
1159 for (octave_idx_type i = 1; i < n; i++) \
1160 r[i] = t = t OP v[i]; \
1161 } \
1162 }
1163
1167
1168#define OP_CUM_NAN_FCN(F, TSRC, TRES, OP, ZERO) \
1169 template <typename T> \
1170 inline void \
1171 F (const TSRC *v, TRES *r, octave_idx_type n, bool nanflag) \
1172 { \
1173 if (n) \
1174 { \
1175 if (nanflag) \
1176 { \
1177 TRES z = ZERO; \
1178 TRES t; \
1179 if (octave::math::isnan (v[0])) \
1180 t = r[0] = z; \
1181 else \
1182 t = r[0] = v[0]; \
1183 for (octave_idx_type i = 1; i < n; i++) \
1184 { \
1185 if (octave::math::isnan (v[i])) \
1186 r[i] = t = t OP z; \
1187 else \
1188 r[i] = t = t OP v[i]; \
1189 } \
1190 } \
1191 else \
1192 { \
1193 TRES t = r[0] = v[0]; \
1194 for (octave_idx_type i = 1; i < n; i++) \
1195 r[i] = t = t OP v[i]; \
1196 } \
1197 } \
1198 }
1199
1202
1203#define OP_CUM_FCN2(F, TSRC, TRES, OP) \
1204 template <typename T> \
1205 inline void \
1206 F (const TSRC *v, TRES *r, octave_idx_type m, octave_idx_type n) \
1207 { \
1208 if (n) \
1209 { \
1210 for (octave_idx_type i = 0; i < m; i++) \
1211 r[i] = v[i]; \
1212 const T *r0 = r; \
1213 for (octave_idx_type j = 1; j < n; j++) \
1214 { \
1215 r += m; v += m; \
1216 for (octave_idx_type i = 0; i < m; i++) \
1217 r[i] = r0[i] OP v[i]; \
1218 r0 += m; \
1219 } \
1220 } \
1221 }
1222
1226
1227#define OP_CUM_NAN_FCN2(F, TSRC, TRES, OP, ZERO) \
1228 template <typename T> \
1229 inline void \
1230 F (const TSRC *v, TRES *r, octave_idx_type m, octave_idx_type n, \
1231 bool nanflag) \
1232 { \
1233 if (n) \
1234 { \
1235 if (nanflag) \
1236 { \
1237 TRES z = ZERO; \
1238 for (octave_idx_type i = 0; i < m; i++) \
1239 { \
1240 if (octave::math::isnan (v[i])) \
1241 r[i] = z; \
1242 else \
1243 r[i] = v[i]; \
1244 } \
1245 const T *r0 = r; \
1246 for (octave_idx_type j = 1; j < n; j++) \
1247 { \
1248 r += m; v += m; \
1249 for (octave_idx_type i = 0; i < m; i++) \
1250 { \
1251 if (octave::math::isnan (v[i])) \
1252 r[i] = r0[i] OP z; \
1253 else \
1254 r[i] = r0[i] OP v[i]; \
1255 } \
1256 r0 += m; \
1257 } \
1258 } \
1259 else \
1260 { \
1261 for (octave_idx_type i = 0; i < m; i++) \
1262 r[i] = v[i]; \
1263 const T *r0 = r; \
1264 for (octave_idx_type j = 1; j < n; j++) \
1265 { \
1266 r += m; v += m; \
1267 for (octave_idx_type i = 0; i < m; i++) \
1268 r[i] = r0[i] OP v[i]; \
1269 r0 += m; \
1270 } \
1271 } \
1272 } \
1273 }
1274
1277
1278#define OP_CUM_FCNN(F, TSRC, TRES) \
1279 template <typename T> \
1280 inline void \
1281 F (const TSRC *v, TRES *r, octave_idx_type l, \
1282 octave_idx_type n, octave_idx_type u) \
1283 { \
1284 if (l == 1) \
1285 { \
1286 for (octave_idx_type i = 0; i < u; i++) \
1287 { \
1288 F (v, r, n); \
1289 v += n; \
1290 r += n; \
1291 } \
1292 } \
1293 else \
1294 { \
1295 for (octave_idx_type i = 0; i < u; i++) \
1296 { \
1297 F (v, r, l, n); \
1298 v += l*n; \
1299 r += l*n; \
1300 } \
1301 } \
1302 }
1303
1307
1308#define OP_CUM_NAN_FCNN(F, TSRC, TRES) \
1309 template <typename T> \
1310 inline void \
1311 F (const TSRC *v, TRES *r, octave_idx_type l, \
1312 octave_idx_type n, octave_idx_type u, bool nanflag) \
1313 { \
1314 if (l == 1) \
1315 { \
1316 for (octave_idx_type i = 0; i < u; i++) \
1317 { \
1318 F (v, r, n, nanflag); \
1319 v += n; \
1320 r += n; \
1321 } \
1322 } \
1323 else \
1324 { \
1325 for (octave_idx_type i = 0; i < u; i++) \
1326 { \
1327 F (v, r, l, n, nanflag); \
1328 v += l*n; \
1329 r += l*n; \
1330 } \
1331 } \
1332 }
1333
1336
1337template <typename T>
1338inline void
1339mx_inline_flip (const T *v, T *r, octave_idx_type n)
1340{
1341 if (n)
1342 {
1343 for (octave_idx_type i = 0; i < n; i++)
1344 r[i] = v[n-1-i];
1345 }
1346}
1347
1348template <typename T>
1349inline void
1351{
1352 if (n)
1353 {
1354 for (octave_idx_type j = 0; j < n; j++)
1355 {
1356 for (octave_idx_type i = 0; i < m; i++)
1357 r[(n-1-j)*m+i] = v[j*m+i];
1358 }
1359 }
1360}
1361
1362template <typename T>
1363inline void
1364mx_inline_flip (const T *v, T *r, octave_idx_type l,
1366{
1367 if (l == 1)
1368 {
1369 for (octave_idx_type i = 0; i < u; i++)
1370 {
1371 mx_inline_flip (v, r, n);
1372 v += n;
1373 r += n;
1374 }
1375 }
1376 else
1377 {
1378 for (octave_idx_type i = 0; i < u; i++)
1379 {
1380 mx_inline_flip (v, r, l, n);
1381 v += l*n;
1382 r += l*n;
1383 }
1384 }
1385}
1386
1387// Templated functions for Fmin and Fmax
1388
1389#define OP_MINMAX_FCN(F, OP) \
1390 template <typename T> \
1391 void F (const T *v, T *r, octave_idx_type n, const bool nanflag, \
1392 const bool realabs) \
1393 { \
1394 if (! n) \
1395 return; \
1396 if (! nanflag) \
1397 { \
1398 for (octave_idx_type i = 0; i < n; i++) \
1399 if (octave::math::isnan (v[i])) \
1400 { \
1401 *r = NAN; \
1402 return; \
1403 } \
1404 } \
1405 T tmp = v[0]; \
1406 octave_idx_type i = 1; \
1407 if (octave::math::isnan (tmp)) \
1408 { \
1409 for (; i < n && octave::math::isnan (v[i]); i++) ; \
1410 if (i < n) \
1411 tmp = v[i]; \
1412 } \
1413 if (realabs) \
1414 { \
1415 for (; i < n; i++) \
1416 if (v[i] OP tmp) \
1417 tmp = v[i]; \
1418 } \
1419 else \
1420 { \
1421 for (; i < n; i++) \
1422 { \
1423 if (std::abs (v[i]) OP std::abs (tmp)) \
1424 tmp = v[i]; \
1425 else if (std::abs (v[i]) == std::abs (tmp) && v[i] OP tmp) \
1426 tmp = v[i]; \
1427 } \
1428 } \
1429 *r = tmp; \
1430 } \
1431 template <typename T> \
1432 void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n, \
1433 const bool nanflag, const bool realabs) \
1434 { \
1435 if (! n) \
1436 return; \
1437 if (! nanflag) \
1438 { \
1439 for (octave_idx_type i = 0; i < n; i++) \
1440 if (octave::math::isnan (v[i])) \
1441 { \
1442 *r = NAN; \
1443 *ri = i; \
1444 return; \
1445 } \
1446 } \
1447 T tmp = v[0]; \
1448 octave_idx_type tmpi = 0; \
1449 octave_idx_type i = 1; \
1450 if (octave::math::isnan (tmp)) \
1451 { \
1452 for (; i < n && octave::math::isnan (v[i]); i++) ; \
1453 if (i < n) \
1454 { \
1455 tmp = v[i]; \
1456 tmpi = i; \
1457 } \
1458 } \
1459 if (realabs) \
1460 { \
1461 for (; i < n; i++) \
1462 if (v[i] OP tmp) \
1463 { \
1464 tmp = v[i]; \
1465 tmpi = i; \
1466 } \
1467 } \
1468 else \
1469 { \
1470 for (; i < n; i++) \
1471 { \
1472 if (std::abs (v[i]) OP std::abs (tmp)) \
1473 { \
1474 tmp = v[i]; \
1475 tmpi = i; \
1476 } \
1477 else if (std::abs (v[i]) == std::abs (tmp) && v[i] OP tmp) \
1478 { \
1479 tmp = v[i]; \
1480 tmpi = i; \
1481 } \
1482 } \
1483 } \
1484 *r = tmp; \
1485 *ri = tmpi; \
1486 }
1487
1490
1491// Row reductions will be slightly complicated. We will proceed with checks
1492// for NaNs until we detect that no row will yield a NaN, in which case we
1493// proceed to a faster code.
1494
1495#define OP_MINMAX_FCN2(F, OP) \
1496 template <typename T> \
1497 inline void \
1498 F (const T *v, T *r, octave_idx_type m, octave_idx_type n, \
1499 const bool nanflag, const bool realabs) \
1500 { \
1501 if (! n) \
1502 return; \
1503 octave_idx_type j = 0; \
1504 for (octave_idx_type i = 0; i < m; i++) \
1505 { \
1506 r[i] = v[i]; \
1507 } \
1508 j++; \
1509 v += m; \
1510 if (realabs) \
1511 { \
1512 while (j < n) \
1513 { \
1514 for (octave_idx_type i = 0; i < m; i++) \
1515 { \
1516 if ((octave::math::isnan (v[i]) || \
1517 octave::math::isnan (r[i])) && ! nanflag) \
1518 r[i] = NAN; \
1519 else if (octave::math::isnan (v[i])); \
1520 else if (octave::math::isnan (r[i]) || v[i] OP r[i]) \
1521 r[i] = v[i]; \
1522 } \
1523 j++; \
1524 v += m; \
1525 } \
1526 } \
1527 else \
1528 { \
1529 while (j < n) \
1530 { \
1531 for (octave_idx_type i = 0; i < m; i++) \
1532 { \
1533 if ((octave::math::isnan (v[i]) || \
1534 octave::math::isnan (r[i])) && ! nanflag) \
1535 r[i] = NAN; \
1536 else if (octave::math::isnan (v[i])); \
1537 else if (octave::math::isnan (r[i]) || \
1538 std::abs (v[i]) OP std::abs (r[i])) \
1539 r[i] = v[i]; \
1540 else if (std::abs (v[i]) == std::abs (r[i]) && \
1541 v[i] OP r[i]) \
1542 r[i] = v[i]; \
1543 } \
1544 j++; \
1545 v += m; \
1546 } \
1547 } \
1548 } \
1549 template <typename T> \
1550 inline void \
1551 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type m, \
1552 octave_idx_type n, const bool nanflag, const bool realabs) \
1553 { \
1554 if (! n) \
1555 return; \
1556 octave_idx_type j = 0; \
1557 for (octave_idx_type i = 0; i < m; i++) \
1558 { \
1559 r[i] = v[i]; \
1560 ri[i] = j; \
1561 } \
1562 j++; \
1563 v += m; \
1564 if (realabs) \
1565 { \
1566 while (j < n) \
1567 { \
1568 for (octave_idx_type i = 0; i < m; i++) \
1569 { \
1570 if (octave::math::isnan (r[i]) && ! nanflag); \
1571 else if (octave::math::isnan (v[i]) && ! nanflag) \
1572 { \
1573 r[i] = NAN; \
1574 ri[i] = j; \
1575 } \
1576 else if (octave::math::isnan (v[i])); \
1577 else if (octave::math::isnan (r[i]) || v[i] OP r[i]) \
1578 { \
1579 r[i] = v[i]; \
1580 ri[i] = j; \
1581 } \
1582 } \
1583 j++; \
1584 v += m; \
1585 } \
1586 } \
1587 else \
1588 { \
1589 while (j < n) \
1590 { \
1591 for (octave_idx_type i = 0; i < m; i++) \
1592 { \
1593 if (octave::math::isnan (r[i]) && ! nanflag); \
1594 else if (octave::math::isnan (v[i]) && ! nanflag) \
1595 { \
1596 r[i] = NAN; \
1597 ri[i] = j; \
1598 } \
1599 else if (octave::math::isnan (v[i])); \
1600 else if (octave::math::isnan (r[i]) || \
1601 std::abs (v[i]) OP std::abs (r[i])) \
1602 { \
1603 r[i] = v[i]; \
1604 ri[i] = j; \
1605 } \
1606 else if (std::abs (v[i]) == std::abs (r[i]) && \
1607 v[i] OP r[i]) \
1608 { \
1609 r[i] = v[i]; \
1610 ri[i] = j; \
1611 } \
1612 } \
1613 j++; \
1614 v += m; \
1615 } \
1616 } \
1617 }
1618
1621
1622#define OP_MINMAX_FCNN(F) \
1623 template <typename T> \
1624 inline void \
1625 F (const T *v, T *r, octave_idx_type l, octave_idx_type n, \
1626 octave_idx_type u, const bool nanflag, const bool realabs) \
1627 { \
1628 if (! n) \
1629 return; \
1630 if (l == 1) \
1631 { \
1632 for (octave_idx_type i = 0; i < u; i++) \
1633 { \
1634 F (v, r, n, nanflag, realabs); \
1635 v += n; \
1636 r++; \
1637 } \
1638 } \
1639 else \
1640 { \
1641 for (octave_idx_type i = 0; i < u; i++) \
1642 { \
1643 F (v, r, l, n, nanflag, realabs); \
1644 v += l*n; \
1645 r += l; \
1646 } \
1647 } \
1648 } \
1649 template <typename T> \
1650 inline void \
1651 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type l, \
1652 octave_idx_type n, octave_idx_type u, const bool nanflag, \
1653 const bool realabs) \
1654 { \
1655 if (! n) return; \
1656 if (l == 1) \
1657 { \
1658 for (octave_idx_type i = 0; i < u; i++) \
1659 { \
1660 F (v, r, ri, n, nanflag, realabs); \
1661 v += n; \
1662 r++; \
1663 ri++; \
1664 } \
1665 } \
1666 else \
1667 { \
1668 for (octave_idx_type i = 0; i < u; i++) \
1669 { \
1670 F (v, r, ri, l, n, nanflag, realabs); \
1671 v += l*n; \
1672 r += l; \
1673 ri += l; \
1674 } \
1675 } \
1676 }
1677
1680
1681// Special implementation for complex arrays, since we need comparisons
1682// in both real and imaginary parts for MATLAB compatibility.
1683
1684#define OP_CMINMAX_FCN(F, OP) \
1685 template <typename T> \
1686 void F (const std::complex<T> *v, std::complex<T> *r, \
1687 octave_idx_type n, const bool nanflag, const bool realabs) \
1688 { \
1689 if (! n) \
1690 return; \
1691 if (! nanflag) \
1692 { \
1693 for (octave_idx_type i = 0; i < n; i++) \
1694 if (octave::math::isnan (v[i])) \
1695 { \
1696 *r = NAN; \
1697 return; \
1698 } \
1699 } \
1700 std::complex<T> tmp = v[0]; \
1701 octave_idx_type i = 1; \
1702 if (octave::math::isnan (tmp)) \
1703 { \
1704 for (; i < n && octave::math::isnan (v[i]); i++) ; \
1705 if (i < n) \
1706 tmp = v[i]; \
1707 } \
1708 if (realabs) \
1709 { \
1710 for (; i < n; i++) \
1711 { \
1712 if (v[i].real () OP tmp.real ()) \
1713 tmp = v[i]; \
1714 else if (v[i].real () == tmp.real () && \
1715 v[i].imag () OP tmp.imag ()) \
1716 tmp = v[i]; \
1717 } \
1718 } \
1719 else \
1720 { \
1721 for (; i < n; i++) \
1722 { \
1723 if (std::abs (v[i]) OP std::abs (tmp)) \
1724 tmp = v[i]; \
1725 else if (std::abs (v[i]) == std::abs (tmp) && \
1726 std::arg (v[i]) OP std::arg (tmp)) \
1727 tmp = v[i]; \
1728 } \
1729 } \
1730 *r = tmp; \
1731 } \
1732 template <typename T> \
1733 void F (const std::complex<T> *v, std::complex<T> *r, \
1734 octave_idx_type *ri, octave_idx_type n, const bool nanflag, \
1735 const bool realabs) \
1736 { \
1737 if (! n) \
1738 return; \
1739 if (! nanflag) \
1740 { \
1741 for (octave_idx_type i = 0; i < n; i++) \
1742 if (octave::math::isnan (v[i])) \
1743 { \
1744 *r = NAN; \
1745 *ri = i; \
1746 return; \
1747 } \
1748 } \
1749 std::complex<T> tmp = v[0]; \
1750 octave_idx_type tmpi = 0; \
1751 octave_idx_type i = 1; \
1752 if (octave::math::isnan (tmp)) \
1753 { \
1754 for (; i < n && octave::math::isnan (v[i]); i++) ; \
1755 if (i < n) \
1756 { \
1757 tmp = v[i]; \
1758 tmpi = i; \
1759 } \
1760 } \
1761 if (realabs) \
1762 { \
1763 for (; i < n; i++) \
1764 { \
1765 if (v[i].real () OP tmp.real ()) \
1766 { \
1767 tmp = v[i]; \
1768 tmpi = i; \
1769 } \
1770 else if (v[i].real () == tmp.real () && \
1771 v[i].imag () OP tmp.imag ()) \
1772 { \
1773 tmp = v[i]; \
1774 tmpi = i; \
1775 } \
1776 } \
1777 } \
1778 else \
1779 { \
1780 for (; i < n; i++) \
1781 { \
1782 if (std::abs (v[i]) OP std::abs (tmp)) \
1783 { \
1784 tmp = v[i]; \
1785 tmpi = i; \
1786 } \
1787 else if (std::abs (v[i]) == std::abs (tmp) && \
1788 std::arg (v[i]) OP std::arg (tmp)) \
1789 { \
1790 tmp = v[i]; \
1791 tmpi = i; \
1792 } \
1793 } \
1794 } \
1795 *r = tmp; \
1796 *ri = tmpi; \
1797 }
1798
1801
1802// Row reductions will be slightly complicated. We will proceed with checks
1803// for NaNs until we detect that no row will yield a NaN, in which case we
1804// proceed to a faster code.
1805
1806#define OP_CMINMAX_FCN2(F, OP) \
1807 template <typename T> \
1808 inline void \
1809 F (const std::complex<T> *v, std::complex<T> *r, octave_idx_type m, \
1810 octave_idx_type n, const bool nanflag, const bool realabs) \
1811 { \
1812 if (! n) \
1813 return; \
1814 octave_idx_type j = 0; \
1815 for (octave_idx_type i = 0; i < m; i++) \
1816 { \
1817 r[i] = v[i]; \
1818 } \
1819 j++; \
1820 v += m; \
1821 if (realabs) \
1822 { \
1823 while (j < n) \
1824 { \
1825 for (octave_idx_type i = 0; i < m; i++) \
1826 { \
1827 if ((octave::math::isnan (v[i]) || \
1828 octave::math::isnan (r[i])) && ! nanflag) \
1829 r[i] = NAN; \
1830 else if (octave::math::isnan (v[i])); \
1831 else if (octave::math::isnan (r[i]) || \
1832 v[i].real () OP r[i].real ()) \
1833 r[i] = v[i]; \
1834 else if (v[i].real () == r[i].real () && \
1835 v[i].imag () OP r[i].imag ()) \
1836 r[i] = v[i]; \
1837 } \
1838 j++; \
1839 v += m; \
1840 } \
1841 } \
1842 else \
1843 { \
1844 while (j < n) \
1845 { \
1846 for (octave_idx_type i = 0; i < m; i++) \
1847 { \
1848 if ((octave::math::isnan (v[i]) || \
1849 octave::math::isnan (r[i])) && ! nanflag) \
1850 r[i] = NAN; \
1851 else if (octave::math::isnan (v[i])); \
1852 else if (octave::math::isnan (r[i]) || \
1853 std::abs (v[i]) OP std::abs (r[i])) \
1854 r[i] = v[i]; \
1855 else if (std::abs (v[i]) == std::abs (r[i]) && \
1856 std::arg (v[i]) OP std::arg (r[i])) \
1857 r[i] = v[i]; \
1858 } \
1859 j++; \
1860 v += m; \
1861 } \
1862 } \
1863 } \
1864 template <typename T> \
1865 inline void \
1866 F (const std::complex<T> *v, std::complex<T> *r, octave_idx_type *ri, \
1867 octave_idx_type m, octave_idx_type n, const bool nanflag, \
1868 const bool realabs) \
1869 { \
1870 if (! n) \
1871 return; \
1872 octave_idx_type j = 0; \
1873 for (octave_idx_type i = 0; i < m; i++) \
1874 { \
1875 r[i] = v[i]; \
1876 ri[i] = j; \
1877 } \
1878 j++; \
1879 v += m; \
1880 if (realabs) \
1881 { \
1882 while (j < n) \
1883 { \
1884 for (octave_idx_type i = 0; i < m; i++) \
1885 { \
1886 if (octave::math::isnan (r[i]) && ! nanflag); \
1887 else if (octave::math::isnan (v[i]) && ! nanflag) \
1888 { \
1889 r[i] = NAN; \
1890 ri[i] = j; \
1891 } \
1892 else if (octave::math::isnan (v[i])); \
1893 else if (octave::math::isnan (r[i]) || \
1894 v[i].real () OP r[i].real ()) \
1895 { \
1896 r[i] = v[i]; \
1897 ri[i] = j; \
1898 } \
1899 else if (v[i].real () == r[i].real () && \
1900 v[i].imag () OP r[i].imag ()) \
1901 { \
1902 r[i] = v[i]; \
1903 ri[i] = j; \
1904 } \
1905 } \
1906 j++; \
1907 v += m; \
1908 } \
1909 } \
1910 else \
1911 { \
1912 while (j < n) \
1913 { \
1914 for (octave_idx_type i = 0; i < m; i++) \
1915 { \
1916 if (octave::math::isnan (r[i]) && ! nanflag); \
1917 else if (octave::math::isnan (v[i]) && ! nanflag) \
1918 { \
1919 r[i] = NAN; \
1920 ri[i] = j; \
1921 } \
1922 else if (octave::math::isnan (v[i])); \
1923 else if (octave::math::isnan (r[i]) || \
1924 std::abs (v[i]) OP std::abs (r[i])) \
1925 { \
1926 r[i] = v[i]; \
1927 ri[i] = j; \
1928 } \
1929 else if (std::abs (v[i]) == std::abs (r[i]) && \
1930 std::arg (v[i]) OP std::arg (r[i])) \
1931 { \
1932 r[i] = v[i]; \
1933 ri[i] = j; \
1934 } \
1935 } \
1936 j++; \
1937 v += m; \
1938 } \
1939 } \
1940 }
1941
1944
1945#define OP_CMINMAX_FCNN(F) \
1946 template <typename T> \
1947 inline void \
1948 F (const T *v, T *r, octave_idx_type l, octave_idx_type n, \
1949 octave_idx_type u, const bool nanflag, const bool realabs) \
1950 { \
1951 if (! n) \
1952 return; \
1953 if (l == 1) \
1954 { \
1955 for (octave_idx_type i = 0; i < u; i++) \
1956 { \
1957 F (v, r, n, nanflag, realabs); \
1958 v += n; \
1959 r++; \
1960 } \
1961 } \
1962 else \
1963 { \
1964 for (octave_idx_type i = 0; i < u; i++) \
1965 { \
1966 F (v, r, l, n, nanflag, realabs); \
1967 v += l*n; \
1968 r += l; \
1969 } \
1970 } \
1971 } \
1972 template <typename T> \
1973 inline void \
1974 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type l, \
1975 octave_idx_type n, octave_idx_type u, bool nanflag, bool realabs) \
1976 { \
1977 if (! n) return; \
1978 if (l == 1) \
1979 { \
1980 for (octave_idx_type i = 0; i < u; i++) \
1981 { \
1982 F (v, r, ri, n, nanflag, realabs); \
1983 v += n; \
1984 r++; \
1985 ri++; \
1986 } \
1987 } \
1988 else \
1989 { \
1990 for (octave_idx_type i = 0; i < u; i++) \
1991 { \
1992 F (v, r, ri, l, n, nanflag, realabs); \
1993 v += l*n; \
1994 r += l; \
1995 ri += l; \
1996 } \
1997 } \
1998 }
1999
2002
2003// Faster implementation for int arrays, since they do not support
2004// NaN values and real/abs comparison methods require specialized mapper.
2005
2006#define OP_INT_MINMAX_FCN(F, OP) \
2007 template <typename T> \
2008 void F (const T *v, T *r, octave_idx_type n, const bool realabs) \
2009 { \
2010 if (! n) \
2011 return; \
2012 T tmp = v[0]; \
2013 octave_idx_type i = 1; \
2014 if (realabs) \
2015 { \
2016 for (; i < n; i++) \
2017 if (v[i] OP tmp) \
2018 tmp = v[i]; \
2019 } \
2020 else \
2021 { \
2022 for (; i < n; i++) \
2023 { \
2024 if (mappers_abs (v[i]) OP mappers_abs (tmp)) \
2025 tmp = v[i]; \
2026 else if (mappers_abs (v[i]) == mappers_abs (tmp) && \
2027 v[i] OP tmp) \
2028 tmp = v[i]; \
2029 } \
2030 } \
2031 *r = tmp; \
2032 } \
2033 template <typename T> \
2034 void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n, \
2035 const bool realabs) \
2036 { \
2037 if (! n) \
2038 return; \
2039 T tmp = v[0]; \
2040 octave_idx_type tmpi = 0; \
2041 octave_idx_type i = 1; \
2042 if (realabs) \
2043 { \
2044 for (; i < n; i++) \
2045 if (v[i] OP tmp) \
2046 { \
2047 tmp = v[i]; \
2048 tmpi = i; \
2049 } \
2050 } \
2051 else \
2052 { \
2053 for (; i < n; i++) \
2054 { \
2055 if (mappers_abs (v[i]) OP mappers_abs (tmp)) \
2056 { \
2057 tmp = v[i]; \
2058 tmpi = i; \
2059 } \
2060 else if (mappers_abs (v[i]) == mappers_abs (tmp) && \
2061 v[i] OP tmp) \
2062 { \
2063 tmp = v[i]; \
2064 tmpi = i; \
2065 } \
2066 } \
2067 } \
2068 *r = tmp; \
2069 *ri = tmpi; \
2070 }
2071
2074
2075#define OP_INT_MINMAX_FCN2(F, OP) \
2076 template <typename T> \
2077 inline void \
2078 F (const T *v, T *r, octave_idx_type m, octave_idx_type n, \
2079 const bool realabs) \
2080 { \
2081 if (! n) \
2082 return; \
2083 octave_idx_type j = 0; \
2084 for (octave_idx_type i = 0; i < m; i++) \
2085 { \
2086 r[i] = v[i]; \
2087 } \
2088 j++; \
2089 v += m; \
2090 if (realabs) \
2091 { \
2092 while (j < n) \
2093 { \
2094 for (octave_idx_type i = 0; i < m; i++) \
2095 if (v[i] OP r[i]) \
2096 r[i] = v[i]; \
2097 j++; \
2098 v += m; \
2099 } \
2100 } \
2101 else \
2102 { \
2103 while (j < n) \
2104 { \
2105 for (octave_idx_type i = 0; i < m; i++) \
2106 { \
2107 if (mappers_abs (v[i]) OP mappers_abs (r[i])) \
2108 r[i] = v[i]; \
2109 else if (mappers_abs (v[i]) == mappers_abs (r[i]) && \
2110 v[i] OP r[i]) \
2111 r[i] = v[i]; \
2112 } \
2113 j++; \
2114 v += m; \
2115 } \
2116 } \
2117 } \
2118 template <typename T> \
2119 inline void \
2120 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type m, \
2121 octave_idx_type n, const bool realabs) \
2122 { \
2123 if (! n) \
2124 return; \
2125 octave_idx_type j = 0; \
2126 for (octave_idx_type i = 0; i < m; i++) \
2127 { \
2128 r[i] = v[i]; \
2129 ri[i] = j; \
2130 } \
2131 j++; \
2132 v += m; \
2133 if (realabs) \
2134 { \
2135 while (j < n) \
2136 { \
2137 for (octave_idx_type i = 0; i < m; i++) \
2138 if (v[i] OP r[i]) \
2139 { \
2140 r[i] = v[i]; \
2141 ri[i] = j; \
2142 } \
2143 j++; \
2144 v += m; \
2145 } \
2146 } \
2147 else \
2148 { \
2149 while (j < n) \
2150 { \
2151 for (octave_idx_type i = 0; i < m; i++) \
2152 { \
2153 if (mappers_abs (v[i]) OP mappers_abs (r[i])) \
2154 { \
2155 r[i] = v[i]; \
2156 ri[i] = j; \
2157 } \
2158 else if (mappers_abs (v[i]) == mappers_abs (r[i]) && \
2159 v[i] OP r[i]) \
2160 { \
2161 r[i] = v[i]; \
2162 ri[i] = j; \
2163 } \
2164 } \
2165 j++; \
2166 v += m; \
2167 } \
2168 } \
2169 }
2170
2173
2174#define OP_INT_MINMAX_FCNN(F) \
2175 template <typename T> \
2176 inline void \
2177 F (const T *v, T *r, octave_idx_type l, octave_idx_type n, \
2178 octave_idx_type u, const bool realabs) \
2179 { \
2180 if (! n) \
2181 return; \
2182 if (l == 1) \
2183 { \
2184 for (octave_idx_type i = 0; i < u; i++) \
2185 { \
2186 F (v, r, n, realabs); \
2187 v += n; \
2188 r++; \
2189 } \
2190 } \
2191 else \
2192 { \
2193 for (octave_idx_type i = 0; i < u; i++) \
2194 { \
2195 F (v, r, l, n, realabs); \
2196 v += l*n; \
2197 r += l; \
2198 } \
2199 } \
2200 } \
2201 template <typename T> \
2202 inline void \
2203 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type l, \
2204 octave_idx_type n, octave_idx_type u, const bool realabs) \
2205 { \
2206 if (! n) return; \
2207 if (l == 1) \
2208 { \
2209 for (octave_idx_type i = 0; i < u; i++) \
2210 { \
2211 F (v, r, ri, n, realabs); \
2212 v += n; \
2213 r++; \
2214 ri++; \
2215 } \
2216 } \
2217 else \
2218 { \
2219 for (octave_idx_type i = 0; i < u; i++) \
2220 { \
2221 F (v, r, ri, l, n, realabs); \
2222 v += l*n; \
2223 r += l; \
2224 ri += l; \
2225 } \
2226 } \
2227 }
2228
2231
2232// Faster implementation for char arrays, since they do not support
2233// NaN values and real/abs comparison methods do apply either.
2234
2235#define OP_CHMINMAX_FCN(F, OP) \
2236 template <typename T> \
2237 void F (const T *v, T *r, octave_idx_type n) \
2238 { \
2239 if (! n) \
2240 return; \
2241 T tmp = v[0]; \
2242 octave_idx_type i = 1; \
2243 for (; i < n; i++) \
2244 if (v[i] OP tmp) \
2245 tmp = v[i]; \
2246 *r = tmp; \
2247 } \
2248 template <typename T> \
2249 void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n) \
2250 { \
2251 if (! n) \
2252 return; \
2253 T tmp = v[0]; \
2254 octave_idx_type tmpi = 0; \
2255 octave_idx_type i = 1; \
2256 for (; i < n; i++) \
2257 if (v[i] OP tmp) \
2258 { \
2259 tmp = v[i]; \
2260 tmpi = i; \
2261 } \
2262 *r = tmp; \
2263 *ri = tmpi; \
2264 }
2265
2268
2269#define OP_CHMINMAX_FCN2(F, OP) \
2270 template <typename T> \
2271 inline void \
2272 F (const T *v, T *r, octave_idx_type m, octave_idx_type n) \
2273 { \
2274 if (! n) \
2275 return; \
2276 octave_idx_type j = 1; \
2277 for (octave_idx_type i = 0; i < m; i++) \
2278 r[i] = v[i]; \
2279 v += m; \
2280 while (j < n) \
2281 { \
2282 for (octave_idx_type i = 0; i < m; i++) \
2283 if (v[i] OP r[i]) \
2284 r[i] = v[i]; \
2285 j++; \
2286 v += m; \
2287 } \
2288 } \
2289 template <typename T> \
2290 inline void \
2291 F (const T *v, T *r, octave_idx_type *ri, \
2292 octave_idx_type m, octave_idx_type n) \
2293 { \
2294 if (! n) \
2295 return; \
2296 octave_idx_type j = 1; \
2297 for (octave_idx_type i = 0; i < m; i++) \
2298 { \
2299 r[i] = v[i]; \
2300 ri[i] = j; \
2301 } \
2302 v += m; \
2303 while (j < n) \
2304 { \
2305 for (octave_idx_type i = 0; i < m; i++) \
2306 if (v[i] OP r[i]) \
2307 { \
2308 r[i] = v[i]; \
2309 ri[i] = j; \
2310 } \
2311 j++; \
2312 v += m; \
2313 } \
2314 }
2315
2318
2319#define OP_CHMINMAX_FCNN(F) \
2320 template <typename T> \
2321 inline void \
2322 F (const T *v, T *r, octave_idx_type l, \
2323 octave_idx_type n,octave_idx_type u) \
2324 { \
2325 if (! n) \
2326 return; \
2327 if (l == 1) \
2328 { \
2329 for (octave_idx_type i = 0; i < u; i++) \
2330 { \
2331 F (v, r, n); \
2332 v += n; \
2333 r++; \
2334 } \
2335 } \
2336 else \
2337 { \
2338 for (octave_idx_type i = 0; i < u; i++) \
2339 { \
2340 F (v, r, l, n); \
2341 v += l*n; \
2342 r += l; \
2343 } \
2344 } \
2345 } \
2346 template <typename T> \
2347 inline void \
2348 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type l, \
2349 octave_idx_type n, octave_idx_type u) \
2350 { \
2351 if (! n) return; \
2352 if (l == 1) \
2353 { \
2354 for (octave_idx_type i = 0; i < u; i++) \
2355 { \
2356 F (v, r, ri, n); \
2357 v += n; \
2358 r++; \
2359 ri++; \
2360 } \
2361 } \
2362 else \
2363 { \
2364 for (octave_idx_type i = 0; i < u; i++) \
2365 { \
2366 F (v, r, ri, l, n); \
2367 v += l*n; \
2368 r += l; \
2369 ri += l; \
2370 } \
2371 } \
2372 }
2373
2376
2377// Templated functions for Fcummin and Fcummax
2378
2379#define OP_CUMMINMAX_FCN(F, OP) \
2380 template <typename T> \
2381 void F (const T *v, T *r, octave_idx_type n, const bool nanflag, \
2382 const bool realabs) \
2383 { \
2384 if (! n) \
2385 return; \
2386 T tmp = v[0]; \
2387 octave_idx_type i = 1; \
2388 octave_idx_type j = 0; \
2389 if (octave::math::isnan (tmp) && ! nanflag) \
2390 { \
2391 for (; j < n; j++) \
2392 r[j] = NAN; \
2393 return; \
2394 } \
2395 else if (octave::math::isnan (tmp)) \
2396 { \
2397 for (; i < n && octave::math::isnan (v[i]); i++) ; \
2398 for (; j < i; j++) \
2399 r[j] = tmp; \
2400 if (i < n) \
2401 tmp = v[i]; \
2402 } \
2403 if (realabs) \
2404 { \
2405 for (; i < n; i++) \
2406 { \
2407 if (octave::math::isnan (v[i]) && ! nanflag) \
2408 { \
2409 for (; j < i; j++) \
2410 r[j] = tmp; \
2411 for (; j < n; j++) \
2412 r[j] = NAN; \
2413 return; \
2414 } \
2415 if (v[i] OP tmp) \
2416 { \
2417 for (; j < i; j++) \
2418 r[j] = tmp; \
2419 tmp = v[i]; \
2420 } \
2421 } \
2422 } \
2423 else \
2424 { \
2425 for (; i < n; i++) \
2426 { \
2427 if (octave::math::isnan (v[i]) && ! nanflag) \
2428 { \
2429 for (; j < i; j++) \
2430 r[j] = tmp; \
2431 for (; j < n; j++) \
2432 r[j] = NAN; \
2433 return; \
2434 } \
2435 if (std::abs (v[i]) OP std::abs (tmp)) \
2436 { \
2437 for (; j < i; j++) \
2438 r[j] = tmp; \
2439 tmp = v[i]; \
2440 } \
2441 else if (std::abs (v[i]) == std::abs (tmp) && v[i] OP tmp) \
2442 { \
2443 for (; j < i; j++) \
2444 r[j] = tmp; \
2445 tmp = v[i]; \
2446 } \
2447 } \
2448 } \
2449 for (; j < i; j++) \
2450 r[j] = tmp; \
2451 } \
2452 template <typename T> \
2453 void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n, \
2454 const bool nanflag, const bool realabs) \
2455 { \
2456 if (! n) \
2457 return; \
2458 T tmp = v[0]; \
2459 octave_idx_type tmpi = 0; \
2460 octave_idx_type i = 1; \
2461 octave_idx_type j = 0; \
2462 if (octave::math::isnan (tmp) && ! nanflag) \
2463 { \
2464 for (; j < n; j++) \
2465 { \
2466 r[j] = NAN; \
2467 ri[j] = tmpi; \
2468 } \
2469 return; \
2470 } \
2471 else if (octave::math::isnan (tmp)) \
2472 { \
2473 for (; i < n && octave::math::isnan (v[i]); i++) ; \
2474 for (; j < i; j++) \
2475 { \
2476 r[j] = tmp; \
2477 ri[j] = tmpi; \
2478 } \
2479 if (i < n) \
2480 { \
2481 tmp = v[i]; \
2482 tmpi = i; \
2483 } \
2484 } \
2485 if (realabs) \
2486 { \
2487 for (; i < n; i++) \
2488 { \
2489 if (octave::math::isnan (v[i]) && ! nanflag) \
2490 { \
2491 for (; j < i; j++) \
2492 { \
2493 r[j] = tmp; \
2494 ri[j] = tmpi; \
2495 } \
2496 for (; j < n; j++) \
2497 { \
2498 r[j] = NAN; \
2499 ri[j] = i; \
2500 } \
2501 return; \
2502 } \
2503 if (v[i] OP tmp) \
2504 { \
2505 for (; j < i; j++) \
2506 { \
2507 r[j] = tmp; \
2508 ri[j] = tmpi; \
2509 } \
2510 tmp = v[i]; \
2511 tmpi = i; \
2512 } \
2513 } \
2514 } \
2515 else \
2516 { \
2517 for (; i < n; i++) \
2518 { \
2519 if (octave::math::isnan (v[i]) && ! nanflag) \
2520 { \
2521 for (; j < i; j++) \
2522 { \
2523 r[j] = tmp; \
2524 ri[j] = tmpi; \
2525 } \
2526 for (; j < n; j++) \
2527 { \
2528 r[j] = NAN; \
2529 ri[j] = i; \
2530 } \
2531 return; \
2532 } \
2533 if (std::abs (v[i]) OP std::abs (tmp)) \
2534 { \
2535 for (; j < i; j++) \
2536 { \
2537 r[j] = tmp; \
2538 ri[j] = tmpi; \
2539 } \
2540 tmp = v[i]; \
2541 tmpi = i; \
2542 } \
2543 else if (std::abs (v[i]) == std::abs (tmp) && v[i] OP tmp) \
2544 { \
2545 for (; j < i; j++) \
2546 { \
2547 r[j] = tmp; \
2548 ri[j] = tmpi; \
2549 } \
2550 tmp = v[i]; \
2551 tmpi = i; \
2552 } \
2553 } \
2554 } \
2555 for (; j < i; j++) \
2556 { \
2557 r[j] = tmp; \
2558 ri[j] = tmpi; \
2559 } \
2560 }
2561
2564
2565#define OP_CUMMINMAX_FCN2(F, OP) \
2566 template <typename T> \
2567 inline void \
2568 F (const T *v, T *r, octave_idx_type m, octave_idx_type n, \
2569 const bool nanflag, const bool realabs) \
2570 { \
2571 if (! n) \
2572 return; \
2573 const T *r0; \
2574 octave_idx_type j = 0; \
2575 for (octave_idx_type i = 0; i < m; i++) \
2576 r[i] = v[i]; \
2577 j++; \
2578 v += m; \
2579 r0 = r; \
2580 r += m; \
2581 if (realabs) \
2582 { \
2583 while (! nanflag && j < n) \
2584 { \
2585 for (octave_idx_type i = 0; i < m; i++) \
2586 { \
2587 if (octave::math::isnan (v[i]) || \
2588 octave::math::isnan (r0[i])) \
2589 r[i] = NAN; \
2590 else if (v[i] OP r0[i]) \
2591 r[i] = v[i]; \
2592 else \
2593 r[i] = r0[i]; \
2594 } \
2595 j++; \
2596 v += m; \
2597 r0 = r; \
2598 r += m; \
2599 } \
2600 while (nanflag && j < n) \
2601 { \
2602 for (octave_idx_type i = 0; i < m; i++) \
2603 { \
2604 if (octave::math::isnan (r0[i]) || v[i] OP r0[i]) \
2605 r[i] = v[i]; \
2606 else \
2607 r[i] = r0[i]; \
2608 } \
2609 j++; \
2610 v += m; \
2611 r0 = r; \
2612 r += m; \
2613 } \
2614 } \
2615 else \
2616 { \
2617 while (! nanflag && j < n) \
2618 { \
2619 for (octave_idx_type i = 0; i < m; i++) \
2620 { \
2621 if (octave::math::isnan (v[i]) || \
2622 octave::math::isnan (r0[i])) \
2623 r[i] = NAN; \
2624 else if (std::abs (v[i]) OP std::abs (r0[i])) \
2625 r[i] = v[i]; \
2626 else if (std::abs (v[i]) == std::abs (r0[i]) && \
2627 v[i] OP r0[i]) \
2628 r[i] = v[i]; \
2629 else \
2630 r[i] = r0[i]; \
2631 } \
2632 j++; \
2633 v += m; \
2634 r0 = r; \
2635 r += m; \
2636 } \
2637 while (nanflag && j < n) \
2638 { \
2639 for (octave_idx_type i = 0; i < m; i++) \
2640 { \
2641 if (octave::math::isnan (r0[i])) \
2642 r[i] = v[i]; \
2643 else if (std::abs (v[i]) OP std::abs (r0[i])) \
2644 r[i] = v[i]; \
2645 else if (std::abs (v[i]) == std::abs (r0[i]) && \
2646 v[i] OP r0[i]) \
2647 r[i] = v[i]; \
2648 else \
2649 r[i] = r0[i]; \
2650 } \
2651 j++; \
2652 v += m; \
2653 r0 = r; \
2654 r += m; \
2655 } \
2656 } \
2657 } \
2658 template <typename T> \
2659 inline void \
2660 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type m, \
2661 octave_idx_type n, const bool nanflag, const bool realabs) \
2662 { \
2663 if (! n) \
2664 return; \
2665 const T *r0; \
2666 const octave_idx_type *r0i; \
2667 octave_idx_type j = 0; \
2668 for (octave_idx_type i = 0; i < m; i++) \
2669 { \
2670 r[i] = v[i]; \
2671 ri[i] = 0; \
2672 } \
2673 j++; \
2674 v += m; \
2675 r0 = r; \
2676 r += m; \
2677 r0i = ri; \
2678 ri += m; \
2679 if (realabs) \
2680 { \
2681 while (! nanflag && j < n) \
2682 { \
2683 for (octave_idx_type i = 0; i < m; i++) \
2684 { \
2685 if (octave::math::isnan (v[i]) && \
2686 octave::math::isnan (r0[i])) \
2687 { \
2688 r[i] = r0[i]; \
2689 ri[i] = r0i[i]; \
2690 } \
2691 else if (octave::math::isnan (v[i])) \
2692 { \
2693 r[i] = NAN; \
2694 ri[i] = j; \
2695 } \
2696 else if (v[i] OP r0[i]) \
2697 { \
2698 r[i] = v[i]; \
2699 ri[i] = j; \
2700 } \
2701 else \
2702 { \
2703 r[i] = r0[i]; \
2704 ri[i] = r0i[i]; \
2705 } \
2706 } \
2707 j++; \
2708 v += m; \
2709 r0 = r; \
2710 r += m; \
2711 r0i = ri; \
2712 ri += m; \
2713 } \
2714 while (nanflag && j < n) \
2715 { \
2716 for (octave_idx_type i = 0; i < m; i++) \
2717 { \
2718 if (octave::math::isnan (r0[i]) && \
2719 octave::math::isnan (v[i])) \
2720 { \
2721 r[i] = NAN; \
2722 ri[i] = r0i[i]; \
2723 } \
2724 else if (octave::math::isnan (r0[i]) || v[i] OP r0[i]) \
2725 { \
2726 r[i] = v[i]; \
2727 ri[i] = j; \
2728 } \
2729 else \
2730 { \
2731 r[i] = r0[i]; \
2732 ri[i] = r0i[i]; \
2733 } \
2734 } \
2735 j++; \
2736 v += m; \
2737 r0 = r; \
2738 r += m; \
2739 r0i = ri; \
2740 ri += m; \
2741 } \
2742 } \
2743 else \
2744 { \
2745 while (! nanflag && j < n) \
2746 { \
2747 for (octave_idx_type i = 0; i < m; i++) \
2748 { \
2749 if (octave::math::isnan (v[i]) && \
2750 octave::math::isnan (r0[i])) \
2751 { \
2752 r[i] = r0[i]; \
2753 ri[i] = r0i[i]; \
2754 } \
2755 else if (octave::math::isnan (v[i])) \
2756 { \
2757 r[i] = NAN; \
2758 ri[i] = j; \
2759 } \
2760 else if (std::abs (v[i]) OP std::abs (r0[i])) \
2761 { \
2762 r[i] = v[i]; \
2763 ri[i] = j; \
2764 } \
2765 else if (std::abs (v[i]) == std::abs (r0[i]) && \
2766 v[i] OP r0[i]) \
2767 { \
2768 r[i] = v[i]; \
2769 ri[i] = j; \
2770 } \
2771 else \
2772 { \
2773 r[i] = r0[i]; \
2774 ri[i] = r0i[i]; \
2775 } \
2776 } \
2777 j++; \
2778 v += m; \
2779 r0 = r; \
2780 r += m; \
2781 r0i = ri; \
2782 ri += m; \
2783 } \
2784 while (nanflag && j < n) \
2785 { \
2786 for (octave_idx_type i = 0; i < m; i++) \
2787 { \
2788 if (octave::math::isnan (r0[i]) && \
2789 octave::math::isnan (v[i])) \
2790 { \
2791 r[i] = NAN; \
2792 ri[i] = r0i[i]; \
2793 } \
2794 if (octave::math::isnan (r0[i])) \
2795 { \
2796 r[i] = v[i]; \
2797 ri[i] = j; \
2798 } \
2799 else if (std::abs (v[i]) OP std::abs (r0[i])) \
2800 { \
2801 r[i] = v[i]; \
2802 ri[i] = j; \
2803 } \
2804 else if (std::abs (v[i]) == std::abs (r0[i]) && \
2805 v[i] OP r0[i]) \
2806 { \
2807 r[i] = v[i]; \
2808 ri[i] = j; \
2809 } \
2810 else \
2811 { \
2812 r[i] = r0[i]; \
2813 ri[i] = r0i[i]; \
2814 } \
2815 } \
2816 j++; \
2817 v += m; \
2818 r0 = r; \
2819 r += m; \
2820 r0i = ri; \
2821 ri += m; \
2822 } \
2823 } \
2824 }
2825
2828
2829#define OP_CUMMINMAX_FCNN(F) \
2830 template <typename T> \
2831 inline void \
2832 F (const T *v, T *r, octave_idx_type l, octave_idx_type n, \
2833 octave_idx_type u, const bool nanflag, const bool realabs) \
2834 { \
2835 if (! n) \
2836 return; \
2837 if (l == 1) \
2838 { \
2839 for (octave_idx_type i = 0; i < u; i++) \
2840 { \
2841 F (v, r, n, nanflag, realabs); \
2842 v += n; \
2843 r += n; \
2844 } \
2845 } \
2846 else \
2847 { \
2848 for (octave_idx_type i = 0; i < u; i++) \
2849 { \
2850 F (v, r, l, n, nanflag, realabs); \
2851 v += l*n; \
2852 r += l*n; \
2853 } \
2854 } \
2855 } \
2856 template <typename T> \
2857 inline void \
2858 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type l, \
2859 octave_idx_type n, octave_idx_type u, const bool nanflag, \
2860 const bool realabs) \
2861 { \
2862 if (! n) \
2863 return; \
2864 if (l == 1) \
2865 { \
2866 for (octave_idx_type i = 0; i < u; i++) \
2867 { \
2868 F (v, r, ri, n, nanflag, realabs); \
2869 v += n; \
2870 r += n; \
2871 ri += n; \
2872 } \
2873 } \
2874 else \
2875 { \
2876 for (octave_idx_type i = 0; i < u; i++) \
2877 { \
2878 F (v, r, ri, l, n, nanflag, realabs); \
2879 v += l*n; \
2880 r += l*n; \
2881 ri += l*n; \
2882 } \
2883 } \
2884 }
2885
2888
2889// Special implementation for complex arrays, since we need comparisons
2890// in both real and imaginary parts for MATLAB compatibility.
2891
2892#define OP_CCUMMINMAX_FCN(F, OP) \
2893 template <typename T> \
2894 void F (const std::complex<T> *v, std::complex<T> *r, \
2895 octave_idx_type n, const bool nanflag, const bool realabs) \
2896 { \
2897 if (! n) \
2898 return; \
2899 std::complex<T> tmp = v[0]; \
2900 octave_idx_type i = 1; \
2901 octave_idx_type j = 0; \
2902 if (octave::math::isnan (tmp) && ! nanflag) \
2903 { \
2904 for (; j < n; j++) \
2905 r[j] = NAN; \
2906 return; \
2907 } \
2908 else if (octave::math::isnan (tmp)) \
2909 { \
2910 for (; i < n && octave::math::isnan (v[i]); i++) ; \
2911 for (; j < i; j++) \
2912 r[j] = tmp; \
2913 if (i < n) \
2914 tmp = v[i]; \
2915 } \
2916 if (realabs) \
2917 { \
2918 for (; i < n; i++) \
2919 { \
2920 if (octave::math::isnan (v[i]) && ! nanflag) \
2921 { \
2922 for (; j < i; j++) \
2923 r[j] = tmp; \
2924 for (; j < n; j++) \
2925 r[j] = NAN; \
2926 return; \
2927 } \
2928 if (v[i].real () OP tmp.real ()) \
2929 { \
2930 for (; j < i; j++) \
2931 r[j] = tmp; \
2932 tmp = v[i]; \
2933 } \
2934 else if (v[i].real () == tmp.real () && \
2935 v[i].imag () OP tmp.imag ()) \
2936 { \
2937 for (; j < i; j++) \
2938 r[j] = tmp; \
2939 tmp = v[i]; \
2940 } \
2941 } \
2942 } \
2943 else \
2944 { \
2945 for (; i < n; i++) \
2946 { \
2947 if (octave::math::isnan (v[i]) && ! nanflag) \
2948 { \
2949 for (; j < i; j++) \
2950 r[j] = tmp; \
2951 for (; j < n; j++) \
2952 r[j] = NAN; \
2953 return; \
2954 } \
2955 if (std::abs (v[i]) OP std::abs (tmp)) \
2956 { \
2957 for (; j < i; j++) \
2958 r[j] = tmp; \
2959 tmp = v[i]; \
2960 } \
2961 else if (std::abs (v[i]) == std::abs (tmp) && \
2962 std::arg (v[i]) OP std::arg (tmp)) \
2963 { \
2964 for (; j < i; j++) \
2965 r[j] = tmp; \
2966 tmp = v[i]; \
2967 } \
2968 } \
2969 } \
2970 for (; j < i; j++) \
2971 r[j] = tmp; \
2972 } \
2973 template <typename T> \
2974 void F (const std::complex<T> *v, std::complex<T> *r, \
2975 octave_idx_type *ri, octave_idx_type n, const bool nanflag, \
2976 const bool realabs) \
2977 { \
2978 if (! n) \
2979 return; \
2980 std::complex<T> tmp = v[0]; \
2981 octave_idx_type tmpi = 0; \
2982 octave_idx_type i = 1; \
2983 octave_idx_type j = 0; \
2984 if (octave::math::isnan (tmp) && ! nanflag) \
2985 { \
2986 for (; j < n; j++) \
2987 { \
2988 r[j] = NAN; \
2989 ri[j] = tmpi; \
2990 } \
2991 return; \
2992 } \
2993 else if (octave::math::isnan (tmp)) \
2994 { \
2995 for (; i < n && octave::math::isnan (v[i]); i++) ; \
2996 for (; j < i; j++) \
2997 { \
2998 r[j] = tmp; \
2999 ri[j] = tmpi; \
3000 } \
3001 if (i < n) \
3002 { \
3003 tmp = v[i]; \
3004 tmpi = i; \
3005 } \
3006 } \
3007 if (realabs) \
3008 { \
3009 for (; i < n; i++) \
3010 { \
3011 if (octave::math::isnan (v[i]) && ! nanflag) \
3012 { \
3013 for (; j < i; j++) \
3014 { \
3015 r[j] = tmp; \
3016 ri[j] = tmpi; \
3017 } \
3018 for (; j < n; j++) \
3019 { \
3020 r[j] = NAN; \
3021 ri[j] = i; \
3022 } \
3023 return; \
3024 } \
3025 if (v[i].real () OP tmp.real ()) \
3026 { \
3027 for (; j < i; j++) \
3028 { \
3029 r[j] = tmp; \
3030 ri[j] = tmpi; \
3031 } \
3032 tmp = v[i]; \
3033 tmpi = i; \
3034 } \
3035 else if (v[i].real () == tmp.real () && \
3036 v[i].imag () OP tmp.imag ()) \
3037 { \
3038 for (; j < i; j++) \
3039 { \
3040 r[j] = tmp; \
3041 ri[j] = tmpi; \
3042 } \
3043 tmp = v[i]; \
3044 tmpi = i; \
3045 } \
3046 } \
3047 } \
3048 else \
3049 { \
3050 for (; i < n; i++) \
3051 { \
3052 if (octave::math::isnan (v[i]) && ! nanflag) \
3053 { \
3054 for (; j < i; j++) \
3055 { \
3056 r[j] = tmp; \
3057 ri[j] = tmpi; \
3058 } \
3059 for (; j < n; j++) \
3060 { \
3061 r[j] = NAN; \
3062 ri[j] = i; \
3063 } \
3064 return; \
3065 } \
3066 if (std::abs (v[i]) OP std::abs (tmp)) \
3067 { \
3068 for (; j < i; j++) \
3069 { \
3070 r[j] = tmp; \
3071 ri[j] = tmpi; \
3072 } \
3073 tmp = v[i]; \
3074 tmpi = i; \
3075 } \
3076 else if (std::abs (v[i]) == std::abs (tmp) && \
3077 std::arg (v[i]) OP std::arg (tmp)) \
3078 { \
3079 for (; j < i; j++) \
3080 { \
3081 r[j] = tmp; \
3082 ri[j] = tmpi; \
3083 } \
3084 tmp = v[i]; \
3085 tmpi = i; \
3086 } \
3087 } \
3088 } \
3089 for (; j < i; j++) \
3090 { \
3091 r[j] = tmp; \
3092 ri[j] = tmpi; \
3093 } \
3094 }
3095
3098
3099#define OP_CCUMMINMAX_FCN2(F, OP) \
3100 template <typename T> \
3101 inline void \
3102 F (const std::complex<T> *v, std::complex<T> *r, octave_idx_type m, \
3103 octave_idx_type n, const bool nanflag, const bool realabs) \
3104 { \
3105 if (! n) \
3106 return; \
3107 const std::complex<T> *r0; \
3108 octave_idx_type j = 0; \
3109 for (octave_idx_type i = 0; i < m; i++) \
3110 r[i] = v[i]; \
3111 j++; \
3112 v += m; \
3113 r0 = r; \
3114 r += m; \
3115 if (realabs) \
3116 { \
3117 while (! nanflag && j < n) \
3118 { \
3119 for (octave_idx_type i = 0; i < m; i++) \
3120 { \
3121 if (octave::math::isnan (v[i]) || \
3122 octave::math::isnan (r0[i])) \
3123 r[i] = NAN; \
3124 else if (v[i].real () OP r0[i].real ()) \
3125 r[i] = v[i]; \
3126 else if (v[i].real () == r0[i].real () && \
3127 v[i].imag () OP r0[i].imag ()) \
3128 r[i] = v[i]; \
3129 else \
3130 r[i] = r0[i]; \
3131 } \
3132 j++; \
3133 v += m; \
3134 r0 = r; \
3135 r += m; \
3136 } \
3137 while (nanflag && j < n) \
3138 { \
3139 for (octave_idx_type i = 0; i < m; i++) \
3140 { \
3141 if (octave::math::isnan (r0[i])) \
3142 r[i] = v[i]; \
3143 else if (v[i].real () OP r0[i].real ()) \
3144 r[i] = v[i]; \
3145 else if (v[i].real () == r0[i].real () && \
3146 v[i].imag () OP r0[i].imag ()) \
3147 r[i] = v[i]; \
3148 else \
3149 r[i] = r0[i]; \
3150 } \
3151 j++; \
3152 v += m; \
3153 r0 = r; \
3154 r += m; \
3155 } \
3156 } \
3157 else \
3158 { \
3159 while (! nanflag && j < n) \
3160 { \
3161 for (octave_idx_type i = 0; i < m; i++) \
3162 { \
3163 if (octave::math::isnan (v[i]) || \
3164 octave::math::isnan (r0[i])) \
3165 r[i] = NAN; \
3166 else if (std::abs (v[i]) OP std::abs (r0[i])) \
3167 r[i] = v[i]; \
3168 else if (std::abs (v[i]) == std::abs (r0[i]) && \
3169 std::arg (v[i]) OP std::arg (r0[i])) \
3170 r[i] = v[i]; \
3171 else \
3172 r[i] = r0[i]; \
3173 } \
3174 j++; \
3175 v += m; \
3176 r0 = r; \
3177 r += m; \
3178 } \
3179 while (nanflag && j < n) \
3180 { \
3181 for (octave_idx_type i = 0; i < m; i++) \
3182 { \
3183 if (octave::math::isnan (r0[i])) \
3184 r[i] = v[i]; \
3185 else if (std::abs (v[i]) OP std::abs (r0[i])) \
3186 r[i] = v[i]; \
3187 else if (std::abs (v[i]) == std::abs (r0[i]) && \
3188 std::arg (v[i]) OP std::arg (r0[i])) \
3189 r[i] = v[i]; \
3190 else \
3191 r[i] = r0[i]; \
3192 } \
3193 j++; \
3194 v += m; \
3195 r0 = r; \
3196 r += m; \
3197 } \
3198 } \
3199 } \
3200 template <typename T> \
3201 inline void \
3202 F (const std::complex<T> *v, std::complex<T> *r, octave_idx_type *ri, \
3203 octave_idx_type m, octave_idx_type n, const bool nanflag, \
3204 const bool realabs) \
3205 { \
3206 if (! n) \
3207 return; \
3208 const std::complex<T> *r0; \
3209 const octave_idx_type *r0i; \
3210 octave_idx_type j = 0; \
3211 for (octave_idx_type i = 0; i < m; i++) \
3212 { \
3213 r[i] = v[i]; \
3214 ri[i] = 0; \
3215 } \
3216 j++; \
3217 v += m; \
3218 r0 = r; \
3219 r += m; \
3220 r0i = ri; \
3221 ri += m; \
3222 if (realabs) \
3223 { \
3224 while (! nanflag && j < n) \
3225 { \
3226 for (octave_idx_type i = 0; i < m; i++) \
3227 { \
3228 if (octave::math::isnan (v[i]) && \
3229 octave::math::isnan (r0[i])) \
3230 { \
3231 r[i] = r0[i]; \
3232 ri[i] = r0i[i]; \
3233 } \
3234 else if (octave::math::isnan (v[i])) \
3235 { \
3236 r[i] = NAN; \
3237 ri[i] = j; \
3238 } \
3239 else if (v[i].real () OP r0[i].real ()) \
3240 { \
3241 r[i] = v[i]; \
3242 ri[i] = j; \
3243 } \
3244 else if (v[i].real () == r0[i].real () && \
3245 v[i].imag () OP r0[i].imag ()) \
3246 { \
3247 r[i] = v[i]; \
3248 ri[i] = j; \
3249 } \
3250 else \
3251 { \
3252 r[i] = r0[i]; \
3253 ri[i] = r0i[i]; \
3254 } \
3255 } \
3256 j++; \
3257 v += m; \
3258 r0 = r; \
3259 r += m; \
3260 r0i = ri; \
3261 ri += m; \
3262 } \
3263 while (nanflag && j < n) \
3264 { \
3265 for (octave_idx_type i = 0; i < m; i++) \
3266 { \
3267 if (octave::math::isnan (r0[i]) && \
3268 octave::math::isnan (v[i])) \
3269 { \
3270 r[i] = NAN; \
3271 ri[i] = r0i[i]; \
3272 } \
3273 else if (octave::math::isnan (r0[i])) \
3274 { \
3275 r[i] = v[i]; \
3276 ri[i] = j; \
3277 } \
3278 else if (v[i].real () OP r0[i].real ()) \
3279 { \
3280 r[i] = v[i]; \
3281 ri[i] = j; \
3282 } \
3283 else if (v[i].real () == r0[i].real () && \
3284 v[i].imag () OP r0[i].imag ()) \
3285 { \
3286 r[i] = v[i]; \
3287 ri[i] = j; \
3288 } \
3289 else \
3290 { \
3291 r[i] = r0[i]; \
3292 ri[i] = r0i[i]; \
3293 } \
3294 } \
3295 j++; \
3296 v += m; \
3297 r0 = r; \
3298 r += m; \
3299 r0i = ri; \
3300 ri += m; \
3301 } \
3302 } \
3303 else \
3304 { \
3305 while (! nanflag && j < n) \
3306 { \
3307 for (octave_idx_type i = 0; i < m; i++) \
3308 { \
3309 if (octave::math::isnan (v[i]) && \
3310 octave::math::isnan (r0[i])) \
3311 { \
3312 r[i] = r0[i]; \
3313 ri[i] = r0i[i]; \
3314 } \
3315 else if (octave::math::isnan (v[i])) \
3316 { \
3317 r[i] = NAN; \
3318 ri[i] = j; \
3319 } \
3320 else if (std::abs (v[i]) OP std::abs (r0[i])) \
3321 { \
3322 r[i] = v[i]; \
3323 ri[i] = j; \
3324 } \
3325 else if (std::abs (v[i]) == std::abs (r0[i]) && \
3326 std::arg (v[i]) OP std::arg (r0[i])) \
3327 { \
3328 r[i] = v[i]; \
3329 ri[i] = j; \
3330 } \
3331 else \
3332 { \
3333 r[i] = r0[i]; \
3334 ri[i] = r0i[i]; \
3335 } \
3336 } \
3337 j++; \
3338 v += m; \
3339 r0 = r; \
3340 r += m; \
3341 r0i = ri; \
3342 ri += m; \
3343 } \
3344 while (nanflag && j < n) \
3345 { \
3346 for (octave_idx_type i = 0; i < m; i++) \
3347 { \
3348 if (octave::math::isnan (r0[i]) && \
3349 octave::math::isnan (v[i])) \
3350 { \
3351 r[i] = NAN; \
3352 ri[i] = r0i[i]; \
3353 } \
3354 else if (octave::math::isnan (r0[i])) \
3355 { \
3356 r[i] = v[i]; \
3357 ri[i] = j; \
3358 } \
3359 else if (std::abs (v[i]) OP std::abs (r0[i])) \
3360 { \
3361 r[i] = v[i]; \
3362 ri[i] = j; \
3363 } \
3364 else if (std::abs (v[i]) == std::abs (r0[i]) && \
3365 std::arg (v[i]) OP std::arg (r0[i])) \
3366 { \
3367 r[i] = v[i]; \
3368 ri[i] = j; \
3369 } \
3370 else \
3371 { \
3372 r[i] = r0[i]; \
3373 ri[i] = r0i[i]; \
3374 } \
3375 } \
3376 j++; \
3377 v += m; \
3378 r0 = r; \
3379 r += m; \
3380 r0i = ri; \
3381 ri += m; \
3382 } \
3383 } \
3384 }
3385
3388
3389#define OP_CCUMMINMAX_FCNN(F) \
3390 template <typename T> \
3391 inline void \
3392 F (const T *v, T *r, octave_idx_type l, octave_idx_type n, \
3393 octave_idx_type u, const bool nanflag, const bool realabs) \
3394 { \
3395 if (! n) \
3396 return; \
3397 if (l == 1) \
3398 { \
3399 for (octave_idx_type i = 0; i < u; i++) \
3400 { \
3401 F (v, r, n, nanflag, realabs); \
3402 v += n; \
3403 r += n; \
3404 } \
3405 } \
3406 else \
3407 { \
3408 for (octave_idx_type i = 0; i < u; i++) \
3409 { \
3410 F (v, r, l, n, nanflag, realabs); \
3411 v += l*n; \
3412 r += l*n; \
3413 } \
3414 } \
3415 } \
3416 template <typename T> \
3417 inline void \
3418 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type l, \
3419 octave_idx_type n, octave_idx_type u, const bool nanflag, \
3420 const bool realabs) \
3421 { \
3422 if (! n) \
3423 return; \
3424 if (l == 1) \
3425 { \
3426 for (octave_idx_type i = 0; i < u; i++) \
3427 { \
3428 F (v, r, ri, n, nanflag, realabs); \
3429 v += n; \
3430 r += n; \
3431 ri += n; \
3432 } \
3433 } \
3434 else \
3435 { \
3436 for (octave_idx_type i = 0; i < u; i++) \
3437 { \
3438 F (v, r, ri, l, n, nanflag, realabs); \
3439 v += l*n; \
3440 r += l*n; \
3441 ri += l*n; \
3442 } \
3443 } \
3444 }
3445
3448
3449// Faster implementation for int arrays, since they do not support
3450// NaN values and real/abs comparison methods require specialized mapper.
3451
3452#define OP_INT_CUMMINMAX_FCN(F, OP) \
3453 template <typename T> \
3454 void F (const T *v, T *r, octave_idx_type n, const bool realabs) \
3455 { \
3456 if (! n) \
3457 return; \
3458 T tmp = v[0]; \
3459 octave_idx_type i = 1; \
3460 octave_idx_type j = 0; \
3461 if (realabs) \
3462 { \
3463 for (; i < n; i++) \
3464 if (v[i] OP tmp) \
3465 { \
3466 for (; j < i; j++) \
3467 r[j] = tmp; \
3468 tmp = v[i]; \
3469 } \
3470 } \
3471 else \
3472 { \
3473 for (; i < n; i++) \
3474 { \
3475 if (mappers_abs (v[i]) OP mappers_abs (tmp)) \
3476 { \
3477 for (; j < i; j++) \
3478 r[j] = tmp; \
3479 tmp = v[i]; \
3480 } \
3481 else if (mappers_abs (v[i]) == mappers_abs (tmp) && \
3482 v[i] OP tmp) \
3483 { \
3484 for (; j < i; j++) \
3485 r[j] = tmp; \
3486 tmp = v[i]; \
3487 } \
3488 } \
3489 } \
3490 for (; j < i; j++) \
3491 r[j] = tmp; \
3492 } \
3493 template <typename T> \
3494 void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n, \
3495 const bool realabs) \
3496 { \
3497 if (! n) \
3498 return; \
3499 T tmp = v[0]; \
3500 octave_idx_type tmpi = 0; \
3501 octave_idx_type i = 1; \
3502 octave_idx_type j = 0; \
3503 if (realabs) \
3504 { \
3505 for (; i < n; i++) \
3506 if (v[i] OP tmp) \
3507 { \
3508 for (; j < i; j++) \
3509 { \
3510 r[j] = tmp; \
3511 ri[j] = tmpi; \
3512 } \
3513 tmp = v[i]; \
3514 tmpi = i; \
3515 } \
3516 } \
3517 else \
3518 { \
3519 for (; i < n; i++) \
3520 { \
3521 if (mappers_abs (v[i]) OP mappers_abs (tmp)) \
3522 { \
3523 for (; j < i; j++) \
3524 { \
3525 r[j] = tmp; \
3526 ri[j] = tmpi; \
3527 } \
3528 tmp = v[i]; \
3529 tmpi = i; \
3530 } \
3531 else if (mappers_abs (v[i]) == mappers_abs (tmp) && \
3532 v[i] OP tmp) \
3533 { \
3534 for (; j < i; j++) \
3535 { \
3536 r[j] = tmp; \
3537 ri[j] = tmpi; \
3538 } \
3539 tmp = v[i]; \
3540 tmpi = i; \
3541 } \
3542 } \
3543 } \
3544 for (; j < i; j++) \
3545 { \
3546 r[j] = tmp; \
3547 ri[j] = tmpi; \
3548 } \
3549 }
3550
3553
3554#define OP_INT_CUMMINMAX_FCN2(F, OP) \
3555 template <typename T> \
3556 inline void \
3557 F (const T *v, T *r, octave_idx_type m, octave_idx_type n, \
3558 const bool realabs) \
3559 { \
3560 if (! n) \
3561 return; \
3562 const T *r0; \
3563 octave_idx_type j = 0; \
3564 for (octave_idx_type i = 0; i < m; i++) \
3565 { \
3566 r[i] = v[i]; \
3567 } \
3568 j++; \
3569 v += m; \
3570 r0 = r; \
3571 r += m; \
3572 if (realabs) \
3573 { \
3574 while (j < n) \
3575 { \
3576 for (octave_idx_type i = 0; i < m; i++) \
3577 { \
3578 if (v[i] OP r0[i]) \
3579 r[i] = v[i]; \
3580 else \
3581 r[i] = r0[i]; \
3582 } \
3583 j++; \
3584 v += m; \
3585 r0 = r; \
3586 r += m; \
3587 } \
3588 } \
3589 else \
3590 { \
3591 while (j < n) \
3592 { \
3593 for (octave_idx_type i = 0; i < m; i++) \
3594 { \
3595 if (mappers_abs (v[i]) OP mappers_abs (r0[i])) \
3596 r[i] = v[i]; \
3597 else if (mappers_abs (v[i]) == mappers_abs (r0[i]) && \
3598 v[i] OP r0[i]) \
3599 r[i] = v[i]; \
3600 else \
3601 r[i] = r0[i]; \
3602 } \
3603 j++; \
3604 v += m; \
3605 r0 = r; \
3606 r += m; \
3607 } \
3608 } \
3609 } \
3610 template <typename T> \
3611 inline void \
3612 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type m, \
3613 octave_idx_type n, const bool realabs) \
3614 { \
3615 if (! n) \
3616 return; \
3617 const T *r0; \
3618 const octave_idx_type *r0i; \
3619 octave_idx_type j = 0; \
3620 for (octave_idx_type i = 0; i < m; i++) \
3621 { \
3622 r[i] = v[i]; ri[i] = 0; \
3623 } \
3624 j++; \
3625 v += m; \
3626 r0 = r; \
3627 r += m; \
3628 r0i = ri; \
3629 ri += m; \
3630 if (realabs) \
3631 { \
3632 while (j < n) \
3633 { \
3634 for (octave_idx_type i = 0; i < m; i++) \
3635 { \
3636 if (v[i] OP r0[i]) \
3637 { \
3638 r[i] = v[i]; \
3639 ri[i] = j; \
3640 } \
3641 else \
3642 { \
3643 r[i] = r0[i]; \
3644 ri[i] = r0i[i]; \
3645 } \
3646 } \
3647 j++; \
3648 v += m; \
3649 r0 = r; \
3650 r += m; \
3651 r0i = ri; \
3652 ri += m; \
3653 } \
3654 } \
3655 else \
3656 { \
3657 while (j < n) \
3658 { \
3659 for (octave_idx_type i = 0; i < m; i++) \
3660 { \
3661 if (mappers_abs (v[i]) OP mappers_abs (r0[i])) \
3662 { \
3663 r[i] = v[i]; \
3664 ri[i] = j; \
3665 } \
3666 else if (mappers_abs (v[i]) == mappers_abs (r0[i]) \
3667 && v[i] OP r0[i]) \
3668 { \
3669 r[i] = v[i]; \
3670 ri[i] = j; \
3671 } \
3672 else \
3673 { \
3674 r[i] = r0[i]; \
3675 ri[i] = r0i[i]; \
3676 } \
3677 } \
3678 j++; \
3679 v += m; \
3680 r0 = r; \
3681 r += m; \
3682 r0i = ri; \
3683 ri += m; \
3684 } \
3685 } \
3686 }
3687
3690
3691#define OP_INT_CUMMINMAX_FCNN(F) \
3692 template <typename T> \
3693 inline void \
3694 F (const T *v, T *r, octave_idx_type l, octave_idx_type n, \
3695 octave_idx_type u, const bool realabs) \
3696 { \
3697 if (! n) \
3698 return; \
3699 if (l == 1) \
3700 { \
3701 for (octave_idx_type i = 0; i < u; i++) \
3702 { \
3703 F (v, r, n, realabs); \
3704 v += n; \
3705 r += n; \
3706 } \
3707 } \
3708 else \
3709 { \
3710 for (octave_idx_type i = 0; i < u; i++) \
3711 { \
3712 F (v, r, l, n, realabs); \
3713 v += l*n; \
3714 r += l*n; \
3715 } \
3716 } \
3717 } \
3718 template <typename T> \
3719 inline void \
3720 F (const T *v, T *r, octave_idx_type *ri, octave_idx_type l, \
3721 octave_idx_type n, octave_idx_type u, const bool realabs) \
3722 { \
3723 if (! n) \
3724 return; \
3725 if (l == 1) \
3726 { \
3727 for (octave_idx_type i = 0; i < u; i++) \
3728 { \
3729 F (v, r, ri, n, realabs); \
3730 v += n; \
3731 r += n; \
3732 ri += n; \
3733 } \
3734 } \
3735 else \
3736 { \
3737 for (octave_idx_type i = 0; i < u; i++) \
3738 { \
3739 F (v, r, ri, l, n, realabs); \
3740 v += l*n; \
3741 r += l*n; \
3742 ri += l*n; \
3743 } \
3744 } \
3745 }
3746
3749
3750template <typename T>
3751void mx_inline_diff (const T *v, T *r, octave_idx_type n,
3752 octave_idx_type order)
3753{
3754 switch (order)
3755 {
3756 case 1:
3757 for (octave_idx_type i = 0; i < n-1; i++)
3758 r[i] = v[i+1] - v[i];
3759 break;
3760 case 2:
3761 if (n > 1)
3762 {
3763 T lst = v[1] - v[0];
3764 for (octave_idx_type i = 0; i < n-2; i++)
3765 {
3766 T dif = v[i+2] - v[i+1];
3767 r[i] = dif - lst;
3768 lst = dif;
3769 }
3770 }
3771 break;
3772 default:
3773 {
3774 OCTAVE_LOCAL_BUFFER (T, buf, n-1);
3775
3776 for (octave_idx_type i = 0; i < n-1; i++)
3777 buf[i] = v[i+1] - v[i];
3778
3779 for (octave_idx_type o = 2; o <= order; o++)
3780 {
3781 for (octave_idx_type i = 0; i < n-o; i++)
3782 buf[i] = buf[i+1] - buf[i];
3783 }
3784
3785 for (octave_idx_type i = 0; i < n-order; i++)
3786 r[i] = buf[i];
3787 }
3788 }
3789}
3790
3791template <typename T>
3792void
3793mx_inline_diff (const T *v, T *r,
3795 octave_idx_type order)
3796{
3797 switch (order)
3798 {
3799 case 1:
3800 for (octave_idx_type i = 0; i < m*(n-1); i++)
3801 r[i] = v[i+m] - v[i];
3802 break;
3803 case 2:
3804 for (octave_idx_type i = 0; i < n-2; i++)
3805 {
3806 for (octave_idx_type j = i*m; j < i*m+m; j++)
3807 r[j] = (v[j+m+m] - v[j+m]) - (v[j+m] - v[j]);
3808 }
3809 break;
3810 default:
3811 {
3812 OCTAVE_LOCAL_BUFFER (T, buf, n-1);
3813
3814 for (octave_idx_type j = 0; j < m; j++)
3815 {
3816 for (octave_idx_type i = 0; i < n-1; i++)
3817 buf[i] = v[i*m+j+m] - v[i*m+j];
3818
3819 for (octave_idx_type o = 2; o <= order; o++)
3820 {
3821 for (octave_idx_type i = 0; i < n-o; i++)
3822 buf[i] = buf[i+1] - buf[i];
3823 }
3824
3825 for (octave_idx_type i = 0; i < n-order; i++)
3826 r[i*m+j] = buf[i];
3827 }
3828 }
3829 }
3830}
3831
3832template <typename T>
3833inline void
3834mx_inline_diff (const T *v, T *r,
3836 octave_idx_type order)
3837{
3838 if (! n) return;
3839 if (l == 1)
3840 {
3841 for (octave_idx_type i = 0; i < u; i++)
3842 {
3843 mx_inline_diff (v, r, n, order);
3844 v += n; r += n-order;
3845 }
3846 }
3847 else
3848 {
3849 for (octave_idx_type i = 0; i < u; i++)
3850 {
3851 mx_inline_diff (v, r, l, n, order);
3852 v += l*n;
3853 r += l*(n-order);
3854 }
3855 }
3856}
3857
3858// Assistant function
3859
3860inline void
3861get_extent_triplet (const dim_vector& dims, int& dim,
3863 octave_idx_type& u)
3864{
3865 octave_idx_type ndims = dims.ndims ();
3866 if (dim >= ndims)
3867 {
3868 l = dims.numel ();
3869 n = 1;
3870 u = 1;
3871 }
3872 else
3873 {
3874 if (dim < 0)
3875 dim = dims.first_non_singleton ();
3876
3877 // calculate extent triplet.
3878 l = 1, n = dims(dim), u = 1;
3879 for (octave_idx_type i = 0; i < dim; i++)
3880 l *= dims(i);
3881 for (octave_idx_type i = dim + 1; i < ndims; i++)
3882 u *= dims(i);
3883 }
3884}
3885
3886// For char and octave_int arrays
3887
3888template <typename R, typename T>
3889inline Array<R>
3890do_mx_red_op (const Array<T>& src, int dim,
3891 void (*mx_red_op) (const T *, R *, octave_idx_type,
3893{
3894 octave_idx_type l, n, u;
3895 dim_vector dims = src.dims ();
3896
3897 // Handle 0x0 empty array
3898 if (dims.ndims () == 2 && dims(0) == 0 && dims(1) == 0 && dim == -1)
3899 dims(1) = 1;
3900
3901 get_extent_triplet (dims, dim, l, n, u);
3902
3903 // Reduction operation reduces the array size.
3904 if (dim < dims.ndims ()) dims(dim) = 1;
3906
3907 Array<R> ret (dims);
3908 mx_red_op (src.data (), ret.rwdata (), l, n, u);
3909
3910 return ret;
3911}
3912
3913// For double, float, and complex arrays
3914
3915template <typename R, typename T>
3916inline Array<R>
3917do_mx_red_op (const Array<T>& src, int dim, bool nanflag,
3918 void (*mx_red_op) (const T *, R *, octave_idx_type,
3920{
3921 octave_idx_type l, n, u;
3922 dim_vector dims = src.dims ();
3923
3924 // Handle 0x0 empty array
3925 if (dims.ndims () == 2 && dims(0) == 0 && dims(1) == 0 && dim == -1)
3926 dims(1) = 1;
3927
3928 get_extent_triplet (dims, dim, l, n, u);
3929
3930 // Reduction operation reduces the array size.
3931 if (dim < dims.ndims ()) dims(dim) = 1;
3933
3934 Array<R> ret (dims);
3935 mx_red_op (src.data (), ret.rwdata (), l, n, u, nanflag);
3936
3937 return ret;
3938}
3939
3940template <typename R, typename T>
3941inline Array<R>
3942do_mx_flip_op (const Array<T>& src, int dim,
3943 void (*mx_flip_op) (const T *, R *, octave_idx_type,
3945{
3946 octave_idx_type l, n, u;
3947 const dim_vector& dims = src.dims ();
3948 get_extent_triplet (dims, dim, l, n, u);
3949
3950 // Cumulative operation doesn't reduce the array size.
3951 Array<R> ret (dims);
3952 mx_flip_op (src.data (), ret.rwdata (), l, n, u);
3953
3954 return ret;
3955}
3956
3957template <typename R, typename T>
3958inline Array<R>
3959do_mx_cum_op (const Array<T>& src, int dim,
3960 void (*mx_cum_op) (const T *, R *, octave_idx_type,
3962{
3963 octave_idx_type l, n, u;
3964 const dim_vector& dims = src.dims ();
3965 get_extent_triplet (dims, dim, l, n, u);
3966
3967 // Cumulative operation doesn't reduce the array size.
3968 Array<R> ret (dims);
3969 mx_cum_op (src.data (), ret.rwdata (), l, n, u);
3970
3971 return ret;
3972}
3973
3974template <typename R, typename T>
3975inline Array<R>
3976do_mx_cum_op (const Array<T>& src, int dim, bool nanflag,
3977 void (*mx_cum_op) (const T *, R *, octave_idx_type,
3979{
3980 octave_idx_type l, n, u;
3981 const dim_vector& dims = src.dims ();
3982 get_extent_triplet (dims, dim, l, n, u);
3983
3984 // Cumulative operation doesn't reduce the array size.
3985 Array<R> ret (dims);
3986 mx_cum_op (src.data (), ret.rwdata (), l, n, u, nanflag);
3987
3988 return ret;
3989}
3990
3991// For double, float, and complex arrays
3992
3993template <typename R>
3994inline Array<R>
3995do_mx_minmax_op (const Array<R>& src, int dim, bool nanflag, bool realabs,
3996 void (*mx_minmax_op) (const R *, R *, octave_idx_type,
3998 bool, bool))
3999{
4000 octave_idx_type l, n, u;
4001 dim_vector dims = src.dims ();
4002 get_extent_triplet (dims, dim, l, n, u);
4003
4004 // If the dimension is zero, we don't do anything.
4005 if (dim < dims.ndims () && dims(dim) != 0) dims(dim) = 1;
4007
4008 Array<R> ret (dims);
4009 mx_minmax_op (src.data (), ret.rwdata (), l, n, u, nanflag, realabs);
4010
4011 return ret;
4012}
4013
4014template <typename R>
4015inline Array<R>
4017 bool nanflag, bool realabs,
4018 void (*mx_minmax_op) (const R *, R *, octave_idx_type *,
4020 octave_idx_type, bool, bool))
4021{
4022 octave_idx_type l, n, u;
4023 dim_vector dims = src.dims ();
4024 get_extent_triplet (dims, dim, l, n, u);
4025
4026 // If the dimension is zero, we don't do anything.
4027 if (dim < dims.ndims () && dims(dim) != 0) dims(dim) = 1;
4029
4030 Array<R> ret (dims);
4031 if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
4032
4033 mx_minmax_op (src.data (), ret.rwdata (), idx.rwdata (),
4034 l, n, u, nanflag, realabs);
4035
4036 return ret;
4037}
4038
4039// For octave_int arrays
4040
4041template <typename R>
4042inline Array<R>
4043do_mx_minmax_op (const Array<R>& src, int dim, bool realabs,
4044 void (*mx_minmax_op) (const R *, R *, octave_idx_type,
4046{
4047 octave_idx_type l, n, u;
4048 dim_vector dims = src.dims ();
4049 get_extent_triplet (dims, dim, l, n, u);
4050
4051 // If the dimension is zero, we don't do anything.
4052 if (dim < dims.ndims () && dims(dim) != 0) dims(dim) = 1;
4054
4055 Array<R> ret (dims);
4056 mx_minmax_op (src.data (), ret.rwdata (), l, n, u, realabs);
4057
4058 return ret;
4059}
4060
4061template <typename R>
4062inline Array<R>
4064 bool realabs,
4065 void (*mx_minmax_op) (const R *, R *, octave_idx_type *,
4067 octave_idx_type, bool))
4068{
4069 octave_idx_type l, n, u;
4070 dim_vector dims = src.dims ();
4071 get_extent_triplet (dims, dim, l, n, u);
4072
4073 // If the dimension is zero, we don't do anything.
4074 if (dim < dims.ndims () && dims(dim) != 0) dims(dim) = 1;
4076
4077 Array<R> ret (dims);
4078 if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
4079
4080 mx_minmax_op (src.data (), ret.rwdata (), idx.rwdata (), l, n, u, realabs);
4081
4082 return ret;
4083}
4084
4085// For char arrays
4086
4087template <typename R>
4088inline Array<R>
4089do_mx_minmax_op (const Array<R>& src, int dim,
4090 void (*mx_minmax_op) (const R *, R *, octave_idx_type,
4092{
4093 octave_idx_type l, n, u;
4094 dim_vector dims = src.dims ();
4095 get_extent_triplet (dims, dim, l, n, u);
4096
4097 // If the dimension is zero, we don't do anything.
4098 if (dim < dims.ndims () && dims(dim) != 0) dims(dim) = 1;
4100
4101 Array<R> ret (dims);
4102 mx_minmax_op (src.data (), ret.rwdata (), l, n, u);
4103
4104 return ret;
4105}
4106
4107template <typename R>
4108inline Array<R>
4110 void (*mx_minmax_op) (const R *, R *, octave_idx_type *,
4113{
4114 octave_idx_type l, n, u;
4115 dim_vector dims = src.dims ();
4116 get_extent_triplet (dims, dim, l, n, u);
4117
4118 // If the dimension is zero, we don't do anything.
4119 if (dim < dims.ndims () && dims(dim) != 0) dims(dim) = 1;
4121
4122 Array<R> ret (dims);
4123 if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
4124
4125 mx_minmax_op (src.data (), ret.rwdata (), idx.rwdata (), l, n, u);
4126
4127 return ret;
4128}
4129
4130// For double, float, and complex arrays
4131
4132template <typename R>
4133inline Array<R>
4134do_mx_cumminmax_op (const Array<R>& src, int dim, bool nanflag, bool realabs,
4135 void (*mx_cumminmax_op) (const R *, R *, octave_idx_type,
4137 bool, bool))
4138{
4139 octave_idx_type l, n, u;
4140 const dim_vector& dims = src.dims ();
4141 get_extent_triplet (dims, dim, l, n, u);
4142
4143 Array<R> ret (dims);
4144 mx_cumminmax_op (src.data (), ret.rwdata (), l, n, u, nanflag, realabs);
4145
4146 return ret;
4147}
4148
4149template <typename R>
4150inline Array<R>
4152 bool nanflag, bool realabs,
4153 void (*mx_cumminmax_op) (const R *, R *, octave_idx_type *,
4155 octave_idx_type, bool, bool))
4156{
4157 octave_idx_type l, n, u;
4158 const dim_vector& dims = src.dims ();
4159 get_extent_triplet (dims, dim, l, n, u);
4160
4161 Array<R> ret (dims);
4162 if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
4163
4164 mx_cumminmax_op (src.data (), ret.rwdata (), idx.rwdata (),
4165 l, n, u, nanflag, realabs);
4166
4167 return ret;
4168}
4169
4170// For octave_int arrays
4171
4172template <typename R>
4173inline Array<R>
4174do_mx_cumminmax_op (const Array<R>& src, int dim, bool realabs,
4175 void (*mx_cumminmax_op) (const R *, R *, octave_idx_type,
4177 bool))
4178{
4179 octave_idx_type l, n, u;
4180 const dim_vector& dims = src.dims ();
4181 get_extent_triplet (dims, dim, l, n, u);
4182
4183 Array<R> ret (dims);
4184 mx_cumminmax_op (src.data (), ret.rwdata (), l, n, u, realabs);
4185
4186 return ret;
4187}
4188
4189template <typename R>
4190inline Array<R>
4192 bool realabs,
4193 void (*mx_cumminmax_op) (const R *, R *, octave_idx_type *,
4195 octave_idx_type, bool))
4196{
4197 octave_idx_type l, n, u;
4198 const dim_vector& dims = src.dims ();
4199 get_extent_triplet (dims, dim, l, n, u);
4200
4201 Array<R> ret (dims);
4202 if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
4203
4204 mx_cumminmax_op (src.data (), ret.rwdata (), idx.rwdata (),
4205 l, n, u, realabs);
4206
4207 return ret;
4208}
4209
4210template <typename R>
4211inline Array<R>
4212do_mx_diff_op (const Array<R>& src, int dim, octave_idx_type order,
4213 void (*mx_diff_op) (const R *, R *,
4216{
4217 octave_idx_type l, n, u;
4218 if (order <= 0)
4219 return src;
4220
4221 dim_vector dims = src.dims ();
4222
4223 get_extent_triplet (dims, dim, l, n, u);
4224 if (dim >= dims.ndims ())
4225 dims.resize (dim+1, 1);
4226
4227 if (dims(dim) <= order)
4228 {
4229 dims(dim) = 0;
4230 return Array<R> (dims);
4231 }
4232 else
4233 {
4234 dims(dim) -= order;
4235 }
4236
4237 Array<R> ret (dims);
4238 mx_diff_op (src.data (), ret.rwdata (), l, n, u, order);
4239
4240 return ret;
4241}
4242
4243// Fast extra-precise summation.
4244// Reference: T. Ogita, S. M. Rump, S. Oishi,
4245// "Accurate Sum And Dot Product",
4246// SIAM J. Sci. Computing, Vol. 26, 2005
4247
4248// Safe implementation (Standard API).
4249// Handles overflows to infinity correctly to prevent NaN pollution in error term.
4250template <typename T>
4251inline void
4252twosum_accum (T& s, T& e, const T& x)
4253{
4254 const T s1 = s + x;
4255
4256 // If sum becomes non-finite (overflow to +/-Inf), do NOT update the
4257 // compensation term (avoids Inf-Inf -> NaN in e).
4258 if (! octave::math::isfinite (s1))
4259 {
4260 s = s1;
4261 e = T (0);
4262 return;
4263 }
4264
4265 const T t = s1 - s;
4266 const T e1 = (s - (s1 - t)) + (x - t);
4267
4268 s = s1;
4269 e += e1;
4270}
4271
4272// Overload with three arguments for vectors (Contiguous / Vector)
4273// Used for column sums or simple vector reductions.
4274template <typename T>
4275inline T
4276mx_inline_xsum (const T *v, octave_idx_type n, bool nanflag)
4277{
4278 T s = 0, e = 0;
4279 const T zero = 0;
4280 bool posinf = false;
4281 bool neginf = false;
4282 bool seen_nan = false;
4283
4284 for (octave_idx_type i = 0; i < n; ++i)
4285 {
4286 const T val = v[i];
4287
4288 if (octave::math::isnan (val))
4289 {
4290 if (nanflag) // omitnan
4291 continue;
4292 else // normal sum: NaN dominates everything
4293 {
4294 seen_nan = true;
4295 break;
4296 }
4297 }
4298
4299 if (! octave::math::isinf (val))
4300 twosum_accum (s, e, val);
4301 else if (val > zero)
4302 posinf = true;
4303 else
4304 neginf = true;
4305 }
4306
4307 if (octave::math::isinf (s))
4308 {
4309 if (s > zero)
4310 posinf = true;
4311 else
4312 neginf = true;
4313 }
4314
4315 if (seen_nan)
4316 return std::numeric_limits<T>::quiet_NaN ();
4317
4318 if (posinf && neginf)
4319 return std::numeric_limits<T>::quiet_NaN ();
4320 else if (posinf)
4321 return std::numeric_limits<T>::infinity ();
4322 else if (neginf)
4323 return -std::numeric_limits<T>::infinity ();
4324
4325 return s + e;
4326}
4327
4328// Overload with five arguments for matrices (strided)
4329template <typename T>
4330inline void
4331mx_inline_xsum (const T *v, T *r,
4332 octave_idx_type l, octave_idx_type n, bool nanflag)
4333{
4334 const T zero = 0;
4335
4336 for (octave_idx_type i = 0; i < l; ++i)
4337 {
4338 T s = 0, e = 0;
4339 bool posinf = false;
4340 bool neginf = false;
4341 bool seen_nan = false;
4342
4343 for (octave_idx_type j = 0; j < n; ++j)
4344 {
4345 const T val = v[i + j*l];
4346
4347 if (octave::math::isnan (val))
4348 {
4349 if (nanflag) // omitnan
4350 continue;
4351 else
4352 {
4353 seen_nan = true;
4354 break;
4355 }
4356 }
4357
4358 if (! octave::math::isinf (val))
4359 twosum_accum (s, e, val);
4360 else if (val > zero)
4361 posinf = true;
4362 else
4363 neginf = true;
4364 }
4365
4366 if (octave::math::isinf (s))
4367 {
4368 if (s > zero)
4369 posinf = true;
4370 else
4371 neginf = true;
4372 }
4373
4374 if (seen_nan || (posinf && neginf))
4375 r[i] = std::numeric_limits<T>::quiet_NaN ();
4376 else if (posinf)
4377 r[i] = std::numeric_limits<T>::infinity ();
4378 else if (neginf)
4379 r[i] = -std::numeric_limits<T>::infinity ();
4380 else
4381 r[i] = s + e;
4382 }
4383}
4384
4386
4387#endif
Array< R > do_bsxfun1_op(const Array< X > &x, const Array< Y > &y, const bool nanflag, void(*op_vv)(std::size_t, R *, const X *, const Y *, const bool), void(*op_sv)(std::size_t, R *, X, const Y *, const bool), void(*op_vs)(std::size_t, R *, const X *, Y, const bool))
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))
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))
Array< R > do_bsxfun2_op(const Array< X > &x, const Array< Y > &y, const bool nanflag, const bool realabs, void(*op_vv)(std::size_t, R *, const X *, const Y *, const bool, const bool), void(*op_sv)(std::size_t, R *, X, const Y *, const bool, const bool), void(*op_vs)(std::size_t, R *, const X *, Y, const bool, const bool))
bool is_valid_inplace_bsxfun(const dim_vector &rdv, const dim_vector &xdv)
Definition bsxfun.h:58
bool is_valid_bsxfun(const dim_vector &xdv, const dim_vector &ydv)
Definition bsxfun.h:39
N Dimensional Array with copy-on-write semantics.
Definition Array-base.h:130
const dim_vector & dims() const
Return a const-reference so that dims ()(i) works efficiently.
Definition Array-base.h:529
const T * data() const
Size of the specified dimension.
Definition Array-base.h:687
T * rwdata()
Size of the specified dimension.
octave_idx_type numel() const
Number of elements in the array.
Definition Array-base.h:440
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:92
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
Definition dim-vector.h:341
void chop_trailing_singletons()
Definition dim-vector.h:162
void resize(int n, int fill_value=0)
Definition dim-vector.h:278
octave_idx_type ndims() const
Number of dimensions.
Definition dim-vector.h:263
int first_non_singleton(int def=0) const
Definition dim-vector.h:481
double double_value() const
ColumnVector real(const ComplexColumnVector &a)
ColumnVector imag(const ComplexColumnVector &a)
void mx_inline_map(std::size_t n, R *r, const X *x)
void mx_inline_intcummax(const T *v, T *r, octave_idx_type n, const bool realabs)
Array< R > do_mx_unary_op(const Array< X > &x, void(*op)(std::size_t, R *, const X *))
void mx_inline_intmin(const T *v, T *r, octave_idx_type n, const bool realabs)
bool mx_inline_any(const T *v, octave_idx_type n)
void mx_inline_div2(std::size_t n, R *r, const X *x)
bool do_mx_check(const Array< T > &a, bool(*op)(std::size_t, const T *))
void mx_inline_sub2(std::size_t n, R *r, const X *x)
subst_template_param< std::complex, T, double >::type mx_inline_dsumsq(const T *v, octave_idx_type n)
void mx_inline_notzero(std::size_t n, bool *r, const X *x)
Definition mx-inlines.cc:82
#define DEFMXCMPOP(F, OP)
void mx_inline_xmin(std::size_t n, T *r, const T *x, const T *y)
void mx_inline_not_and(std::size_t n, bool *r, const X *x, const Y *y)
#define OP_RED_ANYR(ac, el)
bool xis_false(T x)
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)
#define OP_CUM_FCN2(F, TSRC, TRES, OP)
void mx_inline_le(std::size_t n, bool *r, const X *x, const Y *y)
#define OP_ROW_SHORT_CIRCUIT(F, PRED, ZERO)
#define OP_RED_SUMSQC(ac, el)
void mx_inline_cummin(const T *v, T *r, octave_idx_type n, const bool nanflag, const bool realabs)
void mx_inline_intcummin(const T *v, T *r, octave_idx_type n, const bool realabs)
#define OP_INT_CUMMINMAX_FCNN(F)
void mx_inline_add2(std::size_t n, R *r, const X *x)
Array< R > do_ms_binary_op(const Array< X > &x, const Y &y, void(*op)(std::size_t, R *, const X *, Y))
void mx_inline_and_not(std::size_t n, bool *r, const X *x, const Y *y)
void mx_inline_ccummax(const std::complex< T > *v, std::complex< T > *r, octave_idx_type n, const bool nanflag, const bool realabs)
#define OP_RED_FCNN(F, TSRC, TRES)
#define OP_CUMMINMAX_FCN(F, OP)
void mx_inline_sub(std::size_t n, R *r, const X *x, const Y *y)
#define OP_RED_NAN_FCN2(F, TSRC, TRES, OP, ZERO)
void mx_inline_uminus(std::size_t n, R *r, const X *x)
Definition mx-inlines.cc:57
Array< R > & do_mx_inplace_op(Array< R > &r, void(*op)(std::size_t, R *))
void mx_inline_cummax(const T *v, T *r, octave_idx_type n, const bool nanflag, const bool realabs)
#define OP_CMINMAX_FCNN(F)
subst_template_param< std::complex, T, double >::type mx_inline_dprod(const T *v, octave_idx_type n)
void mx_inline_uminus2(std::size_t n, R *r)
Definition mx-inlines.cc:65
void mx_inline_eq(std::size_t n, bool *r, const X *x, const Y *y)
#define OP_RED_ALLC(ac, el)
#define OP_CHMINMAX_FCN2(F, OP)
#define OP_INT_CUMMINMAX_FCN(F, OP)
#define DEFMXBOOLOP(F, NOT1, OP, NOT2)
void mx_inline_ge(std::size_t n, bool *r, const X *x, const Y *y)
bool mx_inline_all_real(std::size_t n, const std::complex< T > *x)
#define PROMOTE_DOUBLE(T)
void mx_inline_min(const T *v, T *r, octave_idx_type n, const bool nanflag, const bool realabs)
void mx_inline_or(std::size_t n, bool *r, const X *x, const Y *y)
void mx_inline_chmin(const T *v, T *r, octave_idx_type n)
void mx_inline_div(std::size_t n, R *r, const X *x, const Y *y)
#define OP_MINMAX_FCN2(F, OP)
Array< R > do_mx_minmax_op(const Array< R > &src, int dim, bool nanflag, bool realabs, void(*mx_minmax_op)(const R *, R *, octave_idx_type, octave_idx_type, octave_idx_type, bool, bool))
void mx_inline_cumprod(const T *v, T *r, octave_idx_type n)
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))
#define OP_INT_CUMMINMAX_FCN2(F, OP)
void mx_inline_cumsum(const T *v, T *r, octave_idx_type n)
void twosum_accum(T &s, T &e, const T &x)
bool mx_inline_any_nan(std::size_t n, const T *x)
void mx_inline_and2(std::size_t n, bool *r, const X *x)
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))
#define DEFMXBINOPEQ(F, OP)
#define OP_CUM_FCN(F, TSRC, TRES, OP)
void mx_inline_iszero(std::size_t n, bool *r, const X *x)
Definition mx-inlines.cc:73
void mx_inline_xmax(std::size_t n, T *r, const T *x, const T *y)
void mx_inline_any_r(const T *v, bool *r, octave_idx_type m, octave_idx_type n)
bool xis_true(T x)
T mx_inline_count(const bool *v, octave_idx_type n)
void mx_inline_and(std::size_t n, bool *r, const X *x, const Y *y)
#define OP_MINMAX_FCNN(F)
#define OP_CCUMMINMAX_FCNN(F)
void mx_inline_add(std::size_t n, R *r, const X *x, const Y *y)
void mx_inline_intmax(const T *v, T *r, octave_idx_type n, const bool realabs)
#define OP_CUM_NAN_FCN(F, TSRC, TRES, OP, ZERO)
void mx_inline_not(std::size_t n, bool *r, const X *x)
#define OP_CUM_FCNN(F, TSRC, TRES)
Array< R > & do_ms_inplace_op(Array< R > &r, const X &x, void(*op)(std::size_t, R *, X))
void mx_inline_cmax(const std::complex< T > *v, std::complex< T > *r, octave_idx_type n, const bool nanflag, const bool realabs)
void mx_inline_not_or(std::size_t n, bool *r, const X *x, const Y *y)
void mx_inline_flip(const T *v, T *r, octave_idx_type n)
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)
Array< R > do_mx_cumminmax_op(const Array< R > &src, int dim, bool nanflag, bool realabs, void(*mx_cumminmax_op)(const R *, R *, octave_idx_type, octave_idx_type, octave_idx_type, bool, bool))
#define OP_CUMMINMAX_FCN2(F, OP)
Array< R > do_sm_binary_op(const X &x, const Array< Y > &y, void(*op)(std::size_t, R *, X, const Y *))
void op_dble_prod(double &ac, float el)
#define OP_CMINMAX_FCN(F, OP)
T cabsq(const std::complex< T > &c)
#define OP_RED_FCN2(F, TSRC, TRES, OP, ZERO)
#define OP_INT_MINMAX_FCNN(F)
void mx_inline_gt(std::size_t n, bool *r, const X *x, const Y *y)
void op_dble_sum(double &ac, float el)
#define OP_INT_MINMAX_FCN(F, OP)
void mx_inline_all_r(const T *v, bool *r, octave_idx_type m, octave_idx_type n)
#define OP_CCUMMINMAX_FCN(F, OP)
void mx_inline_cumcount(const bool *v, T *r, octave_idx_type n)
void mx_inline_mul2(std::size_t n, R *r, const X *x)
#define OP_RED_ANYC(ac, el)
void mx_inline_lt(std::size_t n, bool *r, const X *x, const Y *y)
void mx_inline_cmin(const std::complex< T > *v, std::complex< T > *r, octave_idx_type n, const bool nanflag, const bool realabs)
#define OP_CUM_NAN_FCN2(F, TSRC, TRES, OP, ZERO)
void mx_inline_max(const T *v, T *r, octave_idx_type n, const bool nanflag, const bool realabs)
bool mx_inline_all(const T *v, octave_idx_type n)
bool mx_inline_equal(std::size_t n, const T1 *x, const T2 *y)
void mx_inline_real(std::size_t n, T *r, const std::complex< T > *x)
#define OP_CHMINMAX_FCN(F, OP)
T mx_inline_sumsq(const T *v, octave_idx_type n)
Array< R > do_mx_unary_map(const Array< X > &x)
#define OP_CUMMINMAX_FCNN(F)
bool mx_inline_any_positive(std::size_t n, const T *x)
void mx_inline_imag(std::size_t n, T *r, const std::complex< T > *x)
T mx_inline_sum(const T *v, octave_idx_type n)
#define OP_CUM_NAN_FCNN(F, TSRC, TRES)
void mx_inline_ccummin(const std::complex< T > *v, std::complex< T > *r, octave_idx_type n, const bool nanflag, const bool realabs)
#define OP_RED_NAN_FCN(F, TSRC, TRES, OP, ZERO)
#define OP_CMINMAX_FCN2(F, OP)
#define OP_CHMINMAX_FCNN(F)
#define OP_RED_SUM(ac, el)
T mx_inline_prod(const T *v, octave_idx_type n)
#define OP_MINMAX_FCN(F, OP)
#define OP_RED_FCN(F, TSRC, TRES, OP, ZERO)
void mx_inline_chmax(const T *v, T *r, octave_idx_type n)
void mx_inline_diff(const T *v, T *r, octave_idx_type n, octave_idx_type order)
void mx_inline_or2(std::size_t n, bool *r, const X *x)
bool mx_inline_all_finite(std::size_t n, const T *x)
void mx_inline_mul(std::size_t n, R *r, const X *x, const Y *y)
#define OP_RED_PROD(ac, el)
bool mx_inline_any_negative(std::size_t n, const T *x)
bool logical_value(T x)
#define OP_CCUMMINMAX_FCN2(F, OP)
Array< R > do_mx_flip_op(const Array< T > &src, int dim, void(*mx_flip_op)(const T *, R *, octave_idx_type, octave_idx_type, octave_idx_type))
T mx_inline_xsum(const T *v, octave_idx_type n, bool nanflag)
void get_extent_triplet(const dim_vector &dims, int &dim, octave_idx_type &l, octave_idx_type &n, octave_idx_type &u)
#define OP_INT_MINMAX_FCN2(F, OP)
void mx_inline_pow(std::size_t n, R *r, const X *x, const Y *y)
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))
void mx_inline_fill(std::size_t n, R *r, S s)
Definition mx-inlines.cc:49
void mx_inline_not2(std::size_t n, bool *r)
void mx_inline_ne(std::size_t n, bool *r, const X *x, const Y *y)
#define OP_RED_ALLR(ac, el)
#define OP_RED_NAN_FCNN(F, TSRC, TRES)
subst_template_param< std::complex, T, double >::type mx_inline_dsum(const T *v, octave_idx_type n)
#define OP_RED_SUMSQ(ac, el)
void mx_inline_or_not(std::size_t n, bool *r, const X *x, const Y *y)
#define DEFMXBINOP(F, OP)
Definition mx-inlines.cc:89
std::complex< double > Complex
Definition oct-cmplx.h:33
std::complex< float > FloatComplex
Definition oct-cmplx.h:34
octave_int< T > pow(const octave_int< T > &a, const octave_int< T > &b)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition oct-locbuf.h:44
F77_RET_T const F77_DBLE * x