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 spdiagview
20 //! @{
21
22
23 template<typename eT>
24 inline
~spdiagview()25 spdiagview<eT>::~spdiagview()
26 {
27 arma_extra_debug_sigprint();
28 }
29
30
31 template<typename eT>
32 arma_inline
spdiagview(const SpMat<eT> & in_m,const uword in_row_offset,const uword in_col_offset,const uword in_len)33 spdiagview<eT>::spdiagview(const SpMat<eT>& in_m, const uword in_row_offset, const uword in_col_offset, const uword in_len)
34 : m(in_m)
35 , row_offset(in_row_offset)
36 , col_offset(in_col_offset)
37 , n_rows(in_len)
38 , n_elem(in_len)
39 {
40 arma_extra_debug_sigprint();
41 }
42
43
44
45 //! set a diagonal of our matrix using a diagonal from a foreign matrix
46 template<typename eT>
47 inline
48 void
operator =(const spdiagview<eT> & x)49 spdiagview<eT>::operator= (const spdiagview<eT>& x)
50 {
51 arma_extra_debug_sigprint();
52
53 spdiagview<eT>& d = *this;
54
55 arma_debug_check( (d.n_elem != x.n_elem), "spdiagview: diagonals have incompatible lengths" );
56
57 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
58 const SpMat<eT>& x_m = x.m;
59
60 if( (&d_m == &x_m) || ((d.row_offset == 0) && (d.col_offset == 0)) )
61 {
62 const Mat<eT> tmp(x);
63
64 (*this).operator=(tmp);
65 }
66 else
67 {
68 const uword d_n_elem = d.n_elem;
69 const uword d_row_offset = d.row_offset;
70 const uword d_col_offset = d.col_offset;
71
72 const uword x_row_offset = x.row_offset;
73 const uword x_col_offset = x.col_offset;
74
75 for(uword i=0; i < d_n_elem; ++i)
76 {
77 d_m.at(i + d_row_offset, i + d_col_offset) = x_m.at(i + x_row_offset, i + x_col_offset);
78 }
79 }
80 }
81
82
83
84 template<typename eT>
85 inline
86 void
operator +=(const eT val)87 spdiagview<eT>::operator+=(const eT val)
88 {
89 arma_extra_debug_sigprint();
90
91 if(val == eT(0)) { return; }
92
93 SpMat<eT>& t_m = const_cast< SpMat<eT>& >(m);
94
95 const uword t_n_elem = n_elem;
96 const uword t_row_offset = row_offset;
97 const uword t_col_offset = col_offset;
98
99 for(uword i=0; i < t_n_elem; ++i)
100 {
101 t_m.at(i + t_row_offset, i + t_col_offset) += val;
102 }
103 }
104
105
106
107 template<typename eT>
108 inline
109 void
operator -=(const eT val)110 spdiagview<eT>::operator-=(const eT val)
111 {
112 arma_extra_debug_sigprint();
113
114 if(val == eT(0)) { return; }
115
116 SpMat<eT>& t_m = const_cast< SpMat<eT>& >(m);
117
118 const uword t_n_elem = n_elem;
119 const uword t_row_offset = row_offset;
120 const uword t_col_offset = col_offset;
121
122 for(uword i=0; i < t_n_elem; ++i)
123 {
124 t_m.at(i + t_row_offset, i + t_col_offset) -= val;
125 }
126 }
127
128
129
130 template<typename eT>
131 inline
132 void
operator *=(const eT val)133 spdiagview<eT>::operator*=(const eT val)
134 {
135 arma_extra_debug_sigprint();
136
137 if(val == eT(0)) { (*this).zeros(); return; }
138
139 SpMat<eT>& t_m = const_cast< SpMat<eT>& >(m);
140
141 const uword t_n_elem = n_elem;
142 const uword t_row_offset = row_offset;
143 const uword t_col_offset = col_offset;
144
145 for(uword i=0; i < t_n_elem; ++i)
146 {
147 t_m.at(i + t_row_offset, i + t_col_offset) *= val;
148 }
149 }
150
151
152
153 template<typename eT>
154 inline
155 void
operator /=(const eT val)156 spdiagview<eT>::operator/=(const eT val)
157 {
158 arma_extra_debug_sigprint();
159
160 SpMat<eT>& t_m = const_cast< SpMat<eT>& >(m);
161
162 const uword t_n_elem = n_elem;
163 const uword t_row_offset = row_offset;
164 const uword t_col_offset = col_offset;
165
166 for(uword i=0; i < t_n_elem; ++i)
167 {
168 t_m.at(i + t_row_offset, i + t_col_offset) /= val;
169 }
170 }
171
172
173
174 //! set a diagonal of our matrix using data from a foreign object
175 template<typename eT>
176 template<typename T1>
177 inline
178 void
operator =(const Base<eT,T1> & o)179 spdiagview<eT>::operator= (const Base<eT,T1>& o)
180 {
181 arma_extra_debug_sigprint();
182
183 spdiagview<eT>& d = *this;
184
185 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
186
187 const uword d_n_elem = d.n_elem;
188 const uword d_row_offset = d.row_offset;
189 const uword d_col_offset = d.col_offset;
190
191 if(is_same_type< T1, Gen<Col<eT>, gen_zeros> >::yes)
192 {
193 const Proxy<T1> P(o.get_ref());
194
195 arma_debug_check( (d_n_elem != P.get_n_elem()), "spdiagview: given object has incompatible size" );
196
197 (*this).zeros();
198
199 return;
200 }
201
202 if(is_same_type< T1, Gen<Col<eT>, gen_ones> >::yes)
203 {
204 const Proxy<T1> P(o.get_ref());
205
206 arma_debug_check( (d_n_elem != P.get_n_elem()), "spdiagview: given object has incompatible size" );
207
208 (*this).ones();
209
210 return;
211 }
212
213 const quasi_unwrap<T1> U(o.get_ref());
214 const Mat<eT>& x = U.M;
215
216 const eT* x_mem = x.memptr();
217
218 arma_debug_check
219 (
220 ( (d_n_elem != x.n_elem) || ((x.n_rows != 1) && (x.n_cols != 1)) ),
221 "spdiagview: given object has incompatible size"
222 );
223
224 if( (d_row_offset == 0) && (d_col_offset == 0) )
225 {
226 SpMat<eT> tmp1;
227
228 tmp1.eye(d_m.n_rows, d_m.n_cols);
229
230 bool has_zero = false;
231
232 for(uword i=0; i < d_n_elem; ++i)
233 {
234 const eT val = x_mem[i];
235
236 access::rw(tmp1.values[i]) = val;
237
238 if(val == eT(0)) { has_zero = true; }
239 }
240
241 if(has_zero) { tmp1.remove_zeros(); }
242
243 if(tmp1.n_nonzero == 0) { (*this).zeros(); return; }
244
245 SpMat<eT> tmp2;
246
247 spglue_merge::diagview_merge(tmp2, d_m, tmp1);
248
249 d_m.steal_mem(tmp2);
250 }
251 else
252 {
253 for(uword i=0; i < d_n_elem; ++i)
254 {
255 d_m.at(i + d_row_offset, i + d_col_offset) = x_mem[i];
256 }
257 }
258 }
259
260
261
262 template<typename eT>
263 template<typename T1>
264 inline
265 void
operator +=(const Base<eT,T1> & o)266 spdiagview<eT>::operator+=(const Base<eT,T1>& o)
267 {
268 arma_extra_debug_sigprint();
269
270 spdiagview<eT>& d = *this;
271
272 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
273
274 const uword d_n_elem = d.n_elem;
275 const uword d_row_offset = d.row_offset;
276 const uword d_col_offset = d.col_offset;
277
278 const Proxy<T1> P( o.get_ref() );
279
280 arma_debug_check
281 (
282 ( (d_n_elem != P.get_n_elem()) || ((P.get_n_rows() != 1) && (P.get_n_cols() != 1)) ),
283 "spdiagview: given object has incompatible size"
284 );
285
286 if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) )
287 {
288 const unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
289 const Mat<eT>& x = tmp.M;
290
291 const eT* x_mem = x.memptr();
292
293 for(uword i=0; i < d_n_elem; ++i)
294 {
295 d_m.at(i + d_row_offset, i + d_col_offset) += x_mem[i];
296 }
297 }
298 else
299 {
300 typename Proxy<T1>::ea_type Pea = P.get_ea();
301
302 for(uword i=0; i < d_n_elem; ++i)
303 {
304 d_m.at(i + d_row_offset, i + d_col_offset) += Pea[i];
305 }
306 }
307 }
308
309
310
311 template<typename eT>
312 template<typename T1>
313 inline
314 void
operator -=(const Base<eT,T1> & o)315 spdiagview<eT>::operator-=(const Base<eT,T1>& o)
316 {
317 arma_extra_debug_sigprint();
318
319 spdiagview<eT>& d = *this;
320
321 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
322
323 const uword d_n_elem = d.n_elem;
324 const uword d_row_offset = d.row_offset;
325 const uword d_col_offset = d.col_offset;
326
327 const Proxy<T1> P( o.get_ref() );
328
329 arma_debug_check
330 (
331 ( (d_n_elem != P.get_n_elem()) || ((P.get_n_rows() != 1) && (P.get_n_cols() != 1)) ),
332 "spdiagview: given object has incompatible size"
333 );
334
335 if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) )
336 {
337 const unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
338 const Mat<eT>& x = tmp.M;
339
340 const eT* x_mem = x.memptr();
341
342 for(uword i=0; i < d_n_elem; ++i)
343 {
344 d_m.at(i + d_row_offset, i + d_col_offset) -= x_mem[i];
345 }
346 }
347 else
348 {
349 typename Proxy<T1>::ea_type Pea = P.get_ea();
350
351 for(uword i=0; i < d_n_elem; ++i)
352 {
353 d_m.at(i + d_row_offset, i + d_col_offset) -= Pea[i];
354 }
355 }
356 }
357
358
359
360 template<typename eT>
361 template<typename T1>
362 inline
363 void
operator %=(const Base<eT,T1> & o)364 spdiagview<eT>::operator%=(const Base<eT,T1>& o)
365 {
366 arma_extra_debug_sigprint();
367
368 spdiagview<eT>& d = *this;
369
370 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
371
372 const uword d_n_elem = d.n_elem;
373 const uword d_row_offset = d.row_offset;
374 const uword d_col_offset = d.col_offset;
375
376 const Proxy<T1> P( o.get_ref() );
377
378 arma_debug_check
379 (
380 ( (d_n_elem != P.get_n_elem()) || ((P.get_n_rows() != 1) && (P.get_n_cols() != 1)) ),
381 "spdiagview: given object has incompatible size"
382 );
383
384 if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) )
385 {
386 const unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
387 const Mat<eT>& x = tmp.M;
388
389 const eT* x_mem = x.memptr();
390
391 for(uword i=0; i < d_n_elem; ++i)
392 {
393 d_m.at(i + d_row_offset, i + d_col_offset) *= x_mem[i];
394 }
395 }
396 else
397 {
398 typename Proxy<T1>::ea_type Pea = P.get_ea();
399
400 for(uword i=0; i < d_n_elem; ++i)
401 {
402 d_m.at(i + d_row_offset, i + d_col_offset) *= Pea[i];
403 }
404 }
405 }
406
407
408
409 template<typename eT>
410 template<typename T1>
411 inline
412 void
operator /=(const Base<eT,T1> & o)413 spdiagview<eT>::operator/=(const Base<eT,T1>& o)
414 {
415 arma_extra_debug_sigprint();
416
417 spdiagview<eT>& d = *this;
418
419 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
420
421 const uword d_n_elem = d.n_elem;
422 const uword d_row_offset = d.row_offset;
423 const uword d_col_offset = d.col_offset;
424
425 const Proxy<T1> P( o.get_ref() );
426
427 arma_debug_check
428 (
429 ( (d_n_elem != P.get_n_elem()) || ((P.get_n_rows() != 1) && (P.get_n_cols() != 1)) ),
430 "spdiagview: given object has incompatible size"
431 );
432
433 if( (is_Mat<typename Proxy<T1>::stored_type>::value) || (Proxy<T1>::use_at) )
434 {
435 const unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
436 const Mat<eT>& x = tmp.M;
437
438 const eT* x_mem = x.memptr();
439
440 for(uword i=0; i < d_n_elem; ++i)
441 {
442 d_m.at(i + d_row_offset, i + d_col_offset) /= x_mem[i];
443 }
444 }
445 else
446 {
447 typename Proxy<T1>::ea_type Pea = P.get_ea();
448
449 for(uword i=0; i < d_n_elem; ++i)
450 {
451 d_m.at(i + d_row_offset, i + d_col_offset) /= Pea[i];
452 }
453 }
454 }
455
456
457
458 //! set a diagonal of our matrix using data from a foreign object
459 template<typename eT>
460 template<typename T1>
461 inline
462 void
operator =(const SpBase<eT,T1> & o)463 spdiagview<eT>::operator= (const SpBase<eT,T1>& o)
464 {
465 arma_extra_debug_sigprint();
466
467 const unwrap_spmat<T1> U( o.get_ref() );
468 const SpMat<eT>& x = U.M;
469
470 arma_debug_check
471 (
472 ( (n_elem != x.n_elem) || ((x.n_rows != 1) && (x.n_cols != 1)) ),
473 "spdiagview: given object has incompatible size"
474 );
475
476 const Mat<eT> tmp(x);
477
478 (*this).operator=(tmp);
479 }
480
481
482
483 template<typename eT>
484 template<typename T1>
485 inline
486 void
operator +=(const SpBase<eT,T1> & o)487 spdiagview<eT>::operator+=(const SpBase<eT,T1>& o)
488 {
489 arma_extra_debug_sigprint();
490
491 spdiagview<eT>& d = *this;
492
493 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
494
495 const uword d_n_elem = d.n_elem;
496 const uword d_row_offset = d.row_offset;
497 const uword d_col_offset = d.col_offset;
498
499 const SpProxy<T1> P( o.get_ref() );
500
501 arma_debug_check
502 (
503 ( (d_n_elem != P.get_n_elem()) || ((P.get_n_rows() != 1) && (P.get_n_cols() != 1)) ),
504 "spdiagview: given object has incompatible size"
505 );
506
507 if( SpProxy<T1>::use_iterator || P.is_alias(d_m) )
508 {
509 const SpMat<eT> tmp(P.Q);
510
511 if(tmp.n_cols == 1)
512 {
513 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) += tmp.at(i,0); }
514 }
515 else
516 if(tmp.n_rows == 1)
517 {
518 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) += tmp.at(0,i); }
519 }
520 }
521 else
522 {
523 if(P.get_n_cols() == 1)
524 {
525 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) += P.at(i,0); }
526 }
527 else
528 if(P.get_n_rows() == 1)
529 {
530 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) += P.at(0,i); }
531 }
532 }
533 }
534
535
536
537 template<typename eT>
538 template<typename T1>
539 inline
540 void
operator -=(const SpBase<eT,T1> & o)541 spdiagview<eT>::operator-=(const SpBase<eT,T1>& o)
542 {
543 arma_extra_debug_sigprint();
544
545 spdiagview<eT>& d = *this;
546
547 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
548
549 const uword d_n_elem = d.n_elem;
550 const uword d_row_offset = d.row_offset;
551 const uword d_col_offset = d.col_offset;
552
553 const SpProxy<T1> P( o.get_ref() );
554
555 arma_debug_check
556 (
557 ( (d_n_elem != P.get_n_elem()) || ((P.get_n_rows() != 1) && (P.get_n_cols() != 1)) ),
558 "spdiagview: given object has incompatible size"
559 );
560
561 if( SpProxy<T1>::use_iterator || P.is_alias(d_m) )
562 {
563 const SpMat<eT> tmp(P.Q);
564
565 if(tmp.n_cols == 1)
566 {
567 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) -= tmp.at(i,0); }
568 }
569 else
570 if(tmp.n_rows == 1)
571 {
572 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) -= tmp.at(0,i); }
573 }
574 }
575 else
576 {
577 if(P.get_n_cols() == 1)
578 {
579 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) -= P.at(i,0); }
580 }
581 else
582 if(P.get_n_rows() == 1)
583 {
584 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) -= P.at(0,i); }
585 }
586 }
587 }
588
589
590
591 template<typename eT>
592 template<typename T1>
593 inline
594 void
operator %=(const SpBase<eT,T1> & o)595 spdiagview<eT>::operator%=(const SpBase<eT,T1>& o)
596 {
597 arma_extra_debug_sigprint();
598
599 spdiagview<eT>& d = *this;
600
601 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
602
603 const uword d_n_elem = d.n_elem;
604 const uword d_row_offset = d.row_offset;
605 const uword d_col_offset = d.col_offset;
606
607 const SpProxy<T1> P( o.get_ref() );
608
609 arma_debug_check
610 (
611 ( (d_n_elem != P.get_n_elem()) || ((P.get_n_rows() != 1) && (P.get_n_cols() != 1)) ),
612 "spdiagview: given object has incompatible size"
613 );
614
615 if( SpProxy<T1>::use_iterator || P.is_alias(d_m) )
616 {
617 const SpMat<eT> tmp(P.Q);
618
619 if(tmp.n_cols == 1)
620 {
621 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) *= tmp.at(i,0); }
622 }
623 else
624 if(tmp.n_rows == 1)
625 {
626 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) *= tmp.at(0,i); }
627 }
628 }
629 else
630 {
631 if(P.get_n_cols() == 1)
632 {
633 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) *= P.at(i,0); }
634 }
635 else
636 if(P.get_n_rows() == 1)
637 {
638 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) *= P.at(0,i); }
639 }
640 }
641 }
642
643
644
645 template<typename eT>
646 template<typename T1>
647 inline
648 void
operator /=(const SpBase<eT,T1> & o)649 spdiagview<eT>::operator/=(const SpBase<eT,T1>& o)
650 {
651 arma_extra_debug_sigprint();
652
653 spdiagview<eT>& d = *this;
654
655 SpMat<eT>& d_m = const_cast< SpMat<eT>& >(d.m);
656
657 const uword d_n_elem = d.n_elem;
658 const uword d_row_offset = d.row_offset;
659 const uword d_col_offset = d.col_offset;
660
661 const SpProxy<T1> P( o.get_ref() );
662
663 arma_debug_check
664 (
665 ( (d_n_elem != P.get_n_elem()) || ((P.get_n_rows() != 1) && (P.get_n_cols() != 1)) ),
666 "spdiagview: given object has incompatible size"
667 );
668
669 if( SpProxy<T1>::use_iterator || P.is_alias(d_m) )
670 {
671 const SpMat<eT> tmp(P.Q);
672
673 if(tmp.n_cols == 1)
674 {
675 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) /= tmp.at(i,0); }
676 }
677 else
678 if(tmp.n_rows == 1)
679 {
680 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) /= tmp.at(0,i); }
681 }
682 }
683 else
684 {
685 if(P.get_n_cols() == 1)
686 {
687 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) /= P.at(i,0); }
688 }
689 else
690 if(P.get_n_rows() == 1)
691 {
692 for(uword i=0; i < d_n_elem; ++i) { d_m.at(i + d_row_offset, i + d_col_offset) /= P.at(0,i); }
693 }
694 }
695 }
696
697
698
699 template<typename eT>
700 inline
701 void
extract(SpMat<eT> & out,const spdiagview<eT> & d)702 spdiagview<eT>::extract(SpMat<eT>& out, const spdiagview<eT>& d)
703 {
704 arma_extra_debug_sigprint();
705
706 const SpMat<eT>& d_m = d.m;
707
708 const uword d_n_elem = d.n_elem;
709 const uword d_row_offset = d.row_offset;
710 const uword d_col_offset = d.col_offset;
711
712 Col<eT> cache(d_n_elem, arma_nozeros_indicator());
713 eT* cache_mem = cache.memptr();
714
715 uword d_n_nonzero = 0;
716
717 for(uword i=0; i < d_n_elem; ++i)
718 {
719 const eT val = d_m.at(i + d_row_offset, i + d_col_offset);
720
721 cache_mem[i] = val;
722
723 d_n_nonzero += (val != eT(0)) ? uword(1) : uword(0);
724 }
725
726 out.reserve(d_n_elem, 1, d_n_nonzero);
727
728 uword count = 0;
729 for(uword i=0; i < d_n_elem; ++i)
730 {
731 const eT val = cache_mem[i];
732
733 if(val != eT(0))
734 {
735 access::rw(out.row_indices[count]) = i;
736 access::rw(out.values[count]) = val;
737 ++count;
738 }
739 }
740
741 access::rw(out.col_ptrs[0]) = 0;
742 access::rw(out.col_ptrs[1]) = d_n_nonzero;
743 }
744
745
746
747 //! extract a diagonal and store it as a dense column vector
748 template<typename eT>
749 inline
750 void
extract(Mat<eT> & out,const spdiagview<eT> & in)751 spdiagview<eT>::extract(Mat<eT>& out, const spdiagview<eT>& in)
752 {
753 arma_extra_debug_sigprint();
754
755 // NOTE: we're assuming that the 'out' matrix has already been set to the correct size;
756 // size setting is done by either the Mat contructor or Mat::operator=()
757
758 const SpMat<eT>& in_m = in.m;
759
760 const uword in_n_elem = in.n_elem;
761 const uword in_row_offset = in.row_offset;
762 const uword in_col_offset = in.col_offset;
763
764 eT* out_mem = out.memptr();
765
766 for(uword i=0; i < in_n_elem; ++i)
767 {
768 out_mem[i] = in_m.at(i + in_row_offset, i + in_col_offset);
769 }
770 }
771
772
773
774 template<typename eT>
775 inline
776 SpMat_MapMat_val<eT>
operator [](const uword i)777 spdiagview<eT>::operator[](const uword i)
778 {
779 return (const_cast< SpMat<eT>& >(m)).at(i+row_offset, i+col_offset);
780 }
781
782
783
784 template<typename eT>
785 inline
786 eT
operator [](const uword i) const787 spdiagview<eT>::operator[](const uword i) const
788 {
789 return m.at(i+row_offset, i+col_offset);
790 }
791
792
793
794 template<typename eT>
795 inline
796 SpMat_MapMat_val<eT>
at(const uword i)797 spdiagview<eT>::at(const uword i)
798 {
799 return (const_cast< SpMat<eT>& >(m)).at(i+row_offset, i+col_offset);
800 }
801
802
803
804 template<typename eT>
805 inline
806 eT
at(const uword i) const807 spdiagview<eT>::at(const uword i) const
808 {
809 return m.at(i+row_offset, i+col_offset);
810 }
811
812
813
814 template<typename eT>
815 inline
816 SpMat_MapMat_val<eT>
operator ()(const uword i)817 spdiagview<eT>::operator()(const uword i)
818 {
819 arma_debug_check_bounds( (i >= n_elem), "spdiagview::operator(): out of bounds" );
820
821 return (const_cast< SpMat<eT>& >(m)).at(i+row_offset, i+col_offset);
822 }
823
824
825
826 template<typename eT>
827 inline
828 eT
operator ()(const uword i) const829 spdiagview<eT>::operator()(const uword i) const
830 {
831 arma_debug_check_bounds( (i >= n_elem), "spdiagview::operator(): out of bounds" );
832
833 return m.at(i+row_offset, i+col_offset);
834 }
835
836
837
838 template<typename eT>
839 inline
840 SpMat_MapMat_val<eT>
at(const uword row,const uword)841 spdiagview<eT>::at(const uword row, const uword)
842 {
843 return (const_cast< SpMat<eT>& >(m)).at(row+row_offset, row+col_offset);
844 }
845
846
847
848 template<typename eT>
849 inline
850 eT
at(const uword row,const uword) const851 spdiagview<eT>::at(const uword row, const uword) const
852 {
853 return m.at(row+row_offset, row+col_offset);
854 }
855
856
857
858 template<typename eT>
859 inline
860 SpMat_MapMat_val<eT>
operator ()(const uword row,const uword col)861 spdiagview<eT>::operator()(const uword row, const uword col)
862 {
863 arma_debug_check_bounds( ((row >= n_elem) || (col > 0)), "spdiagview::operator(): out of bounds" );
864
865 return (const_cast< SpMat<eT>& >(m)).at(row+row_offset, row+col_offset);
866 }
867
868
869
870 template<typename eT>
871 inline
872 eT
operator ()(const uword row,const uword col) const873 spdiagview<eT>::operator()(const uword row, const uword col) const
874 {
875 arma_debug_check_bounds( ((row >= n_elem) || (col > 0)), "spdiagview::operator(): out of bounds" );
876
877 return m.at(row+row_offset, row+col_offset);
878 }
879
880
881
882 template<typename eT>
883 inline
884 void
replace(const eT old_val,const eT new_val)885 spdiagview<eT>::replace(const eT old_val, const eT new_val)
886 {
887 arma_extra_debug_sigprint();
888
889 if(old_val == eT(0))
890 {
891 arma_debug_warn_level(1, "spdiagview::replace(): replacement not done, as old_val = 0");
892 }
893 else
894 {
895 Mat<eT> tmp(*this);
896
897 tmp.replace(old_val, new_val);
898
899 (*this).operator=(tmp);
900 }
901 }
902
903
904
905 template<typename eT>
906 inline
907 void
clean(const typename get_pod_type<eT>::result threshold)908 spdiagview<eT>::clean(const typename get_pod_type<eT>::result threshold)
909 {
910 arma_extra_debug_sigprint();
911
912 Mat<eT> tmp(*this);
913
914 tmp.clean(threshold);
915
916 (*this).operator=(tmp);
917 }
918
919
920
921 template<typename eT>
922 inline
923 void
clamp(const eT min_val,const eT max_val)924 spdiagview<eT>::clamp(const eT min_val, const eT max_val)
925 {
926 arma_extra_debug_sigprint();
927
928 SpMat<eT> tmp(*this);
929
930 tmp.clamp(min_val, max_val);
931
932 (*this).operator=(tmp);
933 }
934
935
936
937 template<typename eT>
938 inline
939 void
fill(const eT val)940 spdiagview<eT>::fill(const eT val)
941 {
942 arma_extra_debug_sigprint();
943
944 if( (row_offset == 0) && (col_offset == 0) && (m.sync_state != 1) )
945 {
946 if(val == eT(0))
947 {
948 SpMat<eT> tmp(arma_reserve_indicator(), m.n_rows, m.n_cols, m.n_nonzero); // worst case scenario
949
950 typename SpMat<eT>::const_iterator it = m.begin();
951 typename SpMat<eT>::const_iterator it_end = m.end();
952
953 uword count = 0;
954
955 for(; it != it_end; ++it)
956 {
957 const uword row = it.row();
958 const uword col = it.col();
959
960 if(row != col)
961 {
962 access::rw(tmp.values[count]) = (*it);
963 access::rw(tmp.row_indices[count]) = row;
964 access::rw(tmp.col_ptrs[col + 1])++;
965 ++count;
966 }
967 }
968
969 for(uword i=0; i < tmp.n_cols; ++i)
970 {
971 access::rw(tmp.col_ptrs[i + 1]) += tmp.col_ptrs[i];
972 }
973
974 // quick resize without reallocating memory and copying data
975 access::rw( tmp.n_nonzero) = count;
976 access::rw( tmp.values[count]) = eT(0);
977 access::rw(tmp.row_indices[count]) = uword(0);
978
979 access::rw(m).steal_mem(tmp);
980 }
981 else // val != eT(0)
982 {
983 SpMat<eT> tmp1;
984
985 tmp1.eye(m.n_rows, m.n_cols);
986
987 if(val != eT(1)) { tmp1 *= val; }
988
989 SpMat<eT> tmp2;
990
991 spglue_merge::diagview_merge(tmp2, m, tmp1);
992
993 access::rw(m).steal_mem(tmp2);
994 }
995 }
996 else
997 {
998 SpMat<eT>& x = const_cast< SpMat<eT>& >(m);
999
1000 const uword local_n_elem = n_elem;
1001
1002 for(uword i=0; i < local_n_elem; ++i)
1003 {
1004 x.at(i+row_offset, i+col_offset) = val;
1005 }
1006 }
1007 }
1008
1009
1010
1011 template<typename eT>
1012 inline
1013 void
zeros()1014 spdiagview<eT>::zeros()
1015 {
1016 arma_extra_debug_sigprint();
1017
1018 (*this).fill(eT(0));
1019 }
1020
1021
1022
1023 template<typename eT>
1024 inline
1025 void
ones()1026 spdiagview<eT>::ones()
1027 {
1028 arma_extra_debug_sigprint();
1029
1030 (*this).fill(eT(1));
1031 }
1032
1033
1034
1035 template<typename eT>
1036 inline
1037 void
randu()1038 spdiagview<eT>::randu()
1039 {
1040 arma_extra_debug_sigprint();
1041
1042 SpMat<eT>& x = const_cast< SpMat<eT>& >(m);
1043
1044 const uword local_n_elem = n_elem;
1045
1046 for(uword i=0; i < local_n_elem; ++i)
1047 {
1048 x.at(i+row_offset, i+col_offset) = eT(arma_rng::randu<eT>());
1049 }
1050 }
1051
1052
1053
1054 template<typename eT>
1055 inline
1056 void
randn()1057 spdiagview<eT>::randn()
1058 {
1059 arma_extra_debug_sigprint();
1060
1061 SpMat<eT>& x = const_cast< SpMat<eT>& >(m);
1062
1063 const uword local_n_elem = n_elem;
1064
1065 for(uword i=0; i < local_n_elem; ++i)
1066 {
1067 x.at(i+row_offset, i+col_offset) = eT(arma_rng::randn<eT>());
1068 }
1069 }
1070
1071
1072
1073 //! @}
1074