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