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 op_sum
20 //! @{
21
22
23
24 template<typename T1>
25 arma_hot
26 inline
27 void
apply(Mat<typename T1::elem_type> & out,const Op<T1,op_sum> & in)28 op_sum::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_sum>& in)
29 {
30 arma_extra_debug_sigprint();
31
32 typedef typename T1::elem_type eT;
33
34 const uword dim = in.aux_uword_a;
35 arma_debug_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" );
36
37 const Proxy<T1> P(in.m);
38
39 if(P.is_alias(out) == false)
40 {
41 op_sum::apply_noalias(out, P, dim);
42 }
43 else
44 {
45 Mat<eT> tmp;
46
47 op_sum::apply_noalias(tmp, P, dim);
48
49 out.steal_mem(tmp);
50 }
51 }
52
53
54
55 template<typename T1>
56 arma_hot
57 inline
58 void
apply_noalias(Mat<typename T1::elem_type> & out,const Proxy<T1> & P,const uword dim)59 op_sum::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword dim)
60 {
61 arma_extra_debug_sigprint();
62
63 if(is_Mat<typename Proxy<T1>::stored_type>::value || (arma_config::openmp && Proxy<T1>::use_mp))
64 {
65 op_sum::apply_noalias_unwrap(out, P, dim);
66 }
67 else
68 {
69 op_sum::apply_noalias_proxy(out, P, dim);
70 }
71 }
72
73
74
75 template<typename T1>
76 arma_hot
77 inline
78 void
apply_noalias_unwrap(Mat<typename T1::elem_type> & out,const Proxy<T1> & P,const uword dim)79 op_sum::apply_noalias_unwrap(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword dim)
80 {
81 arma_extra_debug_sigprint();
82
83 typedef typename T1::elem_type eT;
84
85 typedef typename Proxy<T1>::stored_type P_stored_type;
86
87 const unwrap<P_stored_type> tmp(P.Q);
88
89 const typename unwrap<P_stored_type>::stored_type& X = tmp.M;
90
91 const uword X_n_rows = X.n_rows;
92 const uword X_n_cols = X.n_cols;
93
94 if(dim == 0)
95 {
96 out.set_size(1, X_n_cols);
97
98 eT* out_mem = out.memptr();
99
100 for(uword col=0; col < X_n_cols; ++col)
101 {
102 out_mem[col] = arrayops::accumulate( X.colptr(col), X_n_rows );
103 }
104 }
105 else
106 {
107 out.zeros(X_n_rows, 1);
108
109 eT* out_mem = out.memptr();
110
111 for(uword col=0; col < X_n_cols; ++col)
112 {
113 arrayops::inplace_plus( out_mem, X.colptr(col), X_n_rows );
114 }
115 }
116 }
117
118
119
120 template<typename T1>
121 arma_hot
122 inline
123 void
apply_noalias_proxy(Mat<typename T1::elem_type> & out,const Proxy<T1> & P,const uword dim)124 op_sum::apply_noalias_proxy(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword dim)
125 {
126 arma_extra_debug_sigprint();
127
128 typedef typename T1::elem_type eT;
129
130 const uword P_n_rows = P.get_n_rows();
131 const uword P_n_cols = P.get_n_cols();
132
133 if(dim == 0)
134 {
135 out.set_size(1, P_n_cols);
136
137 eT* out_mem = out.memptr();
138
139 for(uword col=0; col < P_n_cols; ++col)
140 {
141 eT val1 = eT(0);
142 eT val2 = eT(0);
143
144 uword i,j;
145 for(i=0, j=1; j < P_n_rows; i+=2, j+=2)
146 {
147 val1 += P.at(i,col);
148 val2 += P.at(j,col);
149 }
150
151 if(i < P_n_rows)
152 {
153 val1 += P.at(i,col);
154 }
155
156 out_mem[col] = (val1 + val2);
157 }
158 }
159 else
160 {
161 out.zeros(P_n_rows, 1);
162
163 eT* out_mem = out.memptr();
164
165 for(uword col=0; col < P_n_cols; ++col)
166 for(uword row=0; row < P_n_rows; ++row)
167 {
168 out_mem[row] += P.at(row,col);
169 }
170 }
171 }
172
173
174
175 //
176 // cubes
177
178
179
180 template<typename T1>
181 arma_hot
182 inline
183 void
apply(Cube<typename T1::elem_type> & out,const OpCube<T1,op_sum> & in)184 op_sum::apply(Cube<typename T1::elem_type>& out, const OpCube<T1,op_sum>& in)
185 {
186 arma_extra_debug_sigprint();
187
188 typedef typename T1::elem_type eT;
189
190 const uword dim = in.aux_uword_a;
191 arma_debug_check( (dim > 2), "sum(): parameter 'dim' must be 0 or 1 or 2" );
192
193 const ProxyCube<T1> P(in.m);
194
195 if(P.is_alias(out) == false)
196 {
197 op_sum::apply_noalias(out, P, dim);
198 }
199 else
200 {
201 Cube<eT> tmp;
202
203 op_sum::apply_noalias(tmp, P, dim);
204
205 out.steal_mem(tmp);
206 }
207 }
208
209
210
211 template<typename T1>
212 arma_hot
213 inline
214 void
apply_noalias(Cube<typename T1::elem_type> & out,const ProxyCube<T1> & P,const uword dim)215 op_sum::apply_noalias(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim)
216 {
217 arma_extra_debug_sigprint();
218
219 if(is_Cube<typename ProxyCube<T1>::stored_type>::value || (arma_config::openmp && ProxyCube<T1>::use_mp))
220 {
221 op_sum::apply_noalias_unwrap(out, P, dim);
222 }
223 else
224 {
225 op_sum::apply_noalias_proxy(out, P, dim);
226 }
227 }
228
229
230
231 template<typename T1>
232 arma_hot
233 inline
234 void
apply_noalias_unwrap(Cube<typename T1::elem_type> & out,const ProxyCube<T1> & P,const uword dim)235 op_sum::apply_noalias_unwrap(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim)
236 {
237 arma_extra_debug_sigprint();
238
239 typedef typename T1::elem_type eT;
240
241 typedef typename ProxyCube<T1>::stored_type P_stored_type;
242
243 const unwrap_cube<P_stored_type> tmp(P.Q);
244
245 const Cube<eT>& X = tmp.M;
246
247 const uword X_n_rows = X.n_rows;
248 const uword X_n_cols = X.n_cols;
249 const uword X_n_slices = X.n_slices;
250
251 if(dim == 0)
252 {
253 out.set_size(1, X_n_cols, X_n_slices);
254
255 for(uword slice=0; slice < X_n_slices; ++slice)
256 {
257 eT* out_mem = out.slice_memptr(slice);
258
259 for(uword col=0; col < X_n_cols; ++col)
260 {
261 out_mem[col] = arrayops::accumulate( X.slice_colptr(slice,col), X_n_rows );
262 }
263 }
264 }
265 else
266 if(dim == 1)
267 {
268 out.zeros(X_n_rows, 1, X_n_slices);
269
270 for(uword slice=0; slice < X_n_slices; ++slice)
271 {
272 eT* out_mem = out.slice_memptr(slice);
273
274 for(uword col=0; col < X_n_cols; ++col)
275 {
276 arrayops::inplace_plus( out_mem, X.slice_colptr(slice,col), X_n_rows );
277 }
278 }
279 }
280 else
281 if(dim == 2)
282 {
283 out.zeros(X_n_rows, X_n_cols, 1);
284
285 eT* out_mem = out.memptr();
286
287 for(uword slice=0; slice < X_n_slices; ++slice)
288 {
289 arrayops::inplace_plus(out_mem, X.slice_memptr(slice), X.n_elem_slice );
290 }
291 }
292 }
293
294
295
296 template<typename T1>
297 arma_hot
298 inline
299 void
apply_noalias_proxy(Cube<typename T1::elem_type> & out,const ProxyCube<T1> & P,const uword dim)300 op_sum::apply_noalias_proxy(Cube<typename T1::elem_type>& out, const ProxyCube<T1>& P, const uword dim)
301 {
302 arma_extra_debug_sigprint();
303
304 typedef typename T1::elem_type eT;
305
306 const uword P_n_rows = P.get_n_rows();
307 const uword P_n_cols = P.get_n_cols();
308 const uword P_n_slices = P.get_n_slices();
309
310 if(dim == 0)
311 {
312 out.set_size(1, P_n_cols, P_n_slices);
313
314 for(uword slice=0; slice < P_n_slices; ++slice)
315 {
316 eT* out_mem = out.slice_memptr(slice);
317
318 for(uword col=0; col < P_n_cols; ++col)
319 {
320 eT val1 = eT(0);
321 eT val2 = eT(0);
322
323 uword i,j;
324 for(i=0, j=1; j < P_n_rows; i+=2, j+=2)
325 {
326 val1 += P.at(i,col,slice);
327 val2 += P.at(j,col,slice);
328 }
329
330 if(i < P_n_rows)
331 {
332 val1 += P.at(i,col,slice);
333 }
334
335 out_mem[col] = (val1 + val2);
336 }
337 }
338 }
339 else
340 if(dim == 1)
341 {
342 out.zeros(P_n_rows, 1, P_n_slices);
343
344 for(uword slice=0; slice < P_n_slices; ++slice)
345 {
346 eT* out_mem = out.slice_memptr(slice);
347
348 for(uword col=0; col < P_n_cols; ++col)
349 for(uword row=0; row < P_n_rows; ++row)
350 {
351 out_mem[row] += P.at(row,col,slice);
352 }
353 }
354 }
355 else
356 if(dim == 2)
357 {
358 out.zeros(P_n_rows, P_n_cols, 1);
359
360 for(uword slice=0; slice < P_n_slices; ++slice)
361 {
362 for(uword col=0; col < P_n_cols; ++col)
363 {
364 eT* out_mem = out.slice_colptr(0,col);
365
366 for(uword row=0; row < P_n_rows; ++row)
367 {
368 out_mem[row] += P.at(row,col,slice);
369 }
370 }
371 }
372 }
373 }
374
375
376
377 //! @}
378