Chris@49
|
1 // Copyright (C) 2011-2012 NICTA (www.nicta.com.au)
|
Chris@49
|
2 // Copyright (C) 2011-2012 Conrad Sanderson
|
Chris@49
|
3 //
|
Chris@49
|
4 // This Source Code Form is subject to the terms of the Mozilla Public
|
Chris@49
|
5 // License, v. 2.0. If a copy of the MPL was not distributed with this
|
Chris@49
|
6 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
Chris@49
|
7
|
Chris@49
|
8
|
Chris@49
|
9 //! \addtogroup op_symmat
|
Chris@49
|
10 //! @{
|
Chris@49
|
11
|
Chris@49
|
12
|
Chris@49
|
13
|
Chris@49
|
14 template<typename T1>
|
Chris@49
|
15 inline
|
Chris@49
|
16 void
|
Chris@49
|
17 op_symmat::apply
|
Chris@49
|
18 (
|
Chris@49
|
19 Mat<typename T1::elem_type>& out,
|
Chris@49
|
20 const Op<T1,op_symmat>& in,
|
Chris@49
|
21 const typename arma_not_cx<typename T1::elem_type>::result* junk
|
Chris@49
|
22 )
|
Chris@49
|
23 {
|
Chris@49
|
24 arma_extra_debug_sigprint();
|
Chris@49
|
25 arma_ignore(junk);
|
Chris@49
|
26
|
Chris@49
|
27 typedef typename T1::elem_type eT;
|
Chris@49
|
28
|
Chris@49
|
29 const unwrap<T1> tmp(in.m);
|
Chris@49
|
30 const Mat<eT>& A = tmp.M;
|
Chris@49
|
31
|
Chris@49
|
32 arma_debug_check( (A.is_square() == false), "symmatu()/symmatl(): given matrix must be square" );
|
Chris@49
|
33
|
Chris@49
|
34 const uword N = A.n_rows;
|
Chris@49
|
35 const bool upper = (in.aux_uword_a == 0);
|
Chris@49
|
36
|
Chris@49
|
37 if(&out != &A)
|
Chris@49
|
38 {
|
Chris@49
|
39 out.copy_size(A);
|
Chris@49
|
40
|
Chris@49
|
41 if(upper)
|
Chris@49
|
42 {
|
Chris@49
|
43 // upper triangular: copy the diagonal and the elements above the diagonal
|
Chris@49
|
44
|
Chris@49
|
45 for(uword i=0; i<N; ++i)
|
Chris@49
|
46 {
|
Chris@49
|
47 const eT* A_data = A.colptr(i);
|
Chris@49
|
48 eT* out_data = out.colptr(i);
|
Chris@49
|
49
|
Chris@49
|
50 arrayops::copy( out_data, A_data, i+1 );
|
Chris@49
|
51 }
|
Chris@49
|
52 }
|
Chris@49
|
53 else
|
Chris@49
|
54 {
|
Chris@49
|
55 // lower triangular: copy the diagonal and the elements below the diagonal
|
Chris@49
|
56
|
Chris@49
|
57 for(uword i=0; i<N; ++i)
|
Chris@49
|
58 {
|
Chris@49
|
59 const eT* A_data = A.colptr(i);
|
Chris@49
|
60 eT* out_data = out.colptr(i);
|
Chris@49
|
61
|
Chris@49
|
62 arrayops::copy( &out_data[i], &A_data[i], N-i );
|
Chris@49
|
63 }
|
Chris@49
|
64 }
|
Chris@49
|
65 }
|
Chris@49
|
66
|
Chris@49
|
67
|
Chris@49
|
68 if(upper)
|
Chris@49
|
69 {
|
Chris@49
|
70 // reflect elements across the diagonal from upper triangle to lower triangle
|
Chris@49
|
71
|
Chris@49
|
72 for(uword col=1; col < N; ++col)
|
Chris@49
|
73 {
|
Chris@49
|
74 const eT* coldata = out.colptr(col);
|
Chris@49
|
75
|
Chris@49
|
76 for(uword row=0; row < col; ++row)
|
Chris@49
|
77 {
|
Chris@49
|
78 out.at(col,row) = coldata[row];
|
Chris@49
|
79 }
|
Chris@49
|
80 }
|
Chris@49
|
81 }
|
Chris@49
|
82 else
|
Chris@49
|
83 {
|
Chris@49
|
84 // reflect elements across the diagonal from lower triangle to upper triangle
|
Chris@49
|
85
|
Chris@49
|
86 for(uword col=0; col < N; ++col)
|
Chris@49
|
87 {
|
Chris@49
|
88 const eT* coldata = out.colptr(col);
|
Chris@49
|
89
|
Chris@49
|
90 for(uword row=(col+1); row < N; ++row)
|
Chris@49
|
91 {
|
Chris@49
|
92 out.at(col,row) = coldata[row];
|
Chris@49
|
93 }
|
Chris@49
|
94 }
|
Chris@49
|
95 }
|
Chris@49
|
96 }
|
Chris@49
|
97
|
Chris@49
|
98
|
Chris@49
|
99
|
Chris@49
|
100 template<typename T1>
|
Chris@49
|
101 inline
|
Chris@49
|
102 void
|
Chris@49
|
103 op_symmat::apply
|
Chris@49
|
104 (
|
Chris@49
|
105 Mat<typename T1::elem_type>& out,
|
Chris@49
|
106 const Op<T1,op_symmat>& in,
|
Chris@49
|
107 const typename arma_cx_only<typename T1::elem_type>::result* junk
|
Chris@49
|
108 )
|
Chris@49
|
109 {
|
Chris@49
|
110 arma_extra_debug_sigprint();
|
Chris@49
|
111 arma_ignore(junk);
|
Chris@49
|
112
|
Chris@49
|
113 typedef typename T1::elem_type eT;
|
Chris@49
|
114
|
Chris@49
|
115 const unwrap<T1> tmp(in.m);
|
Chris@49
|
116 const Mat<eT>& A = tmp.M;
|
Chris@49
|
117
|
Chris@49
|
118 arma_debug_check( (A.is_square() == false), "symmatu()/symmatl(): given matrix must be square" );
|
Chris@49
|
119
|
Chris@49
|
120 const uword N = A.n_rows;
|
Chris@49
|
121 const bool upper = (in.aux_uword_a == 0);
|
Chris@49
|
122
|
Chris@49
|
123 if(&out != &A)
|
Chris@49
|
124 {
|
Chris@49
|
125 out.copy_size(A);
|
Chris@49
|
126
|
Chris@49
|
127 if(upper)
|
Chris@49
|
128 {
|
Chris@49
|
129 // upper triangular: copy the diagonal and the elements above the diagonal
|
Chris@49
|
130
|
Chris@49
|
131 for(uword i=0; i<N; ++i)
|
Chris@49
|
132 {
|
Chris@49
|
133 const eT* A_data = A.colptr(i);
|
Chris@49
|
134 eT* out_data = out.colptr(i);
|
Chris@49
|
135
|
Chris@49
|
136 arrayops::copy( out_data, A_data, i+1 );
|
Chris@49
|
137 }
|
Chris@49
|
138 }
|
Chris@49
|
139 else
|
Chris@49
|
140 {
|
Chris@49
|
141 // lower triangular: copy the diagonal and the elements below the diagonal
|
Chris@49
|
142
|
Chris@49
|
143 for(uword i=0; i<N; ++i)
|
Chris@49
|
144 {
|
Chris@49
|
145 const eT* A_data = A.colptr(i);
|
Chris@49
|
146 eT* out_data = out.colptr(i);
|
Chris@49
|
147
|
Chris@49
|
148 arrayops::copy( &out_data[i], &A_data[i], N-i );
|
Chris@49
|
149 }
|
Chris@49
|
150 }
|
Chris@49
|
151 }
|
Chris@49
|
152
|
Chris@49
|
153
|
Chris@49
|
154 if(upper)
|
Chris@49
|
155 {
|
Chris@49
|
156 // reflect elements across the diagonal from upper triangle to lower triangle
|
Chris@49
|
157
|
Chris@49
|
158 for(uword col=1; col < N; ++col)
|
Chris@49
|
159 {
|
Chris@49
|
160 const eT* coldata = out.colptr(col);
|
Chris@49
|
161
|
Chris@49
|
162 for(uword row=0; row < col; ++row)
|
Chris@49
|
163 {
|
Chris@49
|
164 out.at(col,row) = std::conj(coldata[row]);
|
Chris@49
|
165 }
|
Chris@49
|
166 }
|
Chris@49
|
167 }
|
Chris@49
|
168 else
|
Chris@49
|
169 {
|
Chris@49
|
170 // reflect elements across the diagonal from lower triangle to upper triangle
|
Chris@49
|
171
|
Chris@49
|
172 for(uword col=0; col < N; ++col)
|
Chris@49
|
173 {
|
Chris@49
|
174 const eT* coldata = out.colptr(col);
|
Chris@49
|
175
|
Chris@49
|
176 for(uword row=(col+1); row < N; ++row)
|
Chris@49
|
177 {
|
Chris@49
|
178 out.at(col,row) = std::conj(coldata[row]);
|
Chris@49
|
179 }
|
Chris@49
|
180 }
|
Chris@49
|
181 }
|
Chris@49
|
182 }
|
Chris@49
|
183
|
Chris@49
|
184
|
Chris@49
|
185
|
Chris@49
|
186 //! @}
|