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 spop_diagmat
20 //! @{
21
22
23
24 template<typename T1>
25 inline
26 void
apply(SpMat<typename T1::elem_type> & out,const SpOp<T1,spop_diagmat> & in)27 spop_diagmat::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1, spop_diagmat>& in)
28 {
29 arma_extra_debug_sigprint();
30
31 typedef typename T1::elem_type eT;
32
33 if(in.is_alias(out) == false)
34 {
35 spop_diagmat::apply_noalias(out, in.m);
36 }
37 else
38 {
39 SpMat<eT> tmp;
40
41 spop_diagmat::apply_noalias(tmp, in.m);
42
43 out.steal_mem(tmp);
44 }
45 }
46
47
48
49 template<typename T1>
50 inline
51 void
apply_noalias(SpMat<typename T1::elem_type> & out,const SpBase<typename T1::elem_type,T1> & expr)52 spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpBase<typename T1::elem_type, T1>& expr)
53 {
54 arma_extra_debug_sigprint();
55
56 typedef typename T1::elem_type eT;
57
58 const SpProxy<T1> P(expr.get_ref());
59
60 const uword P_n_rows = P.get_n_rows();
61 const uword P_n_cols = P.get_n_cols();
62 const uword P_n_nz = P.get_n_nonzero();
63
64 const bool P_is_vec = (P_n_rows == 1) || (P_n_cols == 1);
65
66 if(P_is_vec) // generate a diagonal matrix out of a vector
67 {
68 const uword N = (P_n_rows == 1) ? P_n_cols : P_n_rows;
69
70 out.zeros(N, N);
71
72 if(P_n_nz == 0) { return; }
73
74 typename SpProxy<T1>::const_iterator_type it = P.begin();
75
76 if(P_n_cols == 1)
77 {
78 for(uword i=0; i < P_n_nz; ++i)
79 {
80 const uword row = it.row();
81
82 out.at(row,row) = (*it);
83
84 ++it;
85 }
86 }
87 else
88 if(P_n_rows == 1)
89 {
90 for(uword i=0; i < P_n_nz; ++i)
91 {
92 const uword col = it.col();
93
94 out.at(col,col) = (*it);
95
96 ++it;
97 }
98 }
99 }
100 else // generate a diagonal matrix out of a matrix
101 {
102 out.zeros(P_n_rows, P_n_cols);
103
104 const uword N = (std::min)(P_n_rows, P_n_cols);
105
106 if( (is_SpMat<typename SpProxy<T1>::stored_type>::value) && (P_n_nz >= 5*N) )
107 {
108 const unwrap_spmat<typename SpProxy<T1>::stored_type> U(P.Q);
109
110 const SpMat<eT>& X = U.M;
111
112 for(uword i=0; i < N; ++i)
113 {
114 const eT val = X.at(i,i); // use binary search
115
116 if(val != eT(0)) { out.at(i,i) = val; }
117 }
118 }
119 else
120 {
121 if(P_n_nz == 0) { return; }
122
123 typename SpProxy<T1>::const_iterator_type it = P.begin();
124
125 for(uword i=0; i < P_n_nz; ++i)
126 {
127 const uword row = it.row();
128 const uword col = it.col();
129
130 if(row == col) { out.at(row,row) = (*it); }
131
132 ++it;
133 }
134 }
135 }
136 }
137
138
139
140 template<typename T1, typename T2>
141 inline
142 void
apply_noalias(SpMat<typename T1::elem_type> & out,const SpGlue<T1,T2,spglue_plus> & expr)143 spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_plus>& expr)
144 {
145 arma_extra_debug_sigprint();
146
147 typedef typename T1::elem_type eT;
148
149 const unwrap_spmat<T1> UA(expr.A);
150 const unwrap_spmat<T2> UB(expr.B);
151
152 const SpMat<eT>& A = UA.M;
153 const SpMat<eT>& B = UB.M;
154
155 arma_debug_assert_same_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "addition");
156
157 const bool is_vec = (A.n_rows == 1) || (A.n_cols == 1);
158
159 if(is_vec) // generate a diagonal matrix out of a vector
160 {
161 const uword N = (A.n_rows == 1) ? A.n_cols : A.n_rows;
162
163 out.zeros(N,N);
164
165 if(A.n_rows == 1)
166 {
167 for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(0,i) + B.at(0,i); }
168 }
169 else
170 {
171 for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(i,0) + B.at(i,0); }
172 }
173 }
174 else // generate a diagonal matrix out of a matrix
175 {
176 SpMat<eT> AA; spop_diagmat::apply_noalias(AA, A);
177 SpMat<eT> BB; spop_diagmat::apply_noalias(BB, B);
178
179 out = AA + BB;
180 }
181 }
182
183
184
185 template<typename T1, typename T2>
186 inline
187 void
apply_noalias(SpMat<typename T1::elem_type> & out,const SpGlue<T1,T2,spglue_minus> & expr)188 spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_minus>& expr)
189 {
190 arma_extra_debug_sigprint();
191
192 typedef typename T1::elem_type eT;
193
194 const unwrap_spmat<T1> UA(expr.A);
195 const unwrap_spmat<T2> UB(expr.B);
196
197 const SpMat<eT>& A = UA.M;
198 const SpMat<eT>& B = UB.M;
199
200 arma_debug_assert_same_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "subtraction");
201
202 const bool is_vec = (A.n_rows == 1) || (A.n_cols == 1);
203
204 if(is_vec) // generate a diagonal matrix out of a vector
205 {
206 const uword N = (A.n_rows == 1) ? A.n_cols : A.n_rows;
207
208 out.zeros(N,N);
209
210 if(A.n_rows == 1)
211 {
212 for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(0,i) - B.at(0,i); }
213 }
214 else
215 {
216 for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(i,0) - B.at(i,0); }
217 }
218 }
219 else // generate a diagonal matrix out of a matrix
220 {
221 SpMat<eT> AA; spop_diagmat::apply_noalias(AA, A);
222 SpMat<eT> BB; spop_diagmat::apply_noalias(BB, B);
223
224 out = AA - BB;
225 }
226 }
227
228
229
230 template<typename T1, typename T2>
231 inline
232 void
apply_noalias(SpMat<typename T1::elem_type> & out,const SpGlue<T1,T2,spglue_schur> & expr)233 spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_schur>& expr)
234 {
235 arma_extra_debug_sigprint();
236
237 typedef typename T1::elem_type eT;
238
239 const unwrap_spmat<T1> UA(expr.A);
240 const unwrap_spmat<T2> UB(expr.B);
241
242 const SpMat<eT>& A = UA.M;
243 const SpMat<eT>& B = UB.M;
244
245 arma_debug_assert_same_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "element-wise multiplication");
246
247 const bool is_vec = (A.n_rows == 1) || (A.n_cols == 1);
248
249 if(is_vec) // generate a diagonal matrix out of a vector
250 {
251 const uword N = (A.n_rows == 1) ? A.n_cols : A.n_rows;
252
253 out.zeros(N,N);
254
255 if(A.n_rows == 1)
256 {
257 for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(0,i) * B.at(0,i); }
258 }
259 else
260 {
261 for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(i,0) * B.at(i,0); }
262 }
263 }
264 else // generate a diagonal matrix out of a matrix
265 {
266 SpMat<eT> AA; spop_diagmat::apply_noalias(AA, A);
267 SpMat<eT> BB; spop_diagmat::apply_noalias(BB, B);
268
269 out = AA % BB;
270 }
271 }
272
273
274
275 template<typename T1, typename T2>
276 inline
277 void
apply_noalias(SpMat<typename T1::elem_type> & out,const SpGlue<T1,T2,spglue_times> & expr)278 spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_times>& expr)
279 {
280 arma_extra_debug_sigprint();
281
282 typedef typename T1::elem_type eT;
283
284 const unwrap_spmat<T1> UA(expr.A);
285 const unwrap_spmat<T2> UB(expr.B);
286
287 const SpMat<eT>& A = UA.M;
288 const SpMat<eT>& B = UB.M;
289
290 arma_debug_assert_mul_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "matrix multiplication");
291
292 const uword C_n_rows = A.n_rows;
293 const uword C_n_cols = B.n_cols;
294
295 const bool is_vec = (C_n_rows == 1) || (C_n_cols == 1);
296
297 if(is_vec) // generate a diagonal matrix out of a vector
298 {
299 const SpMat<eT> C = A*B;
300
301 spop_diagmat::apply_noalias(out, C);
302 }
303 else // generate a diagonal matrix out of a matrix
304 {
305 const uword N = (std::min)(C_n_rows, C_n_cols);
306
307 if( (A.n_nonzero >= 5*N) || (B.n_nonzero >= 5*N) )
308 {
309 out.zeros(C_n_rows, C_n_cols);
310
311 for(uword k=0; k < N; ++k)
312 {
313 typename SpMat<eT>::const_col_iterator B_it = B.begin_col_no_sync(k);
314 typename SpMat<eT>::const_col_iterator B_it_end = B.end_col_no_sync(k);
315
316 eT acc = eT(0);
317
318 while(B_it != B_it_end)
319 {
320 const eT B_val = (*B_it);
321 const uword i = B_it.row();
322
323 acc += A.at(k,i) * B_val;
324
325 ++B_it;
326 }
327
328 out(k,k) = acc;
329 }
330 }
331 else
332 {
333 const SpMat<eT> C = A*B;
334
335 spop_diagmat::apply_noalias(out, C);
336 }
337 }
338 }
339
340
341
342 //
343 //
344
345
346
347 template<typename T1>
348 inline
349 void
apply(SpMat<typename T1::elem_type> & out,const SpOp<T1,spop_diagmat2> & in)350 spop_diagmat2::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1, spop_diagmat2>& in)
351 {
352 arma_extra_debug_sigprint();
353
354 typedef typename T1::elem_type eT;
355
356 const uword row_offset = in.aux_uword_a;
357 const uword col_offset = in.aux_uword_b;
358
359 const unwrap_spmat<T1> U(in.m);
360
361 if(U.is_alias(out))
362 {
363 SpMat<eT> tmp;
364
365 spop_diagmat2::apply_noalias(tmp, U.M, row_offset, col_offset);
366
367 out.steal_mem(tmp);
368 }
369 else
370 {
371 spop_diagmat2::apply_noalias(out, U.M, row_offset, col_offset);
372 }
373 }
374
375
376
377 template<typename eT>
378 inline
379 void
apply_noalias(SpMat<eT> & out,const SpMat<eT> & X,const uword row_offset,const uword col_offset)380 spop_diagmat2::apply_noalias(SpMat<eT>& out, const SpMat<eT>& X, const uword row_offset, const uword col_offset)
381 {
382 arma_extra_debug_sigprint();
383
384 const uword n_rows = X.n_rows;
385 const uword n_cols = X.n_cols;
386 const uword n_elem = X.n_elem;
387
388 if(n_elem == 0) { out.reset(); return; }
389
390 const bool X_is_vec = (n_rows == 1) || (n_cols == 1);
391
392 if(X_is_vec) // generate a diagonal matrix out of a vector
393 {
394 const uword n_pad = (std::max)(row_offset, col_offset);
395
396 out.zeros(n_elem + n_pad, n_elem + n_pad);
397
398 const uword X_n_nz = X.n_nonzero;
399
400 if(X_n_nz == 0) { return; }
401
402 typename SpMat<eT>::const_iterator it = X.begin();
403
404 if(n_cols == 1)
405 {
406 for(uword i=0; i < X_n_nz; ++i)
407 {
408 const uword row = it.row();
409
410 out.at(row_offset + row, col_offset + row) = (*it);
411
412 ++it;
413 }
414 }
415 else
416 if(n_rows == 1)
417 {
418 for(uword i=0; i < X_n_nz; ++i)
419 {
420 const uword col = it.col();
421
422 out.at(row_offset + col, col_offset + col) = (*it);
423
424 ++it;
425 }
426 }
427 }
428 else // generate a diagonal matrix out of a matrix
429 {
430 arma_debug_check_bounds
431 (
432 ((row_offset > 0) && (row_offset >= n_rows)) || ((col_offset > 0) && (col_offset >= n_cols)),
433 "diagmat(): requested diagonal out of bounds"
434 );
435
436 out.zeros(n_rows, n_cols);
437
438 if(X.n_nonzero == 0) { return; }
439
440 const uword N = (std::min)(n_rows - row_offset, n_cols - col_offset);
441
442 for(uword i=0; i<N; ++i)
443 {
444 const uword row = i + row_offset;
445 const uword col = i + col_offset;
446
447 const eT val = X.at(row,col);
448
449 if(val != eT(0)) { out.at(row,col) = val; }
450 }
451 }
452 }
453
454
455
456 //! @}
457