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