GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
bsxfun-defs.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2009-2018 Jaroslav Hajek
4 Copyright (C) 2009 VZLU Prague
5 
6 This file is part of Octave.
7 
8 Octave is free software: you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <https://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if ! defined (octave_bsxfun_defs_h)
25 #define octave_bsxfun_defs_h 1
26 
27 // This file should *not* include config.h. It is only included in other C++
28 // source files that should have included config.h before including this file.
29 
30 #include <algorithm>
31 
32 #include "dim-vector.h"
33 #include "lo-error.h"
34 #include "mx-inlines.cc"
35 #include "oct-locbuf.h"
36 
37 template <typename R, typename X, typename Y>
39 do_bsxfun_op (const Array<X>& x, const Array<Y>& y,
40  void (*op_vv) (size_t, R *, const X *, const Y *),
41  void (*op_sv) (size_t, R *, X, const Y *),
42  void (*op_vs) (size_t, R *, const X *, Y))
43 {
44  int nd = std::max (x.ndims (), y.ndims ());
45  dim_vector dvx = x.dims ().redim (nd);
46  dim_vector dvy = y.dims ().redim (nd);
47 
48  // Construct the result dimensions.
49  dim_vector dvr;
50  dvr.resize (nd);
51  for (int i = 0; i < nd; i++)
52  {
53  octave_idx_type xk = dvx(i);
54  octave_idx_type yk = dvy(i);
55  if (xk == 1)
56  dvr(i) = yk;
57  else if (yk == 1 || xk == yk)
58  dvr(i) = xk;
59  else
61  ("bsxfun: nonconformant dimensions: %s and %s",
62  x.dims ().str ().c_str (), y.dims ().str ().c_str ());
63  }
64 
65  Array<R> retval (dvr);
66 
67  const X *xvec = x.fortran_vec ();
68  const Y *yvec = y.fortran_vec ();
69  R *rvec = retval.fortran_vec ();
70 
71  // Fold the common leading dimensions.
72  octave_idx_type start, ldr = 1;
73  for (start = 0; start < nd; start++)
74  {
75  if (dvx(start) != dvy(start))
76  break;
77  ldr *= dvr(start);
78  }
79 
80  if (retval.isempty ())
81  ; // do nothing
82  else if (start == nd)
83  op_vv (retval.numel (), rvec, xvec, yvec);
84  else
85  {
86  // Determine the type of the low-level loop.
87  bool xsing = false;
88  bool ysing = false;
89  if (ldr == 1)
90  {
91  xsing = dvx(start) == 1;
92  ysing = dvy(start) == 1;
93  if (xsing || ysing)
94  {
95  ldr *= dvx(start) * dvy(start);
96  start++;
97  }
98  }
99  dim_vector cdvx = dvx.cumulative ();
100  dim_vector cdvy = dvy.cumulative ();
101  // Nullify singleton dims to achieve a spread effect.
102  for (int i = std::max (start, octave_idx_type (1)); i < nd; i++)
103  {
104  if (dvx(i) == 1)
105  cdvx(i-1) = 0;
106  if (dvy(i) == 1)
107  cdvy(i-1) = 0;
108  }
109 
110  octave_idx_type niter = dvr.numel (start);
111  // The index array.
113  for (octave_idx_type iter = 0; iter < niter; iter++)
114  {
115  octave_quit ();
116 
117  // Compute indices.
118  // FIXME: performance impact noticeable?
119  octave_idx_type xidx = cdvx.cum_compute_index (idx);
120  octave_idx_type yidx = cdvy.cum_compute_index (idx);
121  octave_idx_type ridx = dvr.compute_index (idx);
122 
123  // Apply the low-level loop.
124  if (xsing)
125  op_sv (ldr, rvec + ridx, xvec[xidx], yvec + yidx);
126  else if (ysing)
127  op_vs (ldr, rvec + ridx, xvec + xidx, yvec[yidx]);
128  else
129  op_vv (ldr, rvec + ridx, xvec + xidx, yvec + yidx);
130 
131  dvr.increment_index (idx + start, start);
132  }
133  }
134 
135  return retval;
136 }
137 
138 template <typename R, typename X>
139 void
141  void (*op_vv) (size_t, R *, const X *),
142  void (*op_vs) (size_t, R *, X))
143 {
144  dim_vector dvr = r.dims ();
145  dim_vector dvx = x.dims ();
146  octave_idx_type nd = r.ndims ();
147  dvx.redim (nd);
148 
149  const X *xvec = x.fortran_vec ();
150  R *rvec = r.fortran_vec ();
151 
152  // Fold the common leading dimensions.
153  octave_idx_type start, ldr = 1;
154  for (start = 0; start < nd; start++)
155  {
156  if (dvr(start) != dvx(start))
157  break;
158  ldr *= dvr(start);
159  }
160 
161  if (r.isempty ())
162  ; // do nothing
163  else if (start == nd)
164  op_vv (r.numel (), rvec, xvec);
165  else
166  {
167  // Determine the type of the low-level loop.
168  bool xsing = false;
169  if (ldr == 1)
170  {
171  xsing = dvx(start) == 1;
172  if (xsing)
173  {
174  ldr *= dvr(start) * dvx(start);
175  start++;
176  }
177  }
178 
179  dim_vector cdvx = dvx.cumulative ();
180  // Nullify singleton dims to achieve a spread effect.
181  for (int i = std::max (start, octave_idx_type (1)); i < nd; i++)
182  {
183  if (dvx(i) == 1)
184  cdvx(i-1) = 0;
185  }
186 
187  octave_idx_type niter = dvr.numel (start);
188  // The index array.
190  for (octave_idx_type iter = 0; iter < niter; iter++)
191  {
192  octave_quit ();
193 
194  // Compute indices.
195  // FIXME: performance impact noticeable?
196  octave_idx_type xidx = cdvx.cum_compute_index (idx);
197  octave_idx_type ridx = dvr.compute_index (idx);
198 
199  // Apply the low-level loop.
200  if (xsing)
201  op_vs (ldr, rvec + ridx, xvec[xidx]);
202  else
203  op_vv (ldr, rvec + ridx, xvec + xidx);
204 
205  dvr.increment_index (idx + start, start);
206  }
207  }
208 }
209 
210 #define BSXFUN_OP_DEF(OP, ARRAY) \
211  ARRAY bsxfun_ ## OP (const ARRAY& x, const ARRAY& y)
212 
213 #define BSXFUN_OP2_DEF(OP, ARRAY, ARRAY1, ARRAY2) \
214  ARRAY bsxfun_ ## OP (const ARRAY1& x, const ARRAY2& y)
215 
216 #define BSXFUN_REL_DEF(OP, ARRAY) \
217  boolNDArray bsxfun_ ## OP (const ARRAY& x, const ARRAY& y)
218 
219 #define BSXFUN_OP_DEF_MXLOOP(OP, ARRAY, LOOP) \
220  BSXFUN_OP_DEF(OP, ARRAY) \
221  { return do_bsxfun_op<ARRAY::element_type, ARRAY::element_type, ARRAY::element_type> \
222  (x, y, LOOP, LOOP, LOOP); }
223 
224 #define BSXFUN_OP2_DEF_MXLOOP(OP, ARRAY, ARRAY1, ARRAY2, LOOP) \
225  BSXFUN_OP2_DEF(OP, ARRAY, ARRAY1, ARRAY2) \
226  { return do_bsxfun_op<ARRAY::element_type, ARRAY1::element_type, ARRAY2::element_type> \
227  (x, y, LOOP, LOOP, LOOP); }
228 
229 #define BSXFUN_REL_DEF_MXLOOP(OP, ARRAY, LOOP) \
230  BSXFUN_REL_DEF(OP, ARRAY) \
231  { return do_bsxfun_op<bool, ARRAY::element_type, ARRAY::element_type> \
232  (x, y, LOOP, LOOP, LOOP); }
233 
234 #define BSXFUN_STDOP_DEFS_MXLOOP(ARRAY) \
235  BSXFUN_OP_DEF_MXLOOP (add, ARRAY, mx_inline_add) \
236  BSXFUN_OP_DEF_MXLOOP (sub, ARRAY, mx_inline_sub) \
237  BSXFUN_OP_DEF_MXLOOP (mul, ARRAY, mx_inline_mul) \
238  BSXFUN_OP_DEF_MXLOOP (div, ARRAY, mx_inline_div) \
239  BSXFUN_OP_DEF_MXLOOP (min, ARRAY, mx_inline_xmin) \
240  BSXFUN_OP_DEF_MXLOOP (max, ARRAY, mx_inline_xmax)
241 
242 #define BSXFUN_STDREL_DEFS_MXLOOP(ARRAY) \
243  BSXFUN_REL_DEF_MXLOOP (eq, ARRAY, mx_inline_eq) \
244  BSXFUN_REL_DEF_MXLOOP (ne, ARRAY, mx_inline_ne) \
245  BSXFUN_REL_DEF_MXLOOP (lt, ARRAY, mx_inline_lt) \
246  BSXFUN_REL_DEF_MXLOOP (le, ARRAY, mx_inline_le) \
247  BSXFUN_REL_DEF_MXLOOP (gt, ARRAY, mx_inline_gt) \
248  BSXFUN_REL_DEF_MXLOOP (ge, ARRAY, mx_inline_ge)
249 
250 //For bsxfun power with mixed integer/float types
251 #define BSXFUN_POW_MIXED_MXLOOP(INT_TYPE) \
252  BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, INT_TYPE, NDArray, mx_inline_pow) \
253  BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, INT_TYPE, FloatNDArray, mx_inline_pow) \
254  BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, NDArray, INT_TYPE, mx_inline_pow) \
255  BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, FloatNDArray, INT_TYPE, mx_inline_pow)
256 
257 #endif
octave_idx_type compute_index(const octave_idx_type *idx) const
Linear index from an index tuple.
Definition: dim-vector.h:487
bool isempty(void) const
Definition: ov.h:529
#define OCTAVE_LOCAL_BUFFER_INIT(T, buf, size, value)
Definition: oct-locbuf.h:47
bool isempty(void) const
Definition: Array.h:565
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
dim_vector cumulative(void) const
Return cumulative dimensions.
Definition: dim-vector.h:519
void resize(int n, int fill_value=0)
Definition: dim-vector.h:310
const T * fortran_vec(void) const
Definition: Array.h:584
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:442
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:504
octave_idx_type numel(const octave_value_list &idx)
Definition: ov.h:412
octave_value retval
Definition: data.cc:6246
dim_vector redim(int n) const
Force certain dimensionality, preserving numel ().
Definition: dim-vector.cc:233
void do_inplace_bsxfun_op(Array< R > &r, const Array< X > &x, void(*op_vv)(size_t, R *, const X *), void(*op_vs)(size_t, R *, X))
Definition: bsxfun-defs.cc:140
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:125
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:227
octave::sys::time start
Definition: graphics.cc:12337
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
Definition: dim-vector.h:362
the element is set to zero In other the statement xample y
Definition: data.cc:5264
for i
Definition: data.cc:5264
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:366
octave_idx_type cum_compute_index(const octave_idx_type *idx) const
Compute a linear index from an index tuple.
Definition: dim-vector.h:534
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x
Array< R > do_bsxfun_op(const Array< X > &x, const Array< Y > &y, void(*op_vv)(size_t, R *, const X *, const Y *), void(*op_sv)(size_t, R *, X, const Y *), void(*op_vs)(size_t, R *, const X *, Y))
Definition: bsxfun-defs.cc:39
int ndims(void) const
Definition: Array.h:590