GNU Octave  3.8.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
MArray.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1993-2013 John W. Eaton
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 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include "MArray.h"
29 #include "Array-util.h"
30 #include "lo-error.h"
31 
32 #include "MArray-defs.h"
33 #include "mx-inlines.cc"
34 
35 template <class T>
37 {
38  T *array;
39  T val;
40  _idxadds_helper (T *a, T v) : array (a), val (v) { }
42  { array[i] += val; }
43 };
44 
45 template <class T>
47 {
48  T *array;
49  const T *vals;
50  _idxadda_helper (T *a, const T *v) : array (a), vals (v) { }
52  { array[i] += *vals++; }
53 };
54 
55 template <class T>
56 void
57 MArray<T>::idx_add (const idx_vector& idx, T val)
58 {
59  octave_idx_type n = this->length ();
60  octave_idx_type ext = idx.extent (n);
61  if (ext > n)
62  {
63  this->resize1 (ext);
64  n = ext;
65  }
66 
67  octave_quit ();
68 
69  octave_idx_type len = idx.length (n);
70  idx.loop (len, _idxadds_helper<T> (this->fortran_vec (), val));
71 }
72 
73 template <class T>
74 void
75 MArray<T>::idx_add (const idx_vector& idx, const MArray<T>& vals)
76 {
77  octave_idx_type n = this->length ();
78  octave_idx_type ext = idx.extent (n);
79  if (ext > n)
80  {
81  this->resize1 (ext);
82  n = ext;
83  }
84 
85  octave_quit ();
86 
87  octave_idx_type len = std::min (idx.length (n), vals.length ());
88  idx.loop (len, _idxadda_helper<T> (this->fortran_vec (), vals.data ()));
89 }
90 
92  typename ref_param<T>::type)>
94 {
95  T *array;
96  const T *vals;
97  _idxbinop_helper (T *a, const T *v) : array (a), vals (v) { }
98  void operator () (octave_idx_type i)
99  { array[i] = op (array[i], *vals++); }
100 };
101 
102 template <class T>
103 void
104 MArray<T>::idx_min (const idx_vector& idx, const MArray<T>& vals)
105 {
106  octave_idx_type n = this->length ();
107  octave_idx_type ext = idx.extent (n);
108  if (ext > n)
109  {
110  this->resize1 (ext);
111  n = ext;
112  }
113 
114  octave_quit ();
115 
116  octave_idx_type len = std::min (idx.length (n), vals.length ());
117  idx.loop (len, _idxbinop_helper<T, xmin> (this->fortran_vec (),
118  vals.data ()));
119 }
120 
121 template <class T>
122 void
123 MArray<T>::idx_max (const idx_vector& idx, const MArray<T>& vals)
124 {
125  octave_idx_type n = this->length ();
126  octave_idx_type ext = idx.extent (n);
127  if (ext > n)
128  {
129  this->resize1 (ext);
130  n = ext;
131  }
132 
133  octave_quit ();
134 
135  octave_idx_type len = std::min (idx.length (n), vals.length ());
136  idx.loop (len, _idxbinop_helper<T, xmax> (this->fortran_vec (),
137  vals.data ()));
138 }
139 
140 #include <iostream>
141 
142 template <class T>
143 void MArray<T>::idx_add_nd (const idx_vector& idx, const MArray<T>& vals,
144  int dim)
145 {
146  int nd = std::max (this->ndims (), vals.ndims ());
147  if (dim < 0)
148  dim = vals.dims ().first_non_singleton ();
149  else if (dim > nd)
150  nd = dim;
151 
152  // Check dimensions.
153  dim_vector ddv = Array<T>::dims ().redim (nd);
154  dim_vector sdv = vals.dims ().redim (nd);
155 
156  octave_idx_type ext = idx.extent (ddv (dim));
157 
158  if (ext > ddv(dim))
159  {
160  ddv(dim) = ext;
161  Array<T>::resize (ddv);
162  ext = ddv(dim);
163  }
164 
165  octave_idx_type l,n,u,ns;
166  get_extent_triplet (ddv, dim, l, n, u);
167  ns = sdv(dim);
168 
169  sdv(dim) = ddv(dim) = 0;
170  if (ddv != sdv)
171  (*current_liboctave_error_handler)
172  ("accumdim: dimension mismatch");
173 
174  T *dst = Array<T>::fortran_vec ();
175  const T *src = vals.data ();
176  octave_idx_type len = idx.length (ns);
177 
178  if (l == 1)
179  {
180  for (octave_idx_type j = 0; j < u; j++)
181  {
182  octave_quit ();
183 
184  idx.loop (len, _idxadda_helper<T> (dst + j*n, src + j*ns));
185  }
186  }
187  else
188  {
189  for (octave_idx_type j = 0; j < u; j++)
190  {
191  octave_quit ();
192  for (octave_idx_type i = 0; i < len; i++)
193  {
194  octave_idx_type k = idx(i);
195 
196  mx_inline_add2 (l, dst + l*k, src + l*i);
197  }
198 
199  dst += l*n;
200  src += l*ns;
201  }
202  }
203 }
204 
205 // N-dimensional array with math ops.
206 template <class T>
207 void
209 {
210  if (Array<T>::is_shared ())
211  *this = - *this;
212  else
213  do_mx_inplace_op<T> (*this, mx_inline_uminus2);
214 }
215 
216 // Element by element MArray by scalar ops.
217 
218 template <class T>
219 MArray<T>&
220 operator += (MArray<T>& a, const T& s)
221 {
222  if (a.is_shared ())
223  a = a + s;
224  else
225  do_ms_inplace_op<T, T> (a, s, mx_inline_add2);
226  return a;
227 }
228 
229 template <class T>
230 MArray<T>&
231 operator -= (MArray<T>& a, const T& s)
232 {
233  if (a.is_shared ())
234  a = a - s;
235  else
236  do_ms_inplace_op<T, T> (a, s, mx_inline_sub2);
237  return a;
238 }
239 
240 template <class T>
241 MArray<T>&
242 operator *= (MArray<T>& a, const T& s)
243 {
244  if (a.is_shared ())
245  a = a * s;
246  else
247  do_ms_inplace_op<T, T> (a, s, mx_inline_mul2);
248  return a;
249 }
250 
251 template <class T>
252 MArray<T>&
253 operator /= (MArray<T>& a, const T& s)
254 {
255  if (a.is_shared ())
256  a = a / s;
257  else
258  do_ms_inplace_op<T, T> (a, s, mx_inline_div2);
259  return a;
260 }
261 
262 // Element by element MArray by MArray ops.
263 
264 template <class T>
265 MArray<T>&
267 {
268  if (a.is_shared ())
269  a = a + b;
270  else
271  do_mm_inplace_op<T, T> (a, b, mx_inline_add2, mx_inline_add2, "+=");
272  return a;
273 }
274 
275 template <class T>
276 MArray<T>&
278 {
279  if (a.is_shared ())
280  a = a - b;
281  else
282  do_mm_inplace_op<T, T> (a, b, mx_inline_sub2, mx_inline_sub2, "-=");
283  return a;
284 }
285 
286 
287 template <class T>
288 MArray<T>&
290 {
291  if (a.is_shared ())
292  return a = product (a, b);
293  else
294  do_mm_inplace_op<T, T> (a, b, mx_inline_mul2, mx_inline_mul2, ".*=");
295  return a;
296 }
297 
298 template <class T>
299 MArray<T>&
301 {
302  if (a.is_shared ())
303  return a = quotient (a, b);
304  else
305  do_mm_inplace_op<T, T> (a, b, mx_inline_div2, mx_inline_div2, "./=");
306  return a;
307 }
308 
309 // Element by element MArray by scalar ops.
310 
311 #define MARRAY_NDS_OP(OP, FN) \
312  template <class T> \
313  MArray<T> \
314  operator OP (const MArray<T>& a, const T& s) \
315  { \
316  return do_ms_binary_op<T, T, T> (a, s, FN); \
317  }
318 
323 
324 // Element by element scalar by MArray ops.
325 
326 #define MARRAY_SND_OP(OP, FN) \
327  template <class T> \
328  MArray<T> \
329  operator OP (const T& s, const MArray<T>& a) \
330  { \
331  return do_sm_binary_op<T, T, T> (s, a, FN); \
332  }
333 
338 
339 // Element by element MArray by MArray ops.
340 
341 #define MARRAY_NDND_OP(FCN, OP, FN) \
342  template <class T> \
343  MArray<T> \
344  FCN (const MArray<T>& a, const MArray<T>& b) \
345  { \
346  return do_mm_binary_op<T, T, T> (a, b, FN, FN, FN, #FCN); \
347  }
348 
353 
354 template <class T>
355 MArray<T>
356 operator + (const MArray<T>& a)
357 {
358  return a;
359 }
360 
361 template <class T>
362 MArray<T>
364 {
365  return do_mx_unary_op<T, T> (a, mx_inline_uminus);
366 }