GNU Octave 11.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
bsxfun-defs.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 2009-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_bsxfun_defs_h)
27#define octave_bsxfun_defs_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 <algorithm>
33
34#include "dim-vector.h"
35#include "mx-inlines.cc"
36#include "oct-error.h"
37#include "oct-locbuf.h"
38
39template <typename R, typename X, typename Y>
40Array<R>
41do_bsxfun_op (const Array<X>& x, const Array<Y>& y,
42 void (*op_vv) (std::size_t, R *, const X *, const Y *),
43 void (*op_sv) (std::size_t, R *, X, const Y *),
44 void (*op_vs) (std::size_t, R *, const X *, Y))
45{
46 int nd = std::max (x.ndims (), y.ndims ());
47 const dim_vector& dvx = x.dims ().redim (nd);
48 const dim_vector& dvy = y.dims ().redim (nd);
49
50 // Construct the result dimensions.
51 dim_vector dvr;
52 dvr.resize (nd);
53 for (int i = 0; i < nd; i++)
54 {
55 octave_idx_type xk = dvx(i);
56 octave_idx_type yk = dvy(i);
57 if (xk == 1)
58 dvr(i) = yk;
59 else if (yk == 1 || xk == yk)
60 dvr(i) = xk;
61 else
62 (*current_liboctave_error_handler)
63 ("bsxfun: nonconformant dimensions: %s and %s",
64 x.dims ().str ().c_str (), y.dims ().str ().c_str ());
65 }
66
67 Array<R> retval (dvr);
68
69 const X *xvec = x.data ();
70 const Y *yvec = y.data ();
71 R *rvec = retval.rwdata ();
72
73 // Fold the common leading dimensions.
74 octave_idx_type start, ldr = 1;
75 for (start = 0; start < nd; start++)
76 {
77 if (dvx(start) != dvy(start))
78 break;
79 ldr *= dvr(start);
80 }
81
82 if (retval.isempty ())
83 ; // do nothing
84 else if (start == nd)
85 op_vv (retval.numel (), rvec, xvec, yvec);
86 else
87 {
88 // Determine the type of the low-level loop.
89 bool xsing = false;
90 bool ysing = false;
91 if (ldr == 1)
92 {
93 xsing = dvx(start) == 1;
94 ysing = dvy(start) == 1;
95 if (xsing || ysing)
96 {
97 ldr *= dvx(start) * dvy(start);
98 start++;
99 }
100 }
101 dim_vector cdvx = dvx.cumulative ();
102 dim_vector cdvy = dvy.cumulative ();
103 // Nullify singleton dims to achieve a spread effect.
104 for (int i = std::max (start, octave_idx_type (1)); i < nd; i++)
105 {
106 if (dvx(i) == 1)
107 cdvx(i-1) = 0;
108 if (dvy(i) == 1)
109 cdvy(i-1) = 0;
110 }
111
112 octave_idx_type niter = dvr.numel (start);
113 // The index array.
114 OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, idx, nd, 0);
115 for (octave_idx_type iter = 0; iter < niter; iter++)
116 {
117 octave_quit ();
118
119 // Compute indices.
120 // FIXME: performance impact noticeable?
121 octave_idx_type xidx = cdvx.cum_compute_index (idx);
122 octave_idx_type yidx = cdvy.cum_compute_index (idx);
123 octave_idx_type ridx = dvr.compute_index (idx);
124
125 // Apply the low-level loop.
126 if (xsing)
127 op_sv (ldr, rvec + ridx, xvec[xidx], yvec + yidx);
128 else if (ysing)
129 op_vs (ldr, rvec + ridx, xvec + xidx, yvec[yidx]);
130 else
131 op_vv (ldr, rvec + ridx, xvec + xidx, yvec + yidx);
132
133 dvr.increment_index (idx + start, start);
134 }
135 }
136
137 return retval;
138}
139
140template <typename R, typename X>
141void
142do_inplace_bsxfun_op (Array<R>& r, const Array<X>& x,
143 void (*op_vv) (std::size_t, R *, const X *),
144 void (*op_vs) (std::size_t, R *, X))
145{
146 const dim_vector& dvr = r.dims ();
147 dim_vector dvx = x.dims ();
148 octave_idx_type nd = r.ndims ();
149 dvx = dvx.redim (nd);
150
151 const X *xvec = x.data ();
152 R *rvec = r.rwdata ();
153
154 // Fold the common leading dimensions.
155 octave_idx_type start, ldr = 1;
156 for (start = 0; start < nd; start++)
157 {
158 if (dvr(start) != dvx(start))
159 break;
160 ldr *= dvr(start);
161 }
162
163 if (r.isempty ())
164 ; // do nothing
165 else if (start == nd)
166 op_vv (r.numel (), rvec, xvec);
167 else
168 {
169 // Determine the type of the low-level loop.
170 bool xsing = false;
171 if (ldr == 1)
172 {
173 xsing = dvx(start) == 1;
174 if (xsing)
175 {
176 ldr *= dvr(start) * dvx(start);
177 start++;
178 }
179 }
180
181 dim_vector cdvx = dvx.cumulative ();
182 // Nullify singleton dims to achieve a spread effect.
183 for (int i = std::max (start, octave_idx_type (1)); i < nd; i++)
184 {
185 if (dvx(i) == 1)
186 cdvx(i-1) = 0;
187 }
188
189 octave_idx_type niter = dvr.numel (start);
190 // The index array.
191 OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, idx, nd, 0);
192 for (octave_idx_type iter = 0; iter < niter; iter++)
193 {
194 octave_quit ();
195
196 // Compute indices.
197 // FIXME: performance impact noticeable?
198 octave_idx_type xidx = cdvx.cum_compute_index (idx);
199 octave_idx_type ridx = dvr.compute_index (idx);
200
201 // Apply the low-level loop.
202 if (xsing)
203 op_vs (ldr, rvec + ridx, xvec[xidx]);
204 else
205 op_vv (ldr, rvec + ridx, xvec + xidx);
206
207 dvr.increment_index (idx + start, start);
208 }
209 }
210}
211
212// Specialization for min/max with extra nanflag and realabs optional arguments
213
214template <typename R, typename X, typename Y>
215Array<R>
216do_bsxfun1_op (const Array<X>& x, const Array<Y>& y, const bool nanflag,
217 void (*op_vv) (std::size_t, R *, const X *, const Y *,
218 const bool),
219 void (*op_sv) (std::size_t, R *, X, const Y *, const bool),
220 void (*op_vs) (std::size_t, R *, const X *, Y, const bool))
221{
222 int nd = std::max (x.ndims (), y.ndims ());
223 const dim_vector& dvx = x.dims ().redim (nd);
224 const dim_vector& dvy = y.dims ().redim (nd);
225
226 // Construct the result dimensions.
227 dim_vector dvr;
228 dvr.resize (nd);
229 for (int i = 0; i < nd; i++)
230 {
231 octave_idx_type xk = dvx(i);
232 octave_idx_type yk = dvy(i);
233 if (xk == 1)
234 dvr(i) = yk;
235 else if (yk == 1 || xk == yk)
236 dvr(i) = xk;
237 else
238 (*current_liboctave_error_handler)
239 ("bsxfun: nonconformant dimensions: %s and %s",
240 x.dims ().str ().c_str (), y.dims ().str ().c_str ());
241 }
242
243 Array<R> retval (dvr);
244
245 const X *xvec = x.data ();
246 const Y *yvec = y.data ();
247 R *rvec = retval.rwdata ();
248
249 // Fold the common leading dimensions.
250 octave_idx_type start, ldr = 1;
251 for (start = 0; start < nd; start++)
252 {
253 if (dvx(start) != dvy(start))
254 break;
255 ldr *= dvr(start);
256 }
257
258 if (retval.isempty ())
259 ; // do nothing
260 else if (start == nd)
261 op_vv (retval.numel (), rvec, xvec, yvec, nanflag);
262 else
263 {
264 // Determine the type of the low-level loop.
265 bool xsing = false;
266 bool ysing = false;
267 if (ldr == 1)
268 {
269 xsing = dvx(start) == 1;
270 ysing = dvy(start) == 1;
271 if (xsing || ysing)
272 {
273 ldr *= dvx(start) * dvy(start);
274 start++;
275 }
276 }
277 dim_vector cdvx = dvx.cumulative ();
278 dim_vector cdvy = dvy.cumulative ();
279 // Nullify singleton dims to achieve a spread effect.
280 for (int i = std::max (start, octave_idx_type (1)); i < nd; i++)
281 {
282 if (dvx(i) == 1)
283 cdvx(i-1) = 0;
284 if (dvy(i) == 1)
285 cdvy(i-1) = 0;
286 }
287
288 octave_idx_type niter = dvr.numel (start);
289 // The index array.
290 OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, idx, nd, 0);
291 for (octave_idx_type iter = 0; iter < niter; iter++)
292 {
293 octave_quit ();
294
295 // Compute indices.
296 // FIXME: performance impact noticeable?
297 octave_idx_type xidx = cdvx.cum_compute_index (idx);
298 octave_idx_type yidx = cdvy.cum_compute_index (idx);
299 octave_idx_type ridx = dvr.compute_index (idx);
300
301 // Apply the low-level loop.
302 if (xsing)
303 op_sv (ldr, rvec + ridx, xvec[xidx], yvec + yidx, nanflag);
304 else if (ysing)
305 op_vs (ldr, rvec + ridx, xvec + xidx, yvec[yidx], nanflag);
306 else
307 op_vv (ldr, rvec + ridx, xvec + xidx, yvec + yidx, nanflag);
308
309 dvr.increment_index (idx + start, start);
310 }
311 }
312
313 return retval;
314}
315
316template <typename R, typename X, typename Y>
317Array<R>
318do_bsxfun2_op (const Array<X>& x, const Array<Y>& y, const bool nanflag,
319 const bool realabs,
320 void (*op_vv) (std::size_t, R *, const X *, const Y *,
321 const bool, const bool),
322 void (*op_sv) (std::size_t, R *, X, const Y *,
323 const bool, const bool),
324 void (*op_vs) (std::size_t, R *, const X *, Y,
325 const bool, const bool))
326{
327 int nd = std::max (x.ndims (), y.ndims ());
328 const dim_vector& dvx = x.dims ().redim (nd);
329 const dim_vector& dvy = y.dims ().redim (nd);
330
331 // Construct the result dimensions.
332 dim_vector dvr;
333 dvr.resize (nd);
334 for (int i = 0; i < nd; i++)
335 {
336 octave_idx_type xk = dvx(i);
337 octave_idx_type yk = dvy(i);
338 if (xk == 1)
339 dvr(i) = yk;
340 else if (yk == 1 || xk == yk)
341 dvr(i) = xk;
342 else
343 (*current_liboctave_error_handler)
344 ("bsxfun: nonconformant dimensions: %s and %s",
345 x.dims ().str ().c_str (), y.dims ().str ().c_str ());
346 }
347
348 Array<R> retval (dvr);
349
350 const X *xvec = x.data ();
351 const Y *yvec = y.data ();
352 R *rvec = retval.rwdata ();
353
354 // Fold the common leading dimensions.
355 octave_idx_type start, ldr = 1;
356 for (start = 0; start < nd; start++)
357 {
358 if (dvx(start) != dvy(start))
359 break;
360 ldr *= dvr(start);
361 }
362
363 if (retval.isempty ())
364 ; // do nothing
365 else if (start == nd)
366 op_vv (retval.numel (), rvec, xvec, yvec, nanflag, realabs);
367 else
368 {
369 // Determine the type of the low-level loop.
370 bool xsing = false;
371 bool ysing = false;
372 if (ldr == 1)
373 {
374 xsing = dvx(start) == 1;
375 ysing = dvy(start) == 1;
376 if (xsing || ysing)
377 {
378 ldr *= dvx(start) * dvy(start);
379 start++;
380 }
381 }
382 dim_vector cdvx = dvx.cumulative ();
383 dim_vector cdvy = dvy.cumulative ();
384 // Nullify singleton dims to achieve a spread effect.
385 for (int i = std::max (start, octave_idx_type (1)); i < nd; i++)
386 {
387 if (dvx(i) == 1)
388 cdvx(i-1) = 0;
389 if (dvy(i) == 1)
390 cdvy(i-1) = 0;
391 }
392
393 octave_idx_type niter = dvr.numel (start);
394 // The index array.
395 OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, idx, nd, 0);
396 for (octave_idx_type iter = 0; iter < niter; iter++)
397 {
398 octave_quit ();
399
400 // Compute indices.
401 // FIXME: performance impact noticeable?
402 octave_idx_type xidx = cdvx.cum_compute_index (idx);
403 octave_idx_type yidx = cdvy.cum_compute_index (idx);
404 octave_idx_type ridx = dvr.compute_index (idx);
405
406 // Apply the low-level loop.
407 if (xsing)
408 op_sv (ldr, rvec + ridx, xvec[xidx], yvec + yidx, nanflag, realabs);
409 else if (ysing)
410 op_vs (ldr, rvec + ridx, xvec + xidx, yvec[yidx], nanflag, realabs);
411 else
412 op_vv (ldr, rvec + ridx, xvec + xidx, yvec + yidx, nanflag, realabs);
413
414 dvr.increment_index (idx + start, start);
415 }
416 }
417
418 return retval;
419}
420
421#define BSXFUN_OP_DEF(OP, ARRAY) \
422 ARRAY bsxfun_ ## OP (const ARRAY& x, const ARRAY& y)
423
424#define BSXFUN1_OP_DEF(OP, ARRAY) \
425 ARRAY bsxfun_ ## OP (const ARRAY& x, const ARRAY& y, \
426 const bool nanflag)
427
428#define BSXFUN2_OP_DEF(OP, ARRAY) \
429 ARRAY bsxfun_ ## OP (const ARRAY& x, const ARRAY& y, \
430 const bool nanflag, const bool realabs)
431
432#define BSXFUN_OP2_DEF(OP, ARRAY, ARRAY1, ARRAY2) \
433 ARRAY bsxfun_ ## OP (const ARRAY1& x, const ARRAY2& y)
434
435#define BSXFUN_REL_DEF(OP, ARRAY) \
436 boolNDArray bsxfun_ ## OP (const ARRAY& x, const ARRAY& y)
437
438#define BSXFUN_OP_DEF_MXLOOP(OP, ARRAY, LOOP) \
439 BSXFUN_OP_DEF(OP, ARRAY) \
440 { return do_bsxfun_op<ARRAY::element_type, ARRAY::element_type, ARRAY::element_type> \
441 (x, y, LOOP, LOOP, LOOP); }
442
443#define BSXFUN1_OP_DEF_MXLOOP(OP, ARRAY, LOOP) \
444 BSXFUN1_OP_DEF(OP, ARRAY) \
445 { return do_bsxfun1_op<ARRAY::element_type, ARRAY::element_type, ARRAY::element_type> \
446 (x, y, nanflag, LOOP, LOOP, LOOP); }
447
448#define BSXFUN2_OP_DEF_MXLOOP(OP, ARRAY, LOOP) \
449 BSXFUN2_OP_DEF(OP, ARRAY) \
450 { return do_bsxfun2_op<ARRAY::element_type, ARRAY::element_type, ARRAY::element_type> \
451 (x, y, nanflag, realabs, LOOP, LOOP, LOOP); }
452
453#define BSXFUN_OP2_DEF_MXLOOP(OP, ARRAY, ARRAY1, ARRAY2, LOOP) \
454 BSXFUN_OP2_DEF(OP, ARRAY, ARRAY1, ARRAY2) \
455 { return do_bsxfun_op<ARRAY::element_type, ARRAY1::element_type, ARRAY2::element_type> \
456 (x, y, LOOP, LOOP, LOOP); }
457
458#define BSXFUN_REL_DEF_MXLOOP(OP, ARRAY, LOOP) \
459 BSXFUN_REL_DEF(OP, ARRAY) \
460 { return do_bsxfun_op<bool, ARRAY::element_type, ARRAY::element_type> \
461 (x, y, LOOP, LOOP, LOOP); }
462
463#define BSXFUN_STDOP_DEFS_MXLOOP(ARRAY) \
464 BSXFUN_OP_DEF_MXLOOP (add, ARRAY, mx_inline_add) \
465 BSXFUN_OP_DEF_MXLOOP (sub, ARRAY, mx_inline_sub) \
466 BSXFUN_OP_DEF_MXLOOP (mul, ARRAY, mx_inline_mul) \
467 BSXFUN_OP_DEF_MXLOOP (div, ARRAY, mx_inline_div) \
468 BSXFUN_OP_DEF_MXLOOP (min, ARRAY, mx_inline_xmin) \
469 BSXFUN_OP_DEF_MXLOOP (max, ARRAY, mx_inline_xmax) \
470 BSXFUN1_OP_DEF_MXLOOP (min, ARRAY, mx_inline_xmin) \
471 BSXFUN1_OP_DEF_MXLOOP (max, ARRAY, mx_inline_xmax) \
472 BSXFUN2_OP_DEF_MXLOOP (min, ARRAY, mx_inline_xmin) \
473 BSXFUN2_OP_DEF_MXLOOP (max, ARRAY, mx_inline_xmax)
474
475#define BSXFUN_STDREL_DEFS_MXLOOP(ARRAY) \
476 BSXFUN_REL_DEF_MXLOOP (eq, ARRAY, mx_inline_eq) \
477 BSXFUN_REL_DEF_MXLOOP (ne, ARRAY, mx_inline_ne) \
478 BSXFUN_REL_DEF_MXLOOP (lt, ARRAY, mx_inline_lt) \
479 BSXFUN_REL_DEF_MXLOOP (le, ARRAY, mx_inline_le) \
480 BSXFUN_REL_DEF_MXLOOP (gt, ARRAY, mx_inline_gt) \
481 BSXFUN_REL_DEF_MXLOOP (ge, ARRAY, mx_inline_ge)
482
483//For bsxfun power with mixed integer/float types
484#define BSXFUN_POW_MIXED_MXLOOP(INT_TYPE) \
485 BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, INT_TYPE, NDArray, mx_inline_pow) \
486 BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, INT_TYPE, FloatNDArray, mx_inline_pow) \
487 BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, NDArray, INT_TYPE, mx_inline_pow) \
488 BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, FloatNDArray, INT_TYPE, mx_inline_pow)
489
490#endif
do_bsxfun1_op
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))
Definition bsxfun-defs.cc:216
do_bsxfun_op
Array< R > do_bsxfun_op(const Array< X > &x, const Array< Y > &y, void(*op_vv)(std::size_t, R *, const X *, const Y *), void(*op_sv)(std::size_t, R *, X, const Y *), void(*op_vs)(std::size_t, R *, const X *, Y))
Definition bsxfun-defs.cc:41
do_inplace_bsxfun_op
void do_inplace_bsxfun_op(Array< R > &r, const Array< X > &x, void(*op_vv)(std::size_t, R *, const X *), void(*op_vs)(std::size_t, R *, X))
Definition bsxfun-defs.cc:142
do_bsxfun2_op
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))
Definition bsxfun-defs.cc:318
Array
N Dimensional Array with copy-on-write semantics.
Definition Array-base.h:130
Array::dims
const dim_vector & dims() const
Return a const-reference so that dims ()(i) works efficiently.
Definition Array-base.h:529
Array::ndims
int ndims() const
Size of the specified dimension.
Definition Array-base.h:701
Array::isempty
bool isempty() const
Size of the specified dimension.
Definition Array-base.h:674
Array::data
const T * data() const
Size of the specified dimension.
Definition Array-base.h:687
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1790
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array-base.h:440
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:92
dim_vector::compute_index
octave_idx_type compute_index(const octave_idx_type *idx) const
Linear index from an index tuple.
Definition dim-vector.h:493
dim_vector::str
std::string str(char sep='x') const
Definition dim-vector.cc:58
dim_vector::numel
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
Definition dim-vector.h:341
dim_vector::resize
void resize(int n, int fill_value=0)
Definition dim-vector.h:278
dim_vector::increment_index
int increment_index(octave_idx_type *idx, int start=0) const
Increment a multi-dimensional index tuple, optionally starting from an offset position and return the...
Definition dim-vector.h:510
dim_vector::cumulative
dim_vector cumulative() const
Return cumulative dimensions.
Definition dim-vector.h:525
dim_vector::cum_compute_index
octave_idx_type cum_compute_index(const octave_idx_type *idx) const
Compute a linear index from an index tuple.
Definition dim-vector.h:540
dim_vector::redim
dim_vector redim(int n) const
Force certain dimensionality, preserving numel ().
Definition dim-vector.cc:197
octave_idx_type
current_liboctave_error_handler
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition lo-error.c:41
OCTAVE_LOCAL_BUFFER_INIT
#define OCTAVE_LOCAL_BUFFER_INIT(T, buf, size, value)
Definition oct-locbuf.h:50
x
F77_RET_T const F77_DBLE * x
Definition slatec-proto.h:38