1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
|
// Copyright (C) 2010-2013 Conrad Sanderson
// Copyright (C) 2010-2013 NICTA (www.nicta.com.au)
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
//! \addtogroup fn_as_scalar
//! @{
template<uword N>
struct as_scalar_redirect
{
template<typename T1>
inline static typename T1::elem_type apply(const T1& X);
};
template<>
struct as_scalar_redirect<2>
{
template<typename T1, typename T2>
inline static typename T1::elem_type apply(const Glue<T1,T2,glue_times>& X);
};
template<>
struct as_scalar_redirect<3>
{
template<typename T1, typename T2, typename T3>
inline static typename T1::elem_type apply(const Glue< Glue<T1, T2, glue_times>, T3, glue_times>& X);
};
template<uword N>
template<typename T1>
inline
typename T1::elem_type
as_scalar_redirect<N>::apply(const T1& X)
{
arma_extra_debug_sigprint();
// typedef typename T1::elem_type eT;
const Proxy<T1> P(X);
arma_debug_check( (P.get_n_elem() != 1), "as_scalar(): expression doesn't evaluate to exactly one element" );
return (Proxy<T1>::prefer_at_accessor == true) ? P.at(0,0) : P[0];
}
template<typename T1, typename T2>
inline
typename T1::elem_type
as_scalar_redirect<2>::apply(const Glue<T1, T2, glue_times>& X)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
// T1 must result in a matrix with one row
// T2 must result in a matrix with one column
const bool has_all_mat = (is_Mat<T1>::value || is_Mat_trans<T1>::value) && (is_Mat<T2>::value || is_Mat_trans<T2>::value);
const bool prefer_at_accessor = Proxy<T1>::prefer_at_accessor || Proxy<T2>::prefer_at_accessor;
const bool do_partial_unwrap = has_all_mat || prefer_at_accessor;
if(do_partial_unwrap == true)
{
const partial_unwrap<T1> tmp1(X.A);
const partial_unwrap<T2> tmp2(X.B);
typedef typename partial_unwrap<T1>::stored_type TA;
typedef typename partial_unwrap<T2>::stored_type TB;
const TA& A = tmp1.M;
const TB& B = tmp2.M;
const uword A_n_rows = (tmp1.do_trans == false) ? (TA::is_row ? 1 : A.n_rows) : (TA::is_col ? 1 : A.n_cols);
const uword A_n_cols = (tmp1.do_trans == false) ? (TA::is_col ? 1 : A.n_cols) : (TA::is_row ? 1 : A.n_rows);
const uword B_n_rows = (tmp2.do_trans == false) ? (TB::is_row ? 1 : B.n_rows) : (TB::is_col ? 1 : B.n_cols);
const uword B_n_cols = (tmp2.do_trans == false) ? (TB::is_col ? 1 : B.n_cols) : (TB::is_row ? 1 : B.n_rows);
arma_debug_check( (A_n_rows != 1) || (B_n_cols != 1) || (A_n_cols != B_n_rows), "as_scalar(): incompatible dimensions" );
const eT val = op_dot::direct_dot(A.n_elem, A.memptr(), B.memptr());
return (tmp1.do_times || tmp2.do_times) ? (val * tmp1.get_val() * tmp2.get_val()) : val;
}
else
{
const Proxy<T1> PA(X.A);
const Proxy<T2> PB(X.B);
arma_debug_check
(
(PA.get_n_rows() != 1) || (PB.get_n_cols() != 1) || (PA.get_n_cols() != PB.get_n_rows()),
"as_scalar(): incompatible dimensions"
);
return op_dot::apply_proxy(PA,PB);
}
}
template<typename T1, typename T2, typename T3>
inline
typename T1::elem_type
as_scalar_redirect<3>::apply(const Glue< Glue<T1, T2, glue_times>, T3, glue_times >& X)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
// T1 * T2 must result in a matrix with one row
// T3 must result in a matrix with one column
typedef typename strip_inv <T2 >::stored_type T2_stripped_1;
typedef typename strip_diagmat<T2_stripped_1>::stored_type T2_stripped_2;
const strip_inv <T2> strip1(X.A.B);
const strip_diagmat<T2_stripped_1> strip2(strip1.M);
const bool tmp2_do_inv = strip1.do_inv;
const bool tmp2_do_diagmat = strip2.do_diagmat;
if(tmp2_do_diagmat == false)
{
const Mat<eT> tmp(X);
arma_debug_check( (tmp.n_elem != 1), "as_scalar(): expression doesn't evaluate to exactly one element" );
return tmp[0];
}
else
{
const partial_unwrap<T1> tmp1(X.A.A);
const partial_unwrap<T2_stripped_2> tmp2(strip2.M);
const partial_unwrap<T3> tmp3(X.B);
const Mat<eT>& A = tmp1.M;
const Mat<eT>& B = tmp2.M;
const Mat<eT>& C = tmp3.M;
const uword A_n_rows = (tmp1.do_trans == false) ? A.n_rows : A.n_cols;
const uword A_n_cols = (tmp1.do_trans == false) ? A.n_cols : A.n_rows;
const bool B_is_vec = B.is_vec();
const uword B_n_rows = (B_is_vec == true) ? B.n_elem : ( (tmp2.do_trans == false) ? B.n_rows : B.n_cols );
const uword B_n_cols = (B_is_vec == true) ? B.n_elem : ( (tmp2.do_trans == false) ? B.n_cols : B.n_rows );
const uword C_n_rows = (tmp3.do_trans == false) ? C.n_rows : C.n_cols;
const uword C_n_cols = (tmp3.do_trans == false) ? C.n_cols : C.n_rows;
const eT val = tmp1.get_val() * tmp2.get_val() * tmp3.get_val();
arma_debug_check
(
(A_n_rows != 1) ||
(C_n_cols != 1) ||
(A_n_cols != B_n_rows) ||
(B_n_cols != C_n_rows)
,
"as_scalar(): incompatible dimensions"
);
if(B_is_vec == true)
{
if(tmp2_do_inv == true)
{
return val * op_dotext::direct_rowvec_invdiagvec_colvec(A.mem, B, C.mem);
}
else
{
return val * op_dot::direct_dot(A.n_elem, A.mem, B.mem, C.mem);
}
}
else
{
if(tmp2_do_inv == true)
{
return val * op_dotext::direct_rowvec_invdiagmat_colvec(A.mem, B, C.mem);
}
else
{
return val * op_dotext::direct_rowvec_diagmat_colvec(A.mem, B, C.mem);
}
}
}
}
template<typename T1>
inline
typename T1::elem_type
as_scalar_diag(const Base<typename T1::elem_type,T1>& X)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
const unwrap<T1> tmp(X.get_ref());
const Mat<eT>& A = tmp.M;
arma_debug_check( (A.n_elem != 1), "as_scalar(): expression doesn't evaluate to exactly one element" );
return A.mem[0];
}
template<typename T1, typename T2, typename T3>
inline
typename T1::elem_type
as_scalar_diag(const Glue< Glue<T1, T2, glue_times_diag>, T3, glue_times >& X)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
// T1 * T2 must result in a matrix with one row
// T3 must result in a matrix with one column
typedef typename strip_diagmat<T2>::stored_type T2_stripped;
const strip_diagmat<T2> strip(X.A.B);
const partial_unwrap<T1> tmp1(X.A.A);
const partial_unwrap<T2_stripped> tmp2(strip.M);
const partial_unwrap<T3> tmp3(X.B);
const Mat<eT>& A = tmp1.M;
const Mat<eT>& B = tmp2.M;
const Mat<eT>& C = tmp3.M;
const uword A_n_rows = (tmp1.do_trans == false) ? A.n_rows : A.n_cols;
const uword A_n_cols = (tmp1.do_trans == false) ? A.n_cols : A.n_rows;
const bool B_is_vec = B.is_vec();
const uword B_n_rows = (B_is_vec == true) ? B.n_elem : ( (tmp2.do_trans == false) ? B.n_rows : B.n_cols );
const uword B_n_cols = (B_is_vec == true) ? B.n_elem : ( (tmp2.do_trans == false) ? B.n_cols : B.n_rows );
const uword C_n_rows = (tmp3.do_trans == false) ? C.n_rows : C.n_cols;
const uword C_n_cols = (tmp3.do_trans == false) ? C.n_cols : C.n_rows;
const eT val = tmp1.get_val() * tmp2.get_val() * tmp3.get_val();
arma_debug_check
(
(A_n_rows != 1) ||
(C_n_cols != 1) ||
(A_n_cols != B_n_rows) ||
(B_n_cols != C_n_rows)
,
"as_scalar(): incompatible dimensions"
);
if(B_is_vec == true)
{
return val * op_dot::direct_dot(A.n_elem, A.mem, B.mem, C.mem);
}
else
{
return val * op_dotext::direct_rowvec_diagmat_colvec(A.mem, B, C.mem);
}
}
template<typename T1, typename T2>
arma_inline
arma_warn_unused
typename T1::elem_type
as_scalar(const Glue<T1, T2, glue_times>& X, const typename arma_not_cx<typename T1::elem_type>::result* junk = 0)
{
arma_extra_debug_sigprint();
arma_ignore(junk);
if(is_glue_times_diag<T1>::value == false)
{
const sword N_mat = 1 + depth_lhs< glue_times, Glue<T1,T2,glue_times> >::num;
arma_extra_debug_print(arma_boost::format("N_mat = %d") % N_mat);
return as_scalar_redirect<N_mat>::apply(X);
}
else
{
return as_scalar_diag(X);
}
}
template<typename T1>
inline
arma_warn_unused
typename T1::elem_type
as_scalar(const Base<typename T1::elem_type,T1>& X)
{
arma_extra_debug_sigprint();
// typedef typename T1::elem_type eT;
const Proxy<T1> P(X.get_ref());
arma_debug_check( (P.get_n_elem() != 1), "as_scalar(): expression doesn't evaluate to exactly one element" );
return (Proxy<T1>::prefer_at_accessor == true) ? P.at(0,0) : P[0];
}
// ensure the following two functions are aware of each other
template<typename T1, typename eop_type> inline arma_warn_unused typename T1::elem_type as_scalar(const eOp<T1, eop_type>& X);
template<typename T1, typename T2, typename eglue_type> inline arma_warn_unused typename T1::elem_type as_scalar(const eGlue<T1, T2, eglue_type>& X);
template<typename T1, typename eop_type>
inline
arma_warn_unused
typename T1::elem_type
as_scalar(const eOp<T1, eop_type>& X)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
const eT val = as_scalar(X.P.Q);
return eop_core<eop_type>::process(val, X.aux);
}
template<typename T1, typename T2, typename eglue_type>
inline
arma_warn_unused
typename T1::elem_type
as_scalar(const eGlue<T1, T2, eglue_type>& X)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
const eT a = as_scalar(X.P1.Q);
const eT b = as_scalar(X.P2.Q);
// the optimiser will keep only one return statement
if(is_same_type<eglue_type, eglue_plus >::yes) { return a + b; }
else if(is_same_type<eglue_type, eglue_minus>::yes) { return a - b; }
else if(is_same_type<eglue_type, eglue_div >::yes) { return a / b; }
else if(is_same_type<eglue_type, eglue_schur>::yes) { return a * b; }
}
template<typename T1>
inline
arma_warn_unused
typename T1::elem_type
as_scalar(const BaseCube<typename T1::elem_type,T1>& X)
{
arma_extra_debug_sigprint();
// typedef typename T1::elem_type eT;
const ProxyCube<T1> P(X.get_ref());
arma_debug_check( (P.get_n_elem() != 1), "as_scalar(): expression doesn't evaluate to exactly one element" );
return (ProxyCube<T1>::prefer_at_accessor == true) ? P.at(0,0,0) : P[0];
}
template<typename T>
arma_inline
arma_warn_unused
const typename arma_scalar_only<T>::result &
as_scalar(const T& x)
{
return x;
}
template<typename T1>
inline
arma_warn_unused
typename T1::elem_type
as_scalar(const SpBase<typename T1::elem_type, T1>& X)
{
typedef typename T1::elem_type eT;
const unwrap_spmat<T1> tmp(X.get_ref());
const SpMat<eT>& A = tmp.M;
arma_debug_check( (A.n_elem != 1), "as_scalar(): expression doesn't evaluate to exactly one element" );
return A.at(0,0);
}
//! @}
|
