|
max@0
|
1 // Copyright (C) 2008-2011 NICTA (www.nicta.com.au)
|
|
max@0
|
2 // Copyright (C) 2008-2011 Conrad Sanderson
|
|
max@0
|
3 //
|
|
max@0
|
4 // This file is part of the Armadillo C++ library.
|
|
max@0
|
5 // It is provided without any warranty of fitness
|
|
max@0
|
6 // for any purpose. You can redistribute this file
|
|
max@0
|
7 // and/or modify it under the terms of the GNU
|
|
max@0
|
8 // Lesser General Public License (LGPL) as published
|
|
max@0
|
9 // by the Free Software Foundation, either version 3
|
|
max@0
|
10 // of the License or (at your option) any later version.
|
|
max@0
|
11 // (see http://www.opensource.org/licenses for more info)
|
|
max@0
|
12
|
|
max@0
|
13
|
|
max@0
|
14 //! \addtogroup op_diagmat
|
|
max@0
|
15 //! @{
|
|
max@0
|
16
|
|
max@0
|
17
|
|
max@0
|
18
|
|
max@0
|
19 template<typename T1>
|
|
max@0
|
20 inline
|
|
max@0
|
21 void
|
|
max@0
|
22 op_diagmat::apply(Mat<typename T1::elem_type>& out, const Op<T1, op_diagmat>& X)
|
|
max@0
|
23 {
|
|
max@0
|
24 arma_extra_debug_sigprint();
|
|
max@0
|
25
|
|
max@0
|
26 typedef typename T1::elem_type eT;
|
|
max@0
|
27
|
|
max@0
|
28 const unwrap<T1> tmp(X.m);
|
|
max@0
|
29 const Mat<eT>& A = tmp.M;
|
|
max@0
|
30
|
|
max@0
|
31 if(A.is_vec() == true)
|
|
max@0
|
32 {
|
|
max@0
|
33 // generate a diagonal matrix out of a vector
|
|
max@0
|
34
|
|
max@0
|
35 const uword N = A.n_elem;
|
|
max@0
|
36 const eT* A_mem = A.memptr();
|
|
max@0
|
37
|
|
max@0
|
38 if(&out != &A)
|
|
max@0
|
39 {
|
|
max@0
|
40 // no aliasing
|
|
max@0
|
41 out.zeros(N,N);
|
|
max@0
|
42
|
|
max@0
|
43 for(uword i=0; i<N; ++i)
|
|
max@0
|
44 {
|
|
max@0
|
45 out.at(i,i) = A_mem[i];
|
|
max@0
|
46 }
|
|
max@0
|
47 }
|
|
max@0
|
48 else
|
|
max@0
|
49 {
|
|
max@0
|
50 // aliasing
|
|
max@0
|
51
|
|
max@0
|
52 const podarray<eT> tmp(A_mem, N);
|
|
max@0
|
53
|
|
max@0
|
54 const eT* tmp_mem = tmp.memptr();
|
|
max@0
|
55
|
|
max@0
|
56 out.zeros(N,N);
|
|
max@0
|
57
|
|
max@0
|
58 for(uword i=0; i<N; ++i)
|
|
max@0
|
59 {
|
|
max@0
|
60 out.at(i,i) = tmp_mem[i];
|
|
max@0
|
61 }
|
|
max@0
|
62 }
|
|
max@0
|
63 }
|
|
max@0
|
64 else
|
|
max@0
|
65 {
|
|
max@0
|
66 // generate a diagonal matrix out of a matrix
|
|
max@0
|
67
|
|
max@0
|
68 arma_debug_check( (A.is_square() == false), "diagmat(): given matrix is not square" );
|
|
max@0
|
69
|
|
max@0
|
70 const uword N = A.n_rows;
|
|
max@0
|
71
|
|
max@0
|
72 if(&out != &A)
|
|
max@0
|
73 {
|
|
max@0
|
74 // no aliasing
|
|
max@0
|
75
|
|
max@0
|
76 out.zeros(N,N);
|
|
max@0
|
77
|
|
max@0
|
78 for(uword i=0; i<N; ++i)
|
|
max@0
|
79 {
|
|
max@0
|
80 out.at(i,i) = A.at(i,i);
|
|
max@0
|
81 }
|
|
max@0
|
82 }
|
|
max@0
|
83 else
|
|
max@0
|
84 {
|
|
max@0
|
85 // aliasing
|
|
max@0
|
86
|
|
max@0
|
87 for(uword i=0; i<N; ++i)
|
|
max@0
|
88 {
|
|
max@0
|
89 eT* colptr = out.colptr(i);
|
|
max@0
|
90
|
|
max@0
|
91 // clear above the diagonal
|
|
max@0
|
92 arrayops::inplace_set(colptr, eT(0), i);
|
|
max@0
|
93
|
|
max@0
|
94 // clear below the diagonal
|
|
max@0
|
95 arrayops::inplace_set(colptr+(i+1), eT(0), N-1-i);
|
|
max@0
|
96 }
|
|
max@0
|
97 }
|
|
max@0
|
98 }
|
|
max@0
|
99 }
|
|
max@0
|
100
|
|
max@0
|
101
|
|
max@0
|
102
|
|
max@0
|
103 //! @}
|
