1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup subview_cube
20 //! @{
21
22
23 template<typename eT>
24 inline
~subview_cube()25 subview_cube<eT>::~subview_cube()
26 {
27 arma_extra_debug_sigprint_this(this);
28 }
29
30
31
32 template<typename eT>
33 arma_inline
subview_cube(const Cube<eT> & in_m,const uword in_row1,const uword in_col1,const uword in_slice1,const uword in_n_rows,const uword in_n_cols,const uword in_n_slices)34 subview_cube<eT>::subview_cube
35 (
36 const Cube<eT>& in_m,
37 const uword in_row1,
38 const uword in_col1,
39 const uword in_slice1,
40 const uword in_n_rows,
41 const uword in_n_cols,
42 const uword in_n_slices
43 )
44 : m (in_m)
45 , aux_row1 (in_row1)
46 , aux_col1 (in_col1)
47 , aux_slice1 (in_slice1)
48 , n_rows (in_n_rows)
49 , n_cols (in_n_cols)
50 , n_elem_slice(in_n_rows * in_n_cols)
51 , n_slices (in_n_slices)
52 , n_elem (n_elem_slice * in_n_slices)
53 {
54 arma_extra_debug_sigprint_this(this);
55 }
56
57
58
59 template<typename eT>
60 inline
subview_cube(const subview_cube<eT> & in)61 subview_cube<eT>::subview_cube(const subview_cube<eT>& in)
62 : m (in.m )
63 , aux_row1 (in.aux_row1 )
64 , aux_col1 (in.aux_col1 )
65 , aux_slice1 (in.aux_slice1 )
66 , n_rows (in.n_rows )
67 , n_cols (in.n_cols )
68 , n_elem_slice(in.n_elem_slice)
69 , n_slices (in.n_slices )
70 , n_elem (in.n_elem )
71 {
72 arma_extra_debug_sigprint(arma_str::format("this = %x in = %x") % this % &in);
73 }
74
75
76
77 template<typename eT>
78 inline
subview_cube(subview_cube<eT> && in)79 subview_cube<eT>::subview_cube(subview_cube<eT>&& in)
80 : m (in.m )
81 , aux_row1 (in.aux_row1 )
82 , aux_col1 (in.aux_col1 )
83 , aux_slice1 (in.aux_slice1 )
84 , n_rows (in.n_rows )
85 , n_cols (in.n_cols )
86 , n_elem_slice(in.n_elem_slice)
87 , n_slices (in.n_slices )
88 , n_elem (in.n_elem )
89 {
90 arma_extra_debug_sigprint(arma_str::format("this = %x in = %x") % this % &in);
91
92 // for paranoia
93
94 access::rw(in.aux_row1 ) = 0;
95 access::rw(in.aux_col1 ) = 0;
96 access::rw(in.aux_slice1 ) = 0;
97 access::rw(in.n_rows ) = 0;
98 access::rw(in.n_cols ) = 0;
99 access::rw(in.n_elem_slice) = 0;
100 access::rw(in.n_slices ) = 0;
101 access::rw(in.n_elem ) = 0;
102 }
103
104
105
106 template<typename eT>
107 template<typename op_type>
108 inline
109 void
inplace_op(const eT val)110 subview_cube<eT>::inplace_op(const eT val)
111 {
112 arma_extra_debug_sigprint();
113
114 subview_cube<eT>& t = *this;
115
116 const uword t_n_rows = t.n_rows;
117 const uword t_n_cols = t.n_cols;
118 const uword t_n_slices = t.n_slices;
119
120 for(uword s=0; s < t_n_slices; ++s)
121 for(uword c=0; c < t_n_cols; ++c)
122 {
123 if(is_same_type<op_type, op_internal_plus >::yes) { arrayops::inplace_plus ( slice_colptr(s,c), val, t_n_rows ); }
124 if(is_same_type<op_type, op_internal_minus>::yes) { arrayops::inplace_minus( slice_colptr(s,c), val, t_n_rows ); }
125 if(is_same_type<op_type, op_internal_schur>::yes) { arrayops::inplace_mul ( slice_colptr(s,c), val, t_n_rows ); }
126 if(is_same_type<op_type, op_internal_div >::yes) { arrayops::inplace_div ( slice_colptr(s,c), val, t_n_rows ); }
127 }
128 }
129
130
131
132
133
134
135 template<typename eT>
136 template<typename op_type, typename T1>
137 inline
138 void
inplace_op(const BaseCube<eT,T1> & in,const char * identifier)139 subview_cube<eT>::inplace_op(const BaseCube<eT,T1>& in, const char* identifier)
140 {
141 arma_extra_debug_sigprint();
142
143 const ProxyCube<T1> P(in.get_ref());
144
145 subview_cube<eT>& t = *this;
146
147 const uword t_n_rows = t.n_rows;
148 const uword t_n_cols = t.n_cols;
149 const uword t_n_slices = t.n_slices;
150
151 arma_debug_assert_same_size(t, P, identifier);
152
153 const bool use_mp = arma_config::openmp && ProxyCube<T1>::use_mp && mp_gate<eT>::eval(t.n_elem);
154 const bool has_overlap = P.has_overlap(t);
155
156 if(has_overlap) { arma_extra_debug_print("aliasing or overlap detected"); }
157
158 if( (is_Cube<typename ProxyCube<T1>::stored_type>::value) || (use_mp) || (has_overlap) )
159 {
160 const unwrap_cube_check<typename ProxyCube<T1>::stored_type> tmp(P.Q, has_overlap);
161 const Cube<eT>& B = tmp.M;
162
163 if( (is_same_type<op_type, op_internal_equ>::yes) && (t.aux_row1 == 0) && (t_n_rows == t.m.n_rows) )
164 {
165 for(uword s=0; s < t_n_slices; ++s)
166 {
167 arrayops::copy( t.slice_colptr(s,0), B.slice_colptr(s,0), t.n_elem_slice );
168 }
169 }
170 else
171 {
172 for(uword s=0; s < t_n_slices; ++s)
173 for(uword c=0; c < t_n_cols; ++c)
174 {
175 if(is_same_type<op_type, op_internal_equ >::yes) { arrayops::copy ( t.slice_colptr(s,c), B.slice_colptr(s,c), t_n_rows ); }
176 if(is_same_type<op_type, op_internal_plus >::yes) { arrayops::inplace_plus ( t.slice_colptr(s,c), B.slice_colptr(s,c), t_n_rows ); }
177 if(is_same_type<op_type, op_internal_minus>::yes) { arrayops::inplace_minus( t.slice_colptr(s,c), B.slice_colptr(s,c), t_n_rows ); }
178 if(is_same_type<op_type, op_internal_schur>::yes) { arrayops::inplace_mul ( t.slice_colptr(s,c), B.slice_colptr(s,c), t_n_rows ); }
179 if(is_same_type<op_type, op_internal_div >::yes) { arrayops::inplace_div ( t.slice_colptr(s,c), B.slice_colptr(s,c), t_n_rows ); }
180 }
181 }
182 }
183 else // use the Proxy
184 {
185 if(ProxyCube<T1>::use_at)
186 {
187 for(uword s=0; s < t_n_slices; ++s)
188 for(uword c=0; c < t_n_cols; ++c)
189 {
190 eT* t_col_data = t.slice_colptr(s,c);
191
192 for(uword r=0; r < t_n_rows; ++r)
193 {
194 const eT tmp = P.at(r,c,s);
195
196 if(is_same_type<op_type, op_internal_equ >::yes) { (*t_col_data) = tmp; t_col_data++; }
197 if(is_same_type<op_type, op_internal_plus >::yes) { (*t_col_data) += tmp; t_col_data++; }
198 if(is_same_type<op_type, op_internal_minus>::yes) { (*t_col_data) -= tmp; t_col_data++; }
199 if(is_same_type<op_type, op_internal_schur>::yes) { (*t_col_data) *= tmp; t_col_data++; }
200 if(is_same_type<op_type, op_internal_div >::yes) { (*t_col_data) /= tmp; t_col_data++; }
201 }
202 }
203 }
204 else
205 {
206 typename ProxyCube<T1>::ea_type Pea = P.get_ea();
207
208 uword count = 0;
209
210 for(uword s=0; s < t_n_slices; ++s)
211 for(uword c=0; c < t_n_cols; ++c)
212 {
213 eT* t_col_data = t.slice_colptr(s,c);
214
215 for(uword r=0; r < t_n_rows; ++r)
216 {
217 const eT tmp = Pea[count]; count++;
218
219 if(is_same_type<op_type, op_internal_equ >::yes) { (*t_col_data) = tmp; t_col_data++; }
220 if(is_same_type<op_type, op_internal_plus >::yes) { (*t_col_data) += tmp; t_col_data++; }
221 if(is_same_type<op_type, op_internal_minus>::yes) { (*t_col_data) -= tmp; t_col_data++; }
222 if(is_same_type<op_type, op_internal_schur>::yes) { (*t_col_data) *= tmp; t_col_data++; }
223 if(is_same_type<op_type, op_internal_div >::yes) { (*t_col_data) /= tmp; t_col_data++; }
224 }
225 }
226 }
227 }
228 }
229
230
231
232 template<typename eT>
233 template<typename op_type>
234 inline
235 void
inplace_op(const subview_cube<eT> & x,const char * identifier)236 subview_cube<eT>::inplace_op(const subview_cube<eT>& x, const char* identifier)
237 {
238 arma_extra_debug_sigprint();
239
240 if(check_overlap(x))
241 {
242 const Cube<eT> tmp(x);
243
244 if(is_same_type<op_type, op_internal_equ >::yes) { (*this).operator= (tmp); }
245 if(is_same_type<op_type, op_internal_plus >::yes) { (*this).operator+=(tmp); }
246 if(is_same_type<op_type, op_internal_minus>::yes) { (*this).operator-=(tmp); }
247 if(is_same_type<op_type, op_internal_schur>::yes) { (*this).operator%=(tmp); }
248 if(is_same_type<op_type, op_internal_div >::yes) { (*this).operator/=(tmp); }
249
250 return;
251 }
252
253 subview_cube<eT>& t = *this;
254
255 arma_debug_assert_same_size(t, x, identifier);
256
257 const uword t_n_rows = t.n_rows;
258 const uword t_n_cols = t.n_cols;
259 const uword t_n_slices = t.n_slices;
260
261 for(uword s=0; s < t_n_slices; ++s)
262 for(uword c=0; c < t_n_cols; ++c)
263 {
264 if(is_same_type<op_type, op_internal_equ >::yes) { arrayops::copy ( t.slice_colptr(s,c), x.slice_colptr(s,c), t_n_rows ); }
265 if(is_same_type<op_type, op_internal_plus >::yes) { arrayops::inplace_plus ( t.slice_colptr(s,c), x.slice_colptr(s,c), t_n_rows ); }
266 if(is_same_type<op_type, op_internal_minus>::yes) { arrayops::inplace_minus( t.slice_colptr(s,c), x.slice_colptr(s,c), t_n_rows ); }
267 if(is_same_type<op_type, op_internal_schur>::yes) { arrayops::inplace_mul ( t.slice_colptr(s,c), x.slice_colptr(s,c), t_n_rows ); }
268 if(is_same_type<op_type, op_internal_div >::yes) { arrayops::inplace_div ( t.slice_colptr(s,c), x.slice_colptr(s,c), t_n_rows ); }
269 }
270 }
271
272
273
274 template<typename eT>
275 inline
276 void
operator =(const eT val)277 subview_cube<eT>::operator= (const eT val)
278 {
279 arma_extra_debug_sigprint();
280
281 if(n_elem != 1)
282 {
283 arma_debug_assert_same_size(n_rows, n_cols, n_slices, 1, 1, 1, "copy into subcube");
284 }
285
286 Cube<eT>& Q = const_cast< Cube<eT>& >(m);
287
288 Q.at(aux_row1, aux_col1, aux_slice1) = val;
289 }
290
291
292
293 template<typename eT>
294 inline
295 void
operator +=(const eT val)296 subview_cube<eT>::operator+= (const eT val)
297 {
298 arma_extra_debug_sigprint();
299
300 inplace_op<op_internal_plus>(val);
301 }
302
303
304
305 template<typename eT>
306 inline
307 void
operator -=(const eT val)308 subview_cube<eT>::operator-= (const eT val)
309 {
310 arma_extra_debug_sigprint();
311
312 inplace_op<op_internal_minus>(val);
313 }
314
315
316
317 template<typename eT>
318 inline
319 void
operator *=(const eT val)320 subview_cube<eT>::operator*= (const eT val)
321 {
322 arma_extra_debug_sigprint();
323
324 inplace_op<op_internal_schur>(val);
325 }
326
327
328
329 template<typename eT>
330 inline
331 void
operator /=(const eT val)332 subview_cube<eT>::operator/= (const eT val)
333 {
334 arma_extra_debug_sigprint();
335
336 inplace_op<op_internal_div>(val);
337 }
338
339
340
341 template<typename eT>
342 template<typename T1>
343 inline
344 void
operator =(const BaseCube<eT,T1> & in)345 subview_cube<eT>::operator= (const BaseCube<eT,T1>& in)
346 {
347 arma_extra_debug_sigprint();
348
349 inplace_op<op_internal_equ>(in, "copy into subcube");
350 }
351
352
353
354 template<typename eT>
355 template<typename T1>
356 inline
357 void
operator +=(const BaseCube<eT,T1> & in)358 subview_cube<eT>::operator+= (const BaseCube<eT,T1>& in)
359 {
360 arma_extra_debug_sigprint();
361
362 inplace_op<op_internal_plus>(in, "addition");
363 }
364
365
366
367 template<typename eT>
368 template<typename T1>
369 inline
370 void
operator -=(const BaseCube<eT,T1> & in)371 subview_cube<eT>::operator-= (const BaseCube<eT,T1>& in)
372 {
373 arma_extra_debug_sigprint();
374
375 inplace_op<op_internal_minus>(in, "subtraction");
376 }
377
378
379
380 template<typename eT>
381 template<typename T1>
382 inline
383 void
operator %=(const BaseCube<eT,T1> & in)384 subview_cube<eT>::operator%= (const BaseCube<eT,T1>& in)
385 {
386 arma_extra_debug_sigprint();
387
388 inplace_op<op_internal_schur>(in, "element-wise multiplication");
389 }
390
391
392
393 template<typename eT>
394 template<typename T1>
395 inline
396 void
operator /=(const BaseCube<eT,T1> & in)397 subview_cube<eT>::operator/= (const BaseCube<eT,T1>& in)
398 {
399 arma_extra_debug_sigprint();
400
401 inplace_op<op_internal_div>(in, "element-wise division");
402 }
403
404
405
406 //! x.subcube(...) = y.subcube(...)
407 template<typename eT>
408 inline
409 void
operator =(const subview_cube<eT> & x)410 subview_cube<eT>::operator= (const subview_cube<eT>& x)
411 {
412 arma_extra_debug_sigprint();
413
414 inplace_op<op_internal_equ>(x, "copy into subcube");
415 }
416
417
418
419 template<typename eT>
420 inline
421 void
operator +=(const subview_cube<eT> & x)422 subview_cube<eT>::operator+= (const subview_cube<eT>& x)
423 {
424 arma_extra_debug_sigprint();
425
426 inplace_op<op_internal_plus>(x, "addition");
427 }
428
429
430
431 template<typename eT>
432 inline
433 void
operator -=(const subview_cube<eT> & x)434 subview_cube<eT>::operator-= (const subview_cube<eT>& x)
435 {
436 arma_extra_debug_sigprint();
437
438 inplace_op<op_internal_minus>(x, "subtraction");
439 }
440
441
442
443 template<typename eT>
444 inline
445 void
operator %=(const subview_cube<eT> & x)446 subview_cube<eT>::operator%= (const subview_cube<eT>& x)
447 {
448 arma_extra_debug_sigprint();
449
450 inplace_op<op_internal_schur>(x, "element-wise multiplication");
451 }
452
453
454
455 template<typename eT>
456 inline
457 void
operator /=(const subview_cube<eT> & x)458 subview_cube<eT>::operator/= (const subview_cube<eT>& x)
459 {
460 arma_extra_debug_sigprint();
461
462 inplace_op<op_internal_div>(x, "element-wise division");
463 }
464
465
466
467 template<typename eT>
468 template<typename T1>
469 inline
470 void
operator =(const Base<eT,T1> & in)471 subview_cube<eT>::operator= (const Base<eT,T1>& in)
472 {
473 arma_extra_debug_sigprint();
474
475 const quasi_unwrap<T1> tmp(in.get_ref());
476
477 const Mat<eT>& x = tmp.M;
478 subview_cube<eT>& t = *this;
479
480 const uword t_n_rows = t.n_rows;
481 const uword t_n_cols = t.n_cols;
482 const uword t_n_slices = t.n_slices;
483
484 const uword x_n_rows = x.n_rows;
485 const uword x_n_cols = x.n_cols;
486
487 if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
488 {
489 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
490
491 const uword t_aux_row1 = t.aux_row1;
492 const uword t_aux_col1 = t.aux_col1;
493 const uword t_aux_slice1 = t.aux_slice1;
494
495 const eT* x_mem = x.memptr();
496
497 uword i,j;
498 for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
499 {
500 const eT tmp_i = x_mem[i];
501 const eT tmp_j = x_mem[j];
502
503 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) = tmp_i;
504 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) = tmp_j;
505 }
506
507 if(i < t_n_slices)
508 {
509 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) = x_mem[i];
510 }
511 }
512 else
513 if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
514 {
515 // interpret the matrix as a cube with one slice
516
517 for(uword col = 0; col < t_n_cols; ++col)
518 {
519 arrayops::copy( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
520 }
521 }
522 else
523 if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
524 {
525 for(uword i=0; i < t_n_slices; ++i)
526 {
527 arrayops::copy( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
528 }
529 }
530 else
531 if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
532 {
533 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
534
535 const uword t_aux_row1 = t.aux_row1;
536 const uword t_aux_col1 = t.aux_col1;
537 const uword t_aux_slice1 = t.aux_slice1;
538
539 for(uword slice=0; slice < t_n_slices; ++slice)
540 {
541 const uword mod_slice = t_aux_slice1 + slice;
542
543 const eT* x_colptr = x.colptr(slice);
544
545 uword i,j;
546 for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
547 {
548 const eT tmp_i = x_colptr[i];
549 const eT tmp_j = x_colptr[j];
550
551 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) = tmp_i;
552 Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) = tmp_j;
553 }
554
555 if(i < t_n_cols)
556 {
557 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) = x_colptr[i];
558 }
559 }
560 }
561 else
562 {
563 if(arma_config::debug)
564 {
565 arma_stop_logic_error( arma_incompat_size_string(t, x, "copy into subcube") );
566 }
567 }
568 }
569
570
571
572 template<typename eT>
573 template<typename T1>
574 inline
575 void
operator +=(const Base<eT,T1> & in)576 subview_cube<eT>::operator+= (const Base<eT,T1>& in)
577 {
578 arma_extra_debug_sigprint();
579
580 const quasi_unwrap<T1> tmp(in.get_ref());
581
582 const Mat<eT>& x = tmp.M;
583 subview_cube<eT>& t = *this;
584
585 const uword t_n_rows = t.n_rows;
586 const uword t_n_cols = t.n_cols;
587 const uword t_n_slices = t.n_slices;
588
589 const uword x_n_rows = x.n_rows;
590 const uword x_n_cols = x.n_cols;
591
592 if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
593 {
594 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
595
596 const uword t_aux_row1 = t.aux_row1;
597 const uword t_aux_col1 = t.aux_col1;
598 const uword t_aux_slice1 = t.aux_slice1;
599
600 const eT* x_mem = x.memptr();
601
602 uword i,j;
603 for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
604 {
605 const eT tmp_i = x_mem[i];
606 const eT tmp_j = x_mem[j];
607
608 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) += tmp_i;
609 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) += tmp_j;
610 }
611
612 if(i < t_n_slices)
613 {
614 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) += x_mem[i];
615 }
616 }
617 else
618 if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
619 {
620 for(uword col = 0; col < t_n_cols; ++col)
621 {
622 arrayops::inplace_plus( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
623 }
624 }
625 else
626 if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
627 {
628 for(uword i=0; i < t_n_slices; ++i)
629 {
630 arrayops::inplace_plus( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
631 }
632 }
633 else
634 if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
635 {
636 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
637
638 const uword t_aux_row1 = t.aux_row1;
639 const uword t_aux_col1 = t.aux_col1;
640 const uword t_aux_slice1 = t.aux_slice1;
641
642 for(uword slice=0; slice < t_n_slices; ++slice)
643 {
644 const uword mod_slice = t_aux_slice1 + slice;
645
646 const eT* x_colptr = x.colptr(slice);
647
648 uword i,j;
649 for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
650 {
651 const eT tmp_i = x_colptr[i];
652 const eT tmp_j = x_colptr[j];
653
654 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) += tmp_i;
655 Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) += tmp_j;
656 }
657
658 if(i < t_n_cols)
659 {
660 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) += x_colptr[i];
661 }
662 }
663 }
664 else
665 {
666 if(arma_config::debug)
667 {
668 arma_stop_logic_error( arma_incompat_size_string(t, x, "addition") );
669 }
670 }
671 }
672
673
674
675 template<typename eT>
676 template<typename T1>
677 inline
678 void
operator -=(const Base<eT,T1> & in)679 subview_cube<eT>::operator-= (const Base<eT,T1>& in)
680 {
681 arma_extra_debug_sigprint();
682
683 const quasi_unwrap<T1> tmp(in.get_ref());
684
685 const Mat<eT>& x = tmp.M;
686 subview_cube<eT>& t = *this;
687
688 const uword t_n_rows = t.n_rows;
689 const uword t_n_cols = t.n_cols;
690 const uword t_n_slices = t.n_slices;
691
692 const uword x_n_rows = x.n_rows;
693 const uword x_n_cols = x.n_cols;
694
695 if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
696 {
697 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
698
699 const uword t_aux_row1 = t.aux_row1;
700 const uword t_aux_col1 = t.aux_col1;
701 const uword t_aux_slice1 = t.aux_slice1;
702
703 const eT* x_mem = x.memptr();
704
705 uword i,j;
706 for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
707 {
708 const eT tmp_i = x_mem[i];
709 const eT tmp_j = x_mem[j];
710
711 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) -= tmp_i;
712 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) -= tmp_j;
713 }
714
715 if(i < t_n_slices)
716 {
717 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) -= x_mem[i];
718 }
719 }
720 else
721 if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
722 {
723 for(uword col = 0; col < t_n_cols; ++col)
724 {
725 arrayops::inplace_minus( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
726 }
727 }
728 else
729 if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
730 {
731 for(uword i=0; i < t_n_slices; ++i)
732 {
733 arrayops::inplace_minus( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
734 }
735 }
736 else
737 if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
738 {
739 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
740
741 const uword t_aux_row1 = t.aux_row1;
742 const uword t_aux_col1 = t.aux_col1;
743 const uword t_aux_slice1 = t.aux_slice1;
744
745 for(uword slice=0; slice < t_n_slices; ++slice)
746 {
747 const uword mod_slice = t_aux_slice1 + slice;
748
749 const eT* x_colptr = x.colptr(slice);
750
751 uword i,j;
752 for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
753 {
754 const eT tmp_i = x_colptr[i];
755 const eT tmp_j = x_colptr[j];
756
757 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) -= tmp_i;
758 Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) -= tmp_j;
759 }
760
761 if(i < t_n_cols)
762 {
763 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) -= x_colptr[i];
764 }
765 }
766 }
767 else
768 {
769 if(arma_config::debug)
770 {
771 arma_stop_logic_error( arma_incompat_size_string(t, x, "subtraction") );
772 }
773 }
774 }
775
776
777
778 template<typename eT>
779 template<typename T1>
780 inline
781 void
operator %=(const Base<eT,T1> & in)782 subview_cube<eT>::operator%= (const Base<eT,T1>& in)
783 {
784 arma_extra_debug_sigprint();
785
786 const quasi_unwrap<T1> tmp(in.get_ref());
787
788 const Mat<eT>& x = tmp.M;
789 subview_cube<eT>& t = *this;
790
791 const uword t_n_rows = t.n_rows;
792 const uword t_n_cols = t.n_cols;
793 const uword t_n_slices = t.n_slices;
794
795 const uword x_n_rows = x.n_rows;
796 const uword x_n_cols = x.n_cols;
797
798 if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
799 {
800 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
801
802 const uword t_aux_row1 = t.aux_row1;
803 const uword t_aux_col1 = t.aux_col1;
804 const uword t_aux_slice1 = t.aux_slice1;
805
806 const eT* x_mem = x.memptr();
807
808 uword i,j;
809 for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
810 {
811 const eT tmp_i = x_mem[i];
812 const eT tmp_j = x_mem[j];
813
814 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) *= tmp_i;
815 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) *= tmp_j;
816 }
817
818 if(i < t_n_slices)
819 {
820 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) *= x_mem[i];
821 }
822 }
823 else
824 if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
825 {
826 for(uword col = 0; col < t_n_cols; ++col)
827 {
828 arrayops::inplace_mul( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
829 }
830 }
831 else
832 if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
833 {
834 for(uword i=0; i < t_n_slices; ++i)
835 {
836 arrayops::inplace_mul( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
837 }
838 }
839 else
840 if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
841 {
842 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
843
844 const uword t_aux_row1 = t.aux_row1;
845 const uword t_aux_col1 = t.aux_col1;
846 const uword t_aux_slice1 = t.aux_slice1;
847
848 for(uword slice=0; slice < t_n_slices; ++slice)
849 {
850 const uword mod_slice = t_aux_slice1 + slice;
851
852 const eT* x_colptr = x.colptr(slice);
853
854 uword i,j;
855 for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
856 {
857 const eT tmp_i = x_colptr[i];
858 const eT tmp_j = x_colptr[j];
859
860 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) *= tmp_i;
861 Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) *= tmp_j;
862 }
863
864 if(i < t_n_cols)
865 {
866 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) *= x_colptr[i];
867 }
868 }
869 }
870 else
871 {
872 if(arma_config::debug)
873 {
874 arma_stop_logic_error( arma_incompat_size_string(t, x, "element-wise multiplication") );
875 }
876 }
877 }
878
879
880
881 template<typename eT>
882 template<typename T1>
883 inline
884 void
operator /=(const Base<eT,T1> & in)885 subview_cube<eT>::operator/= (const Base<eT,T1>& in)
886 {
887 arma_extra_debug_sigprint();
888
889 const quasi_unwrap<T1> tmp(in.get_ref());
890
891 const Mat<eT>& x = tmp.M;
892 subview_cube<eT>& t = *this;
893
894 const uword t_n_rows = t.n_rows;
895 const uword t_n_cols = t.n_cols;
896 const uword t_n_slices = t.n_slices;
897
898 const uword x_n_rows = x.n_rows;
899 const uword x_n_cols = x.n_cols;
900
901 if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
902 {
903 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
904
905 const uword t_aux_row1 = t.aux_row1;
906 const uword t_aux_col1 = t.aux_col1;
907 const uword t_aux_slice1 = t.aux_slice1;
908
909 const eT* x_mem = x.memptr();
910
911 uword i,j;
912 for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
913 {
914 const eT tmp_i = x_mem[i];
915 const eT tmp_j = x_mem[j];
916
917 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) /= tmp_i;
918 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) /= tmp_j;
919 }
920
921 if(i < t_n_slices)
922 {
923 Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) /= x_mem[i];
924 }
925 }
926 else
927 if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
928 {
929 for(uword col = 0; col < t_n_cols; ++col)
930 {
931 arrayops::inplace_div( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
932 }
933 }
934 else
935 if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
936 {
937 for(uword i=0; i < t_n_slices; ++i)
938 {
939 arrayops::inplace_div( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
940 }
941 }
942 else
943 if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
944 {
945 Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
946
947 const uword t_aux_row1 = t.aux_row1;
948 const uword t_aux_col1 = t.aux_col1;
949 const uword t_aux_slice1 = t.aux_slice1;
950
951 for(uword slice=0; slice < t_n_slices; ++slice)
952 {
953 const uword mod_slice = t_aux_slice1 + slice;
954
955 const eT* x_colptr = x.colptr(slice);
956
957 uword i,j;
958 for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
959 {
960 const eT tmp_i = x_colptr[i];
961 const eT tmp_j = x_colptr[j];
962
963 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) /= tmp_i;
964 Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) /= tmp_j;
965 }
966
967 if(i < t_n_cols)
968 {
969 Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) /= x_colptr[i];
970 }
971 }
972 }
973 else
974 {
975 if(arma_config::debug)
976 {
977 arma_stop_logic_error( arma_incompat_size_string(t, x, "element-wise division") );
978 }
979 }
980 }
981
982
983
984 template<typename eT>
985 template<typename gen_type>
986 inline
987 void
operator =(const GenCube<eT,gen_type> & in)988 subview_cube<eT>::operator= (const GenCube<eT,gen_type>& in)
989 {
990 arma_extra_debug_sigprint();
991
992 arma_debug_assert_same_size(n_rows, n_cols, n_slices, in.n_rows, in.n_cols, in.n_slices, "copy into subcube");
993
994 in.apply(*this);
995 }
996
997
998
999 //! apply a functor to each element
1000 template<typename eT>
1001 template<typename functor>
1002 inline
1003 void
for_each(functor F)1004 subview_cube<eT>::for_each(functor F)
1005 {
1006 arma_extra_debug_sigprint();
1007
1008 Cube<eT>& Q = const_cast< Cube<eT>& >(m);
1009
1010 const uword start_col = aux_col1;
1011 const uword start_row = aux_row1;
1012 const uword start_slice = aux_slice1;
1013
1014 const uword end_col_plus1 = start_col + n_cols;
1015 const uword end_row_plus1 = start_row + n_rows;
1016 const uword end_slice_plus1 = start_slice + n_slices;
1017
1018 for(uword uslice = start_slice; uslice < end_slice_plus1; ++uslice)
1019 for(uword ucol = start_col; ucol < end_col_plus1; ++ucol )
1020 for(uword urow = start_row; urow < end_row_plus1; ++urow )
1021 {
1022 F( Q.at(urow, ucol, uslice) );
1023 }
1024 }
1025
1026
1027
1028 template<typename eT>
1029 template<typename functor>
1030 inline
1031 void
for_each(functor F) const1032 subview_cube<eT>::for_each(functor F) const
1033 {
1034 arma_extra_debug_sigprint();
1035
1036 const Cube<eT>& Q = m;
1037
1038 const uword start_col = aux_col1;
1039 const uword start_row = aux_row1;
1040 const uword start_slice = aux_slice1;
1041
1042 const uword end_col_plus1 = start_col + n_cols;
1043 const uword end_row_plus1 = start_row + n_rows;
1044 const uword end_slice_plus1 = start_slice + n_slices;
1045
1046 for(uword uslice = start_slice; uslice < end_slice_plus1; ++uslice)
1047 for(uword ucol = start_col; ucol < end_col_plus1; ++ucol )
1048 for(uword urow = start_row; urow < end_row_plus1; ++urow )
1049 {
1050 F( Q.at(urow, ucol, uslice) );
1051 }
1052 }
1053
1054
1055
1056 //! transform each element in the subview using a functor
1057 template<typename eT>
1058 template<typename functor>
1059 inline
1060 void
transform(functor F)1061 subview_cube<eT>::transform(functor F)
1062 {
1063 arma_extra_debug_sigprint();
1064
1065 Cube<eT>& Q = const_cast< Cube<eT>& >(m);
1066
1067 const uword start_col = aux_col1;
1068 const uword start_row = aux_row1;
1069 const uword start_slice = aux_slice1;
1070
1071 const uword end_col_plus1 = start_col + n_cols;
1072 const uword end_row_plus1 = start_row + n_rows;
1073 const uword end_slice_plus1 = start_slice + n_slices;
1074
1075 for(uword uslice = start_slice; uslice < end_slice_plus1; ++uslice)
1076 for(uword ucol = start_col; ucol < end_col_plus1; ++ucol )
1077 for(uword urow = start_row; urow < end_row_plus1; ++urow )
1078 {
1079 Q.at(urow, ucol, uslice) = eT( F( Q.at(urow, ucol, uslice) ) );
1080 }
1081 }
1082
1083
1084
1085 //! imbue (fill) the subview with values provided by a functor
1086 template<typename eT>
1087 template<typename functor>
1088 inline
1089 void
imbue(functor F)1090 subview_cube<eT>::imbue(functor F)
1091 {
1092 arma_extra_debug_sigprint();
1093
1094 Cube<eT>& Q = const_cast< Cube<eT>& >(m);
1095
1096 const uword start_col = aux_col1;
1097 const uword start_row = aux_row1;
1098 const uword start_slice = aux_slice1;
1099
1100 const uword end_col_plus1 = start_col + n_cols;
1101 const uword end_row_plus1 = start_row + n_rows;
1102 const uword end_slice_plus1 = start_slice + n_slices;
1103
1104 for(uword uslice = start_slice; uslice < end_slice_plus1; ++uslice)
1105 for(uword ucol = start_col; ucol < end_col_plus1; ++ucol )
1106 for(uword urow = start_row; urow < end_row_plus1; ++urow )
1107 {
1108 Q.at(urow, ucol, uslice) = eT( F() );
1109 }
1110 }
1111
1112
1113
1114 //! apply a lambda function to each slice, where each slice is interpreted as a matrix
1115 template<typename eT>
1116 inline
1117 void
each_slice(const std::function<void (Mat<eT> &)> & F)1118 subview_cube<eT>::each_slice(const std::function< void(Mat<eT>&) >& F)
1119 {
1120 arma_extra_debug_sigprint();
1121
1122 Mat<eT> tmp1(n_rows, n_cols, arma_nozeros_indicator());
1123 Mat<eT> tmp2('j', tmp1.memptr(), n_rows, n_cols);
1124
1125 for(uword slice_id=0; slice_id < n_slices; ++slice_id)
1126 {
1127 for(uword col_id=0; col_id < n_cols; ++col_id)
1128 {
1129 arrayops::copy( tmp1.colptr(col_id), slice_colptr(slice_id, col_id), n_rows );
1130 }
1131
1132 F(tmp2);
1133
1134 for(uword col_id=0; col_id < n_cols; ++col_id)
1135 {
1136 arrayops::copy( slice_colptr(slice_id, col_id), tmp1.colptr(col_id), n_rows );
1137 }
1138 }
1139 }
1140
1141
1142
1143 template<typename eT>
1144 inline
1145 void
each_slice(const std::function<void (const Mat<eT> &)> & F) const1146 subview_cube<eT>::each_slice(const std::function< void(const Mat<eT>&) >& F) const
1147 {
1148 arma_extra_debug_sigprint();
1149
1150 Mat<eT> tmp1(n_rows, n_cols, arma_nozeros_indicator());
1151 const Mat<eT> tmp2('j', tmp1.memptr(), n_rows, n_cols);
1152
1153 for(uword slice_id=0; slice_id < n_slices; ++slice_id)
1154 {
1155 for(uword col_id=0; col_id < n_cols; ++col_id)
1156 {
1157 arrayops::copy( tmp1.colptr(col_id), slice_colptr(slice_id, col_id), n_rows );
1158 }
1159
1160 F(tmp2);
1161 }
1162 }
1163
1164
1165
1166 template<typename eT>
1167 inline
1168 void
replace(const eT old_val,const eT new_val)1169 subview_cube<eT>::replace(const eT old_val, const eT new_val)
1170 {
1171 arma_extra_debug_sigprint();
1172
1173 const uword local_n_rows = n_rows;
1174 const uword local_n_cols = n_cols;
1175 const uword local_n_slices = n_slices;
1176
1177 for(uword slice = 0; slice < local_n_slices; ++slice)
1178 {
1179 for(uword col = 0; col < local_n_cols; ++col)
1180 {
1181 arrayops::replace(slice_colptr(slice,col), local_n_rows, old_val, new_val);
1182 }
1183 }
1184 }
1185
1186
1187
1188 template<typename eT>
1189 inline
1190 void
clean(const typename get_pod_type<eT>::result threshold)1191 subview_cube<eT>::clean(const typename get_pod_type<eT>::result threshold)
1192 {
1193 arma_extra_debug_sigprint();
1194
1195 const uword local_n_rows = n_rows;
1196 const uword local_n_cols = n_cols;
1197 const uword local_n_slices = n_slices;
1198
1199 for(uword slice = 0; slice < local_n_slices; ++slice)
1200 {
1201 for(uword col = 0; col < local_n_cols; ++col)
1202 {
1203 arrayops::clean( slice_colptr(slice,col), local_n_rows, threshold );
1204 }
1205 }
1206 }
1207
1208
1209
1210 template<typename eT>
1211 inline
1212 void
clamp(const eT min_val,const eT max_val)1213 subview_cube<eT>::clamp(const eT min_val, const eT max_val)
1214 {
1215 arma_extra_debug_sigprint();
1216
1217 if(is_cx<eT>::no)
1218 {
1219 arma_debug_check( (access::tmp_real(min_val) > access::tmp_real(max_val)), "subview_cube::clamp(): min_val must be less than max_val" );
1220 }
1221 else
1222 {
1223 arma_debug_check( (access::tmp_real(min_val) > access::tmp_real(max_val)), "subview_cube::clamp(): real(min_val) must be less than real(max_val)" );
1224 arma_debug_check( (access::tmp_imag(min_val) > access::tmp_imag(max_val)), "subview_cube::clamp(): imag(min_val) must be less than imag(max_val)" );
1225 }
1226
1227 const uword local_n_rows = n_rows;
1228 const uword local_n_cols = n_cols;
1229 const uword local_n_slices = n_slices;
1230
1231 for(uword slice = 0; slice < local_n_slices; ++slice)
1232 {
1233 for(uword col = 0; col < local_n_cols; ++col)
1234 {
1235 arrayops::clamp( slice_colptr(slice,col), local_n_rows, min_val, max_val );
1236 }
1237 }
1238 }
1239
1240
1241
1242 template<typename eT>
1243 inline
1244 void
fill(const eT val)1245 subview_cube<eT>::fill(const eT val)
1246 {
1247 arma_extra_debug_sigprint();
1248
1249 const uword local_n_rows = n_rows;
1250 const uword local_n_cols = n_cols;
1251 const uword local_n_slices = n_slices;
1252
1253 for(uword slice = 0; slice < local_n_slices; ++slice)
1254 {
1255 for(uword col = 0; col < local_n_cols; ++col)
1256 {
1257 arrayops::inplace_set( slice_colptr(slice,col), val, local_n_rows );
1258 }
1259 }
1260 }
1261
1262
1263
1264 template<typename eT>
1265 inline
1266 void
zeros()1267 subview_cube<eT>::zeros()
1268 {
1269 arma_extra_debug_sigprint();
1270
1271 const uword local_n_rows = n_rows;
1272 const uword local_n_cols = n_cols;
1273 const uword local_n_slices = n_slices;
1274
1275 for(uword slice = 0; slice < local_n_slices; ++slice)
1276 {
1277 for(uword col = 0; col < local_n_cols; ++col)
1278 {
1279 arrayops::fill_zeros( slice_colptr(slice,col), local_n_rows );
1280 }
1281 }
1282 }
1283
1284
1285
1286 template<typename eT>
1287 inline
1288 void
ones()1289 subview_cube<eT>::ones()
1290 {
1291 arma_extra_debug_sigprint();
1292
1293 fill(eT(1));
1294 }
1295
1296
1297
1298 template<typename eT>
1299 inline
1300 void
randu()1301 subview_cube<eT>::randu()
1302 {
1303 arma_extra_debug_sigprint();
1304
1305 const uword local_n_rows = n_rows;
1306 const uword local_n_cols = n_cols;
1307 const uword local_n_slices = n_slices;
1308
1309 for(uword slice = 0; slice < local_n_slices; ++slice)
1310 {
1311 for(uword col = 0; col < local_n_cols; ++col)
1312 {
1313 arma_rng::randu<eT>::fill( slice_colptr(slice,col), local_n_rows );
1314 }
1315 }
1316 }
1317
1318
1319
1320 template<typename eT>
1321 inline
1322 void
randn()1323 subview_cube<eT>::randn()
1324 {
1325 arma_extra_debug_sigprint();
1326
1327 const uword local_n_rows = n_rows;
1328 const uword local_n_cols = n_cols;
1329 const uword local_n_slices = n_slices;
1330
1331 for(uword slice = 0; slice < local_n_slices; ++slice)
1332 {
1333 for(uword col = 0; col < local_n_cols; ++col)
1334 {
1335 arma_rng::randn<eT>::fill( slice_colptr(slice,col), local_n_rows );
1336 }
1337 }
1338 }
1339
1340
1341
1342 template<typename eT>
1343 inline
1344 arma_warn_unused
1345 bool
is_finite() const1346 subview_cube<eT>::is_finite() const
1347 {
1348 arma_extra_debug_sigprint();
1349
1350 const uword local_n_rows = n_rows;
1351 const uword local_n_cols = n_cols;
1352 const uword local_n_slices = n_slices;
1353
1354 for(uword slice = 0; slice < local_n_slices; ++slice)
1355 {
1356 for(uword col = 0; col < local_n_cols; ++col)
1357 {
1358 if(arrayops::is_finite(slice_colptr(slice,col), local_n_rows) == false) { return false; }
1359 }
1360 }
1361
1362 return true;
1363 }
1364
1365
1366
1367 template<typename eT>
1368 inline
1369 arma_warn_unused
1370 bool
is_zero(const typename get_pod_type<eT>::result tol) const1371 subview_cube<eT>::is_zero(const typename get_pod_type<eT>::result tol) const
1372 {
1373 arma_extra_debug_sigprint();
1374
1375 const uword local_n_rows = n_rows;
1376 const uword local_n_cols = n_cols;
1377 const uword local_n_slices = n_slices;
1378
1379 for(uword slice = 0; slice < local_n_slices; ++slice)
1380 {
1381 for(uword col = 0; col < local_n_cols; ++col)
1382 {
1383 if(arrayops::is_zero(slice_colptr(slice,col), local_n_rows, tol) == false) { return false; }
1384 }
1385 }
1386
1387 return true;
1388 }
1389
1390
1391
1392 template<typename eT>
1393 inline
1394 arma_warn_unused
1395 bool
has_inf() const1396 subview_cube<eT>::has_inf() const
1397 {
1398 arma_extra_debug_sigprint();
1399
1400 const uword local_n_rows = n_rows;
1401 const uword local_n_cols = n_cols;
1402 const uword local_n_slices = n_slices;
1403
1404 for(uword slice = 0; slice < local_n_slices; ++slice)
1405 {
1406 for(uword col = 0; col < local_n_cols; ++col)
1407 {
1408 if(arrayops::has_inf(slice_colptr(slice,col), local_n_rows)) { return true; }
1409 }
1410 }
1411
1412 return false;
1413 }
1414
1415
1416
1417 template<typename eT>
1418 inline
1419 arma_warn_unused
1420 bool
has_nan() const1421 subview_cube<eT>::has_nan() const
1422 {
1423 arma_extra_debug_sigprint();
1424
1425 const uword local_n_rows = n_rows;
1426 const uword local_n_cols = n_cols;
1427 const uword local_n_slices = n_slices;
1428
1429 for(uword slice = 0; slice < local_n_slices; ++slice)
1430 {
1431 for(uword col = 0; col < local_n_cols; ++col)
1432 {
1433 if(arrayops::has_nan(slice_colptr(slice,col), local_n_rows)) { return true; }
1434 }
1435 }
1436
1437 return false;
1438 }
1439
1440
1441
1442 template<typename eT>
1443 inline
1444 eT
at_alt(const uword i) const1445 subview_cube<eT>::at_alt(const uword i) const
1446 {
1447 return operator[](i);
1448 }
1449
1450
1451
1452 template<typename eT>
1453 inline
1454 eT&
operator [](const uword i)1455 subview_cube<eT>::operator[](const uword i)
1456 {
1457 const uword in_slice = i / n_elem_slice;
1458 const uword offset = in_slice * n_elem_slice;
1459 const uword j = i - offset;
1460
1461 const uword in_col = j / n_rows;
1462 const uword in_row = j % n_rows;
1463
1464 const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
1465
1466 return access::rw( (const_cast< Cube<eT>& >(m)).mem[index] );
1467 }
1468
1469
1470
1471 template<typename eT>
1472 inline
1473 eT
operator [](const uword i) const1474 subview_cube<eT>::operator[](const uword i) const
1475 {
1476 const uword in_slice = i / n_elem_slice;
1477 const uword offset = in_slice * n_elem_slice;
1478 const uword j = i - offset;
1479
1480 const uword in_col = j / n_rows;
1481 const uword in_row = j % n_rows;
1482
1483 const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
1484
1485 return m.mem[index];
1486 }
1487
1488
1489
1490 template<typename eT>
1491 inline
1492 eT&
operator ()(const uword i)1493 subview_cube<eT>::operator()(const uword i)
1494 {
1495 arma_debug_check_bounds( (i >= n_elem), "subview_cube::operator(): index out of bounds" );
1496
1497 const uword in_slice = i / n_elem_slice;
1498 const uword offset = in_slice * n_elem_slice;
1499 const uword j = i - offset;
1500
1501 const uword in_col = j / n_rows;
1502 const uword in_row = j % n_rows;
1503
1504 const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
1505
1506 return access::rw( (const_cast< Cube<eT>& >(m)).mem[index] );
1507 }
1508
1509
1510
1511 template<typename eT>
1512 inline
1513 eT
operator ()(const uword i) const1514 subview_cube<eT>::operator()(const uword i) const
1515 {
1516 arma_debug_check_bounds( (i >= n_elem), "subview_cube::operator(): index out of bounds" );
1517
1518 const uword in_slice = i / n_elem_slice;
1519 const uword offset = in_slice * n_elem_slice;
1520 const uword j = i - offset;
1521
1522 const uword in_col = j / n_rows;
1523 const uword in_row = j % n_rows;
1524
1525 const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
1526
1527 return m.mem[index];
1528 }
1529
1530
1531
1532 template<typename eT>
1533 arma_inline
1534 eT&
operator ()(const uword in_row,const uword in_col,const uword in_slice)1535 subview_cube<eT>::operator()(const uword in_row, const uword in_col, const uword in_slice)
1536 {
1537 arma_debug_check_bounds( ( (in_row >= n_rows) || (in_col >= n_cols) || (in_slice >= n_slices) ), "subview_cube::operator(): location out of bounds" );
1538
1539 const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
1540
1541 return access::rw( (const_cast< Cube<eT>& >(m)).mem[index] );
1542 }
1543
1544
1545
1546 template<typename eT>
1547 arma_inline
1548 eT
operator ()(const uword in_row,const uword in_col,const uword in_slice) const1549 subview_cube<eT>::operator()(const uword in_row, const uword in_col, const uword in_slice) const
1550 {
1551 arma_debug_check_bounds( ( (in_row >= n_rows) || (in_col >= n_cols) || (in_slice >= n_slices) ), "subview_cube::operator(): location out of bounds" );
1552
1553 const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
1554
1555 return m.mem[index];
1556 }
1557
1558
1559
1560 template<typename eT>
1561 arma_inline
1562 eT&
at(const uword in_row,const uword in_col,const uword in_slice)1563 subview_cube<eT>::at(const uword in_row, const uword in_col, const uword in_slice)
1564 {
1565 const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
1566
1567 return access::rw( (const_cast< Cube<eT>& >(m)).mem[index] );
1568 }
1569
1570
1571
1572 template<typename eT>
1573 arma_inline
1574 eT
at(const uword in_row,const uword in_col,const uword in_slice) const1575 subview_cube<eT>::at(const uword in_row, const uword in_col, const uword in_slice) const
1576 {
1577 const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
1578
1579 return m.mem[index];
1580 }
1581
1582
1583
1584 template<typename eT>
1585 arma_inline
1586 eT*
slice_colptr(const uword in_slice,const uword in_col)1587 subview_cube<eT>::slice_colptr(const uword in_slice, const uword in_col)
1588 {
1589 return & access::rw((const_cast< Cube<eT>& >(m)).mem[ (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 ]);
1590 }
1591
1592
1593
1594 template<typename eT>
1595 arma_inline
1596 const eT*
slice_colptr(const uword in_slice,const uword in_col) const1597 subview_cube<eT>::slice_colptr(const uword in_slice, const uword in_col) const
1598 {
1599 return & m.mem[ (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 ];
1600 }
1601
1602
1603
1604 template<typename eT>
1605 template<typename eT2>
1606 inline
1607 bool
check_overlap(const subview_cube<eT2> & x) const1608 subview_cube<eT>::check_overlap(const subview_cube<eT2>& x) const
1609 {
1610 if(is_same_type<eT,eT2>::value == false) { return false; }
1611
1612 const subview_cube<eT>& t = (*this);
1613
1614 if(void_ptr(&(t.m)) != void_ptr(&(x.m))) { return false; }
1615
1616 if( (t.n_elem == 0) || (x.n_elem == 0) ) { return false; }
1617
1618 const uword t_row_start = t.aux_row1;
1619 const uword t_row_end_p1 = t_row_start + t.n_rows;
1620
1621 const uword t_col_start = t.aux_col1;
1622 const uword t_col_end_p1 = t_col_start + t.n_cols;
1623
1624 const uword t_slice_start = t.aux_slice1;
1625 const uword t_slice_end_p1 = t_slice_start + t.n_slices;
1626
1627
1628 const uword x_row_start = x.aux_row1;
1629 const uword x_row_end_p1 = x_row_start + x.n_rows;
1630
1631 const uword x_col_start = x.aux_col1;
1632 const uword x_col_end_p1 = x_col_start + x.n_cols;
1633
1634 const uword x_slice_start = x.aux_slice1;
1635 const uword x_slice_end_p1 = x_slice_start + x.n_slices;
1636
1637
1638 const bool outside_rows = ( (x_row_start >= t_row_end_p1 ) || (t_row_start >= x_row_end_p1 ) );
1639 const bool outside_cols = ( (x_col_start >= t_col_end_p1 ) || (t_col_start >= x_col_end_p1 ) );
1640 const bool outside_slices = ( (x_slice_start >= t_slice_end_p1) || (t_slice_start >= x_slice_end_p1) );
1641
1642 return ( (outside_rows == false) && (outside_cols == false) && (outside_slices == false) );
1643 }
1644
1645
1646
1647 template<typename eT>
1648 inline
1649 bool
check_overlap(const Mat<eT> & x) const1650 subview_cube<eT>::check_overlap(const Mat<eT>& x) const
1651 {
1652 const subview_cube<eT>& t = *this;
1653
1654 const uword t_aux_slice1 = t.aux_slice1;
1655 const uword t_aux_slice2_plus_1 = t_aux_slice1 + t.n_slices;
1656
1657 for(uword slice = t_aux_slice1; slice < t_aux_slice2_plus_1; ++slice)
1658 {
1659 if(t.m.mat_ptrs[slice] != nullptr)
1660 {
1661 const Mat<eT>& y = *(t.m.mat_ptrs[slice]);
1662
1663 if( x.memptr() == y.memptr() ) { return true; }
1664 }
1665 }
1666
1667 return false;
1668 }
1669
1670
1671
1672 //! cube X = Y.subcube(...)
1673 template<typename eT>
1674 inline
1675 void
extract(Cube<eT> & out,const subview_cube<eT> & in)1676 subview_cube<eT>::extract(Cube<eT>& out, const subview_cube<eT>& in)
1677 {
1678 arma_extra_debug_sigprint();
1679
1680 // NOTE: we're assuming that the cube has already been set to the correct size and there is no aliasing;
1681 // size setting and alias checking is done by either the Cube contructor or operator=()
1682
1683 const uword n_rows = in.n_rows;
1684 const uword n_cols = in.n_cols;
1685 const uword n_slices = in.n_slices;
1686
1687 arma_extra_debug_print(arma_str::format("out.n_rows = %d out.n_cols = %d out.n_slices = %d in.m.n_rows = %d in.m.n_cols = %d in.m.n_slices = %d") % out.n_rows % out.n_cols % out.n_slices % in.m.n_rows % in.m.n_cols % in.m.n_slices);
1688
1689 if( (in.aux_row1 == 0) && (n_rows == in.m.n_rows) )
1690 {
1691 for(uword s=0; s < n_slices; ++s)
1692 {
1693 arrayops::copy( out.slice_colptr(s,0), in.slice_colptr(s,0), in.n_elem_slice );
1694 }
1695
1696 return;
1697 }
1698
1699 for(uword s=0; s < n_slices; ++s)
1700 for(uword c=0; c < n_cols; ++c)
1701 {
1702 arrayops::copy( out.slice_colptr(s,c), in.slice_colptr(s,c), n_rows );
1703 }
1704 }
1705
1706
1707
1708 //! cube X += Y.subcube(...)
1709 template<typename eT>
1710 inline
1711 void
plus_inplace(Cube<eT> & out,const subview_cube<eT> & in)1712 subview_cube<eT>::plus_inplace(Cube<eT>& out, const subview_cube<eT>& in)
1713 {
1714 arma_extra_debug_sigprint();
1715
1716 arma_debug_assert_same_size(out, in, "addition");
1717
1718 const uword n_rows = out.n_rows;
1719 const uword n_cols = out.n_cols;
1720 const uword n_slices = out.n_slices;
1721
1722 for(uword slice = 0; slice<n_slices; ++slice)
1723 {
1724 for(uword col = 0; col<n_cols; ++col)
1725 {
1726 arrayops::inplace_plus( out.slice_colptr(slice,col), in.slice_colptr(slice,col), n_rows );
1727 }
1728 }
1729 }
1730
1731
1732
1733 //! cube X -= Y.subcube(...)
1734 template<typename eT>
1735 inline
1736 void
minus_inplace(Cube<eT> & out,const subview_cube<eT> & in)1737 subview_cube<eT>::minus_inplace(Cube<eT>& out, const subview_cube<eT>& in)
1738 {
1739 arma_extra_debug_sigprint();
1740
1741 arma_debug_assert_same_size(out, in, "subtraction");
1742
1743 const uword n_rows = out.n_rows;
1744 const uword n_cols = out.n_cols;
1745 const uword n_slices = out.n_slices;
1746
1747 for(uword slice = 0; slice<n_slices; ++slice)
1748 {
1749 for(uword col = 0; col<n_cols; ++col)
1750 {
1751 arrayops::inplace_minus( out.slice_colptr(slice,col), in.slice_colptr(slice,col), n_rows );
1752 }
1753 }
1754 }
1755
1756
1757
1758 //! cube X %= Y.subcube(...)
1759 template<typename eT>
1760 inline
1761 void
schur_inplace(Cube<eT> & out,const subview_cube<eT> & in)1762 subview_cube<eT>::schur_inplace(Cube<eT>& out, const subview_cube<eT>& in)
1763 {
1764 arma_extra_debug_sigprint();
1765
1766 arma_debug_assert_same_size(out, in, "element-wise multiplication");
1767
1768 const uword n_rows = out.n_rows;
1769 const uword n_cols = out.n_cols;
1770 const uword n_slices = out.n_slices;
1771
1772 for(uword slice = 0; slice<n_slices; ++slice)
1773 {
1774 for(uword col = 0; col<n_cols; ++col)
1775 {
1776 arrayops::inplace_mul( out.slice_colptr(slice,col), in.slice_colptr(slice,col), n_rows );
1777 }
1778 }
1779 }
1780
1781
1782
1783 //! cube X /= Y.subcube(...)
1784 template<typename eT>
1785 inline
1786 void
div_inplace(Cube<eT> & out,const subview_cube<eT> & in)1787 subview_cube<eT>::div_inplace(Cube<eT>& out, const subview_cube<eT>& in)
1788 {
1789 arma_extra_debug_sigprint();
1790
1791 arma_debug_assert_same_size(out, in, "element-wise division");
1792
1793 const uword n_rows = out.n_rows;
1794 const uword n_cols = out.n_cols;
1795 const uword n_slices = out.n_slices;
1796
1797 for(uword slice = 0; slice<n_slices; ++slice)
1798 {
1799 for(uword col = 0; col<n_cols; ++col)
1800 {
1801 arrayops::inplace_div( out.slice_colptr(slice,col), in.slice_colptr(slice,col), n_rows );
1802 }
1803 }
1804 }
1805
1806
1807
1808 //! mat X = Y.subcube(...)
1809 template<typename eT>
1810 inline
1811 void
extract(Mat<eT> & out,const subview_cube<eT> & in)1812 subview_cube<eT>::extract(Mat<eT>& out, const subview_cube<eT>& in)
1813 {
1814 arma_extra_debug_sigprint();
1815
1816 arma_debug_assert_cube_as_mat(out, in, "copy into matrix", false);
1817
1818 const uword in_n_rows = in.n_rows;
1819 const uword in_n_cols = in.n_cols;
1820 const uword in_n_slices = in.n_slices;
1821
1822 const uword out_vec_state = out.vec_state;
1823
1824 if(in_n_slices == 1)
1825 {
1826 out.set_size(in_n_rows, in_n_cols);
1827
1828 for(uword col=0; col < in_n_cols; ++col)
1829 {
1830 arrayops::copy( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
1831 }
1832 }
1833 else
1834 {
1835 if(out_vec_state == 0)
1836 {
1837 if(in_n_cols == 1)
1838 {
1839 out.set_size(in_n_rows, in_n_slices);
1840
1841 for(uword i=0; i < in_n_slices; ++i)
1842 {
1843 arrayops::copy( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
1844 }
1845 }
1846 else
1847 if(in_n_rows == 1)
1848 {
1849 const Cube<eT>& Q = in.m;
1850
1851 const uword in_aux_row1 = in.aux_row1;
1852 const uword in_aux_col1 = in.aux_col1;
1853 const uword in_aux_slice1 = in.aux_slice1;
1854
1855 out.set_size(in_n_cols, in_n_slices);
1856
1857 for(uword slice=0; slice < in_n_slices; ++slice)
1858 {
1859 const uword mod_slice = in_aux_slice1 + slice;
1860
1861 eT* out_colptr = out.colptr(slice);
1862
1863 uword i,j;
1864 for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
1865 {
1866 const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
1867 const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
1868
1869 out_colptr[i] = tmp_i;
1870 out_colptr[j] = tmp_j;
1871 }
1872
1873 if(i < in_n_cols)
1874 {
1875 out_colptr[i] = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
1876 }
1877 }
1878 }
1879 }
1880 else
1881 {
1882 out.set_size(in_n_slices);
1883
1884 eT* out_mem = out.memptr();
1885
1886 const Cube<eT>& Q = in.m;
1887
1888 const uword in_aux_row1 = in.aux_row1;
1889 const uword in_aux_col1 = in.aux_col1;
1890 const uword in_aux_slice1 = in.aux_slice1;
1891
1892 for(uword i=0; i<in_n_slices; ++i)
1893 {
1894 out_mem[i] = Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
1895 }
1896 }
1897 }
1898 }
1899
1900
1901
1902 //! mat X += Y.subcube(...)
1903 template<typename eT>
1904 inline
1905 void
plus_inplace(Mat<eT> & out,const subview_cube<eT> & in)1906 subview_cube<eT>::plus_inplace(Mat<eT>& out, const subview_cube<eT>& in)
1907 {
1908 arma_extra_debug_sigprint();
1909
1910 arma_debug_assert_cube_as_mat(out, in, "addition", true);
1911
1912 const uword in_n_rows = in.n_rows;
1913 const uword in_n_cols = in.n_cols;
1914 const uword in_n_slices = in.n_slices;
1915
1916 const uword out_n_rows = out.n_rows;
1917 const uword out_n_cols = out.n_cols;
1918 const uword out_vec_state = out.vec_state;
1919
1920 if(in_n_slices == 1)
1921 {
1922 if( (arma_config::debug) && ((out_n_rows != in_n_rows) || (out_n_cols != in_n_cols)) )
1923 {
1924 std::ostringstream tmp;
1925
1926 tmp
1927 << "in-place addition: "
1928 << out_n_rows << 'x' << out_n_cols << " output matrix is incompatible with "
1929 << in_n_rows << 'x' << in_n_cols << 'x' << in_n_slices << " cube interpreted as "
1930 << in_n_rows << 'x' << in_n_cols << " matrix";
1931
1932 arma_stop_logic_error(tmp.str());
1933 }
1934
1935 for(uword col=0; col < in_n_cols; ++col)
1936 {
1937 arrayops::inplace_plus( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
1938 }
1939 }
1940 else
1941 {
1942 if(out_vec_state == 0)
1943 {
1944 if( (in_n_rows == out_n_rows) && (in_n_cols == 1) && (in_n_slices == out_n_cols) )
1945 {
1946 for(uword i=0; i < in_n_slices; ++i)
1947 {
1948 arrayops::inplace_plus( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
1949 }
1950 }
1951 else
1952 if( (in_n_rows == 1) && (in_n_cols == out_n_rows) && (in_n_slices == out_n_cols) )
1953 {
1954 const Cube<eT>& Q = in.m;
1955
1956 const uword in_aux_row1 = in.aux_row1;
1957 const uword in_aux_col1 = in.aux_col1;
1958 const uword in_aux_slice1 = in.aux_slice1;
1959
1960 for(uword slice=0; slice < in_n_slices; ++slice)
1961 {
1962 const uword mod_slice = in_aux_slice1 + slice;
1963
1964 eT* out_colptr = out.colptr(slice);
1965
1966 uword i,j;
1967 for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
1968 {
1969 const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
1970 const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
1971
1972 out_colptr[i] += tmp_i;
1973 out_colptr[j] += tmp_j;
1974 }
1975
1976 if(i < in_n_cols)
1977 {
1978 out_colptr[i] += Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
1979 }
1980 }
1981 }
1982 }
1983 else
1984 {
1985 eT* out_mem = out.memptr();
1986
1987 const Cube<eT>& Q = in.m;
1988
1989 const uword in_aux_row1 = in.aux_row1;
1990 const uword in_aux_col1 = in.aux_col1;
1991 const uword in_aux_slice1 = in.aux_slice1;
1992
1993 for(uword i=0; i<in_n_slices; ++i)
1994 {
1995 out_mem[i] += Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
1996 }
1997 }
1998 }
1999 }
2000
2001
2002
2003 //! mat X -= Y.subcube(...)
2004 template<typename eT>
2005 inline
2006 void
minus_inplace(Mat<eT> & out,const subview_cube<eT> & in)2007 subview_cube<eT>::minus_inplace(Mat<eT>& out, const subview_cube<eT>& in)
2008 {
2009 arma_extra_debug_sigprint();
2010
2011 arma_debug_assert_cube_as_mat(out, in, "subtraction", true);
2012
2013 const uword in_n_rows = in.n_rows;
2014 const uword in_n_cols = in.n_cols;
2015 const uword in_n_slices = in.n_slices;
2016
2017 const uword out_n_rows = out.n_rows;
2018 const uword out_n_cols = out.n_cols;
2019 const uword out_vec_state = out.vec_state;
2020
2021 if(in_n_slices == 1)
2022 {
2023 if( (arma_config::debug) && ((out_n_rows != in_n_rows) || (out_n_cols != in_n_cols)) )
2024 {
2025 std::ostringstream tmp;
2026
2027 tmp
2028 << "in-place subtraction: "
2029 << out_n_rows << 'x' << out_n_cols << " output matrix is incompatible with "
2030 << in_n_rows << 'x' << in_n_cols << 'x' << in_n_slices << " cube interpreted as "
2031 << in_n_rows << 'x' << in_n_cols << " matrix";
2032
2033 arma_stop_logic_error(tmp.str());
2034 }
2035
2036 for(uword col=0; col < in_n_cols; ++col)
2037 {
2038 arrayops::inplace_minus( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
2039 }
2040 }
2041 else
2042 {
2043 if(out_vec_state == 0)
2044 {
2045 if( (in_n_rows == out_n_rows) && (in_n_cols == 1) && (in_n_slices == out_n_cols) )
2046 {
2047 for(uword i=0; i < in_n_slices; ++i)
2048 {
2049 arrayops::inplace_minus( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
2050 }
2051 }
2052 else
2053 if( (in_n_rows == 1) && (in_n_cols == out_n_rows) && (in_n_slices == out_n_cols) )
2054 {
2055 const Cube<eT>& Q = in.m;
2056
2057 const uword in_aux_row1 = in.aux_row1;
2058 const uword in_aux_col1 = in.aux_col1;
2059 const uword in_aux_slice1 = in.aux_slice1;
2060
2061 for(uword slice=0; slice < in_n_slices; ++slice)
2062 {
2063 const uword mod_slice = in_aux_slice1 + slice;
2064
2065 eT* out_colptr = out.colptr(slice);
2066
2067 uword i,j;
2068 for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
2069 {
2070 const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
2071 const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
2072
2073 out_colptr[i] -= tmp_i;
2074 out_colptr[j] -= tmp_j;
2075 }
2076
2077 if(i < in_n_cols)
2078 {
2079 out_colptr[i] -= Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
2080 }
2081 }
2082 }
2083 }
2084 else
2085 {
2086 eT* out_mem = out.memptr();
2087
2088 const Cube<eT>& Q = in.m;
2089
2090 const uword in_aux_row1 = in.aux_row1;
2091 const uword in_aux_col1 = in.aux_col1;
2092 const uword in_aux_slice1 = in.aux_slice1;
2093
2094 for(uword i=0; i<in_n_slices; ++i)
2095 {
2096 out_mem[i] -= Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
2097 }
2098 }
2099 }
2100 }
2101
2102
2103
2104 //! mat X %= Y.subcube(...)
2105 template<typename eT>
2106 inline
2107 void
schur_inplace(Mat<eT> & out,const subview_cube<eT> & in)2108 subview_cube<eT>::schur_inplace(Mat<eT>& out, const subview_cube<eT>& in)
2109 {
2110 arma_extra_debug_sigprint();
2111
2112 arma_debug_assert_cube_as_mat(out, in, "element-wise multiplication", true);
2113
2114 const uword in_n_rows = in.n_rows;
2115 const uword in_n_cols = in.n_cols;
2116 const uword in_n_slices = in.n_slices;
2117
2118 const uword out_n_rows = out.n_rows;
2119 const uword out_n_cols = out.n_cols;
2120 const uword out_vec_state = out.vec_state;
2121
2122 if(in_n_slices == 1)
2123 {
2124 if( (arma_config::debug) && ((out_n_rows != in_n_rows) || (out_n_cols != in_n_cols)) )
2125 {
2126 std::ostringstream tmp;
2127
2128 tmp
2129 << "in-place element-wise multiplication: "
2130 << out_n_rows << 'x' << out_n_cols << " output matrix is incompatible with "
2131 << in_n_rows << 'x' << in_n_cols << 'x' << in_n_slices << " cube interpreted as "
2132 << in_n_rows << 'x' << in_n_cols << " matrix";
2133
2134 arma_stop_logic_error(tmp.str());
2135 }
2136
2137 for(uword col=0; col < in_n_cols; ++col)
2138 {
2139 arrayops::inplace_mul( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
2140 }
2141 }
2142 else
2143 {
2144 if(out_vec_state == 0)
2145 {
2146 if( (in_n_rows == out_n_rows) && (in_n_cols == 1) && (in_n_slices == out_n_cols) )
2147 {
2148 for(uword i=0; i < in_n_slices; ++i)
2149 {
2150 arrayops::inplace_mul( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
2151 }
2152 }
2153 else
2154 if( (in_n_rows == 1) && (in_n_cols == out_n_rows) && (in_n_slices == out_n_cols) )
2155 {
2156 const Cube<eT>& Q = in.m;
2157
2158 const uword in_aux_row1 = in.aux_row1;
2159 const uword in_aux_col1 = in.aux_col1;
2160 const uword in_aux_slice1 = in.aux_slice1;
2161
2162 for(uword slice=0; slice < in_n_slices; ++slice)
2163 {
2164 const uword mod_slice = in_aux_slice1 + slice;
2165
2166 eT* out_colptr = out.colptr(slice);
2167
2168 uword i,j;
2169 for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
2170 {
2171 const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
2172 const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
2173
2174 out_colptr[i] *= tmp_i;
2175 out_colptr[j] *= tmp_j;
2176 }
2177
2178 if(i < in_n_cols)
2179 {
2180 out_colptr[i] *= Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
2181 }
2182 }
2183 }
2184 }
2185 else
2186 {
2187 eT* out_mem = out.memptr();
2188
2189 const Cube<eT>& Q = in.m;
2190
2191 const uword in_aux_row1 = in.aux_row1;
2192 const uword in_aux_col1 = in.aux_col1;
2193 const uword in_aux_slice1 = in.aux_slice1;
2194
2195 for(uword i=0; i<in_n_slices; ++i)
2196 {
2197 out_mem[i] *= Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
2198 }
2199 }
2200 }
2201 }
2202
2203
2204
2205 //! mat X /= Y.subcube(...)
2206 template<typename eT>
2207 inline
2208 void
div_inplace(Mat<eT> & out,const subview_cube<eT> & in)2209 subview_cube<eT>::div_inplace(Mat<eT>& out, const subview_cube<eT>& in)
2210 {
2211 arma_extra_debug_sigprint();
2212
2213 arma_debug_assert_cube_as_mat(out, in, "element-wise division", true);
2214
2215 const uword in_n_rows = in.n_rows;
2216 const uword in_n_cols = in.n_cols;
2217 const uword in_n_slices = in.n_slices;
2218
2219 const uword out_n_rows = out.n_rows;
2220 const uword out_n_cols = out.n_cols;
2221 const uword out_vec_state = out.vec_state;
2222
2223 if(in_n_slices == 1)
2224 {
2225 if( (arma_config::debug) && ((out_n_rows != in_n_rows) || (out_n_cols != in_n_cols)) )
2226 {
2227 std::ostringstream tmp;
2228
2229 tmp
2230 << "in-place element-wise division: "
2231 << out_n_rows << 'x' << out_n_cols << " output matrix is incompatible with "
2232 << in_n_rows << 'x' << in_n_cols << 'x' << in_n_slices << " cube interpreted as "
2233 << in_n_rows << 'x' << in_n_cols << " matrix";
2234
2235 arma_stop_logic_error(tmp.str());
2236 }
2237
2238 for(uword col=0; col < in_n_cols; ++col)
2239 {
2240 arrayops::inplace_div( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
2241 }
2242 }
2243 else
2244 {
2245 if(out_vec_state == 0)
2246 {
2247 if( (in_n_rows == out_n_rows) && (in_n_cols == 1) && (in_n_slices == out_n_cols) )
2248 {
2249 for(uword i=0; i < in_n_slices; ++i)
2250 {
2251 arrayops::inplace_div( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
2252 }
2253 }
2254 else
2255 if( (in_n_rows == 1) && (in_n_cols == out_n_rows) && (in_n_slices == out_n_cols) )
2256 {
2257 const Cube<eT>& Q = in.m;
2258
2259 const uword in_aux_row1 = in.aux_row1;
2260 const uword in_aux_col1 = in.aux_col1;
2261 const uword in_aux_slice1 = in.aux_slice1;
2262
2263 for(uword slice=0; slice < in_n_slices; ++slice)
2264 {
2265 const uword mod_slice = in_aux_slice1 + slice;
2266
2267 eT* out_colptr = out.colptr(slice);
2268
2269 uword i,j;
2270 for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
2271 {
2272 const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
2273 const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
2274
2275 out_colptr[i] /= tmp_i;
2276 out_colptr[j] /= tmp_j;
2277 }
2278
2279 if(i < in_n_cols)
2280 {
2281 out_colptr[i] /= Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
2282 }
2283 }
2284 }
2285 }
2286 else
2287 {
2288 eT* out_mem = out.memptr();
2289
2290 const Cube<eT>& Q = in.m;
2291
2292 const uword in_aux_row1 = in.aux_row1;
2293 const uword in_aux_col1 = in.aux_col1;
2294 const uword in_aux_slice1 = in.aux_slice1;
2295
2296 for(uword i=0; i<in_n_slices; ++i)
2297 {
2298 out_mem[i] /= Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
2299 }
2300 }
2301 }
2302 }
2303
2304
2305
2306 template<typename eT>
2307 inline
2308 typename subview_cube<eT>::iterator
begin()2309 subview_cube<eT>::begin()
2310 {
2311 return iterator(*this, aux_row1, aux_col1, aux_slice1);
2312 }
2313
2314
2315
2316 template<typename eT>
2317 inline
2318 typename subview_cube<eT>::const_iterator
begin() const2319 subview_cube<eT>::begin() const
2320 {
2321 return const_iterator(*this, aux_row1, aux_col1, aux_slice1);
2322 }
2323
2324
2325
2326 template<typename eT>
2327 inline
2328 typename subview_cube<eT>::const_iterator
cbegin() const2329 subview_cube<eT>::cbegin() const
2330 {
2331 return const_iterator(*this, aux_row1, aux_col1, aux_slice1);
2332 }
2333
2334
2335
2336 template<typename eT>
2337 inline
2338 typename subview_cube<eT>::iterator
end()2339 subview_cube<eT>::end()
2340 {
2341 return iterator(*this, aux_row1, aux_col1, aux_slice1 + n_slices);
2342 }
2343
2344
2345
2346 template<typename eT>
2347 inline
2348 typename subview_cube<eT>::const_iterator
end() const2349 subview_cube<eT>::end() const
2350 {
2351 return const_iterator(*this, aux_row1, aux_col1, aux_slice1 + n_slices);
2352 }
2353
2354
2355
2356 template<typename eT>
2357 inline
2358 typename subview_cube<eT>::const_iterator
cend() const2359 subview_cube<eT>::cend() const
2360 {
2361 return const_iterator(*this, aux_row1, aux_col1, aux_slice1 + n_slices);
2362 }
2363
2364
2365
2366 //
2367 //
2368 //
2369
2370
2371
2372 template<typename eT>
2373 inline
iterator()2374 subview_cube<eT>::iterator::iterator()
2375 : M (nullptr)
2376 , current_ptr (nullptr)
2377 , current_row (0 )
2378 , current_col (0 )
2379 , current_slice(0 )
2380 , aux_row1 (0 )
2381 , aux_col1 (0 )
2382 , aux_row2_p1 (0 )
2383 , aux_col2_p1 (0 )
2384 {
2385 arma_extra_debug_sigprint();
2386 // Technically this iterator is invalid (it does not point to a valid element)
2387 }
2388
2389
2390
2391 template<typename eT>
2392 inline
iterator(const iterator & X)2393 subview_cube<eT>::iterator::iterator(const iterator& X)
2394 : M (X.M )
2395 , current_ptr (X.current_ptr )
2396 , current_row (X.current_row )
2397 , current_col (X.current_col )
2398 , current_slice(X.current_slice)
2399 , aux_row1 (X.aux_row1 )
2400 , aux_col1 (X.aux_col1 )
2401 , aux_row2_p1 (X.aux_row2_p1 )
2402 , aux_col2_p1 (X.aux_col2_p1 )
2403 {
2404 arma_extra_debug_sigprint();
2405 }
2406
2407
2408
2409 template<typename eT>
2410 inline
iterator(subview_cube<eT> & in_sv,const uword in_row,const uword in_col,const uword in_slice)2411 subview_cube<eT>::iterator::iterator(subview_cube<eT>& in_sv, const uword in_row, const uword in_col, const uword in_slice)
2412 : M (&(const_cast< Cube<eT>& >(in_sv.m)))
2413 , current_ptr (&(M->at(in_row,in_col,in_slice)) )
2414 , current_row (in_row )
2415 , current_col (in_col )
2416 , current_slice(in_slice )
2417 , aux_row1 (in_sv.aux_row1 )
2418 , aux_col1 (in_sv.aux_col1 )
2419 , aux_row2_p1 (in_sv.aux_row1 + in_sv.n_rows )
2420 , aux_col2_p1 (in_sv.aux_col1 + in_sv.n_cols )
2421 {
2422 arma_extra_debug_sigprint();
2423 }
2424
2425
2426
2427 template<typename eT>
2428 inline
2429 arma_warn_unused
2430 eT&
operator *()2431 subview_cube<eT>::iterator::operator*()
2432 {
2433 return (*current_ptr);
2434 }
2435
2436
2437
2438 template<typename eT>
2439 inline
2440 typename subview_cube<eT>::iterator&
operator ++()2441 subview_cube<eT>::iterator::operator++()
2442 {
2443 current_row++;
2444
2445 if(current_row == aux_row2_p1)
2446 {
2447 current_row = aux_row1;
2448 current_col++;
2449
2450 if(current_col == aux_col2_p1)
2451 {
2452 current_col = aux_col1;
2453 current_slice++;
2454 }
2455
2456 current_ptr = &( (*M).at(current_row,current_col,current_slice) );
2457 }
2458 else
2459 {
2460 current_ptr++;
2461 }
2462
2463 return *this;
2464 }
2465
2466
2467
2468 template<typename eT>
2469 inline
2470 arma_warn_unused
2471 typename subview_cube<eT>::iterator
operator ++(int)2472 subview_cube<eT>::iterator::operator++(int)
2473 {
2474 typename subview_cube<eT>::iterator temp(*this);
2475
2476 ++(*this);
2477
2478 return temp;
2479 }
2480
2481
2482
2483 template<typename eT>
2484 inline
2485 arma_warn_unused
2486 bool
operator ==(const iterator & rhs) const2487 subview_cube<eT>::iterator::operator==(const iterator& rhs) const
2488 {
2489 return (current_ptr == rhs.current_ptr);
2490 }
2491
2492
2493
2494 template<typename eT>
2495 inline
2496 arma_warn_unused
2497 bool
operator !=(const iterator & rhs) const2498 subview_cube<eT>::iterator::operator!=(const iterator& rhs) const
2499 {
2500 return (current_ptr != rhs.current_ptr);
2501 }
2502
2503
2504
2505 template<typename eT>
2506 inline
2507 arma_warn_unused
2508 bool
operator ==(const const_iterator & rhs) const2509 subview_cube<eT>::iterator::operator==(const const_iterator& rhs) const
2510 {
2511 return (current_ptr == rhs.current_ptr);
2512 }
2513
2514
2515
2516 template<typename eT>
2517 inline
2518 arma_warn_unused
2519 bool
operator !=(const const_iterator & rhs) const2520 subview_cube<eT>::iterator::operator!=(const const_iterator& rhs) const
2521 {
2522 return (current_ptr != rhs.current_ptr);
2523 }
2524
2525
2526
2527 //
2528 //
2529 //
2530
2531
2532
2533 template<typename eT>
2534 inline
const_iterator()2535 subview_cube<eT>::const_iterator::const_iterator()
2536 : M (nullptr)
2537 , current_ptr (nullptr)
2538 , current_row (0 )
2539 , current_col (0 )
2540 , current_slice(0 )
2541 , aux_row1 (0 )
2542 , aux_col1 (0 )
2543 , aux_row2_p1 (0 )
2544 , aux_col2_p1 (0 )
2545 {
2546 arma_extra_debug_sigprint();
2547 // Technically this iterator is invalid (it does not point to a valid element)
2548 }
2549
2550
2551
2552 template<typename eT>
2553 inline
const_iterator(const iterator & X)2554 subview_cube<eT>::const_iterator::const_iterator(const iterator& X)
2555 : M (X.M )
2556 , current_ptr (X.current_ptr )
2557 , current_row (X.current_row )
2558 , current_col (X.current_col )
2559 , current_slice(X.current_slice)
2560 , aux_row1 (X.aux_row1 )
2561 , aux_col1 (X.aux_col1 )
2562 , aux_row2_p1 (X.aux_row2_p1 )
2563 , aux_col2_p1 (X.aux_col2_p1 )
2564 {
2565 arma_extra_debug_sigprint();
2566 }
2567
2568
2569
2570 template<typename eT>
2571 inline
const_iterator(const const_iterator & X)2572 subview_cube<eT>::const_iterator::const_iterator(const const_iterator& X)
2573 : M (X.M )
2574 , current_ptr (X.current_ptr )
2575 , current_row (X.current_row )
2576 , current_col (X.current_col )
2577 , current_slice(X.current_slice)
2578 , aux_row1 (X.aux_row1 )
2579 , aux_col1 (X.aux_col1 )
2580 , aux_row2_p1 (X.aux_row2_p1 )
2581 , aux_col2_p1 (X.aux_col2_p1 )
2582 {
2583 arma_extra_debug_sigprint();
2584 }
2585
2586
2587
2588 template<typename eT>
2589 inline
const_iterator(const subview_cube<eT> & in_sv,const uword in_row,const uword in_col,const uword in_slice)2590 subview_cube<eT>::const_iterator::const_iterator(const subview_cube<eT>& in_sv, const uword in_row, const uword in_col, const uword in_slice)
2591 : M (&(in_sv.m) )
2592 , current_ptr (&(M->at(in_row,in_col,in_slice)))
2593 , current_row (in_row )
2594 , current_col (in_col )
2595 , current_slice(in_slice )
2596 , aux_row1 (in_sv.aux_row1 )
2597 , aux_col1 (in_sv.aux_col1 )
2598 , aux_row2_p1 (in_sv.aux_row1 + in_sv.n_rows )
2599 , aux_col2_p1 (in_sv.aux_col1 + in_sv.n_cols )
2600 {
2601 arma_extra_debug_sigprint();
2602 }
2603
2604
2605
2606 template<typename eT>
2607 inline
2608 arma_warn_unused
2609 const eT&
operator *()2610 subview_cube<eT>::const_iterator::operator*()
2611 {
2612 return (*current_ptr);
2613 }
2614
2615
2616
2617 template<typename eT>
2618 inline
2619 typename subview_cube<eT>::const_iterator&
operator ++()2620 subview_cube<eT>::const_iterator::operator++()
2621 {
2622 current_row++;
2623
2624 if(current_row == aux_row2_p1)
2625 {
2626 current_row = aux_row1;
2627 current_col++;
2628
2629 if(current_col == aux_col2_p1)
2630 {
2631 current_col = aux_col1;
2632 current_slice++;
2633 }
2634
2635 current_ptr = &( (*M).at(current_row,current_col,current_slice) );
2636 }
2637 else
2638 {
2639 current_ptr++;
2640 }
2641
2642 return *this;
2643 }
2644
2645
2646
2647 template<typename eT>
2648 inline
2649 arma_warn_unused
2650 typename subview_cube<eT>::const_iterator
operator ++(int)2651 subview_cube<eT>::const_iterator::operator++(int)
2652 {
2653 typename subview_cube<eT>::const_iterator temp(*this);
2654
2655 ++(*this);
2656
2657 return temp;
2658 }
2659
2660
2661
2662 template<typename eT>
2663 inline
2664 arma_warn_unused
2665 bool
operator ==(const iterator & rhs) const2666 subview_cube<eT>::const_iterator::operator==(const iterator& rhs) const
2667 {
2668 return (current_ptr == rhs.current_ptr);
2669 }
2670
2671
2672
2673 template<typename eT>
2674 inline
2675 arma_warn_unused
2676 bool
operator !=(const iterator & rhs) const2677 subview_cube<eT>::const_iterator::operator!=(const iterator& rhs) const
2678 {
2679 return (current_ptr != rhs.current_ptr);
2680 }
2681
2682
2683
2684 template<typename eT>
2685 inline
2686 arma_warn_unused
2687 bool
operator ==(const const_iterator & rhs) const2688 subview_cube<eT>::const_iterator::operator==(const const_iterator& rhs) const
2689 {
2690 return (current_ptr == rhs.current_ptr);
2691 }
2692
2693
2694
2695 template<typename eT>
2696 inline
2697 arma_warn_unused
2698 bool
operator !=(const const_iterator & rhs) const2699 subview_cube<eT>::const_iterator::operator!=(const const_iterator& rhs) const
2700 {
2701 return (current_ptr != rhs.current_ptr);
2702 }
2703
2704
2705
2706 //! @}
2707