GNU Octave  4.0.0
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-2015 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 #include <algorithm>
28 #include <iostream>
29 
30 #include "dim-vector.h"
31 #include "oct-locbuf.h"
32 #include "lo-error.h"
33 
34 #include "mx-inlines.cc"
35 
36 template <class R, class X, class Y>
38 do_bsxfun_op (const Array<X>& x, const Array<Y>& y,
39  void (*op_vv) (size_t, R *, const X *, const Y *),
40  void (*op_sv) (size_t, R *, X, const Y *),
41  void (*op_vs) (size_t, R *, const X *, Y))
42 {
43  int nd = std::max (x.ndims (), y.ndims ());
44  dim_vector dvx = x.dims ().redim (nd);
45  dim_vector dvy = y.dims ().redim (nd);
46 
47  // Construct the result dimensions.
48  dim_vector dvr;
49  dvr.resize (nd);
50  for (int i = 0; i < nd; i++)
51  {
52  octave_idx_type xk = dvx(i);
53  octave_idx_type yk = dvy(i);
54  if (xk == 1)
55  dvr(i) = yk;
56  else if (yk == 1 || xk == yk)
57  dvr(i) = xk;
58  else
59  {
60  (*current_liboctave_error_handler)
61  ("bsxfun: nonconformant dimensions: %s and %s",
62  x.dims ().str ().c_str (), y.dims ().str ().c_str ());
63  break;
64  }
65  }
66 
67  Array<R> retval (dvr);
68 
69  const X *xvec = x.fortran_vec ();
70  const Y *yvec = y.fortran_vec ();
71  R *rvec = retval.fortran_vec ();
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.is_empty ())
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.
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 
140 template <class R, class X>
141 void
143  void (*op_vv) (size_t, R *, const X *),
144  void (*op_vs) (size_t, R *, X))
145 {
146  dim_vector dvr = r.dims ();
147  dim_vector dvx = x.dims ();
148  octave_idx_type nd = r.ndims ();
149  dvx.redim (nd);
150 
151  const X* xvec = x.fortran_vec ();
152  R* rvec = r.fortran_vec ();
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.is_empty ())
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.
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 #define BSXFUN_OP_DEF(OP, ARRAY) \
213 ARRAY bsxfun_ ## OP (const ARRAY& x, const ARRAY& y)
214 
215 #define BSXFUN_OP2_DEF(OP, ARRAY, ARRAY1, ARRAY2) \
216 ARRAY bsxfun_ ## OP (const ARRAY1& x, const ARRAY2& y)
217 
218 #define BSXFUN_REL_DEF(OP, ARRAY) \
219 boolNDArray bsxfun_ ## OP (const ARRAY& x, const ARRAY& y)
220 
221 #define BSXFUN_OP_DEF_MXLOOP(OP, ARRAY, LOOP) \
222  BSXFUN_OP_DEF(OP, ARRAY) \
223  { return do_bsxfun_op<ARRAY::element_type, ARRAY::element_type, ARRAY::element_type> \
224  (x, y, LOOP, LOOP, LOOP); }
225 
226 #define BSXFUN_OP2_DEF_MXLOOP(OP, ARRAY, ARRAY1, ARRAY2, LOOP) \
227  BSXFUN_OP2_DEF(OP, ARRAY, ARRAY1, ARRAY2) \
228  { return do_bsxfun_op<ARRAY::element_type, ARRAY1::element_type, ARRAY2::element_type> \
229  (x, y, LOOP, LOOP, LOOP); }
230 
231 #define BSXFUN_REL_DEF_MXLOOP(OP, ARRAY, LOOP) \
232  BSXFUN_REL_DEF(OP, ARRAY) \
233  { return do_bsxfun_op<bool, ARRAY::element_type, ARRAY::element_type> \
234  (x, y, LOOP, LOOP, LOOP); }
235 
236 #define BSXFUN_STDOP_DEFS_MXLOOP(ARRAY) \
237  BSXFUN_OP_DEF_MXLOOP (add, ARRAY, mx_inline_add) \
238  BSXFUN_OP_DEF_MXLOOP (sub, ARRAY, mx_inline_sub) \
239  BSXFUN_OP_DEF_MXLOOP (mul, ARRAY, mx_inline_mul) \
240  BSXFUN_OP_DEF_MXLOOP (div, ARRAY, mx_inline_div) \
241  BSXFUN_OP_DEF_MXLOOP (min, ARRAY, mx_inline_xmin) \
242  BSXFUN_OP_DEF_MXLOOP (max, ARRAY, mx_inline_xmax)
243 
244 #define BSXFUN_STDREL_DEFS_MXLOOP(ARRAY) \
245  BSXFUN_REL_DEF_MXLOOP (eq, ARRAY, mx_inline_eq) \
246  BSXFUN_REL_DEF_MXLOOP (ne, ARRAY, mx_inline_ne) \
247  BSXFUN_REL_DEF_MXLOOP (lt, ARRAY, mx_inline_lt) \
248  BSXFUN_REL_DEF_MXLOOP (le, ARRAY, mx_inline_le) \
249  BSXFUN_REL_DEF_MXLOOP (gt, ARRAY, mx_inline_gt) \
250  BSXFUN_REL_DEF_MXLOOP (ge, ARRAY, mx_inline_ge)
251 
252 //For bsxfun power with mixed integer/float types
253 #define BSXFUN_POW_MIXED_MXLOOP(INT_TYPE) \
254  BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, INT_TYPE, NDArray, mx_inline_pow) \
255  BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, INT_TYPE, FloatNDArray, mx_inline_pow)\
256  BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, NDArray, INT_TYPE, mx_inline_pow) \
257  BSXFUN_OP2_DEF_MXLOOP (pow, INT_TYPE, FloatNDArray, INT_TYPE, mx_inline_pow)
258 
259 #endif
bool is_empty(void) const
Definition: Array.h:472
std::string str(char sep= 'x') const
Definition: dim-vector.cc:63
#define OCTAVE_LOCAL_BUFFER_INIT(T, buf, size, value)
Definition: oct-locbuf.h:206
int ndims(void) const
Definition: Array.h:487
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:275
dim_vector cumulative(void) const
Return cumulative dimensions.
Definition: dim-vector.h:489
void resize(int n, int fill_value=0)
Definition: dim-vector.h:287
octave_idx_type cum_compute_index(const octave_idx_type *idx) const
Compute a linear index from an index tuple.
Definition: dim-vector.h:504
int increment_index(octave_idx_type *idx, int start=0) const
Definition: dim-vector.h:474
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
Definition: dim-vector.h:361
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:337
dim_vector redim(int n) const
Definition: dim-vector.cc:266
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:142
Handles the reference counting for all the derived classes.
Definition: Array.h:45
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:233
octave_idx_type compute_index(const octave_idx_type *idx) const
Compute a linear index from an index tuple.
Definition: dim-vector.h:448
const T * fortran_vec(void) const
Definition: Array.h:481
F77_RET_T const double * 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:38