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