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 spglue_plus
20 //! @{
21
22
23
24 template<typename T1, typename T2>
25 arma_hot
26 inline
27 void
apply(SpMat<typename T1::elem_type> & out,const SpGlue<T1,T2,spglue_plus> & X)28 spglue_plus::apply(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_plus>& X)
29 {
30 arma_extra_debug_sigprint();
31
32 typedef typename T1::elem_type eT;
33
34 const SpProxy<T1> pa(X.A);
35 const SpProxy<T2> pb(X.B);
36
37 const bool is_alias = pa.is_alias(out) || pb.is_alias(out);
38
39 if(is_alias == false)
40 {
41 spglue_plus::apply_noalias(out, pa, pb);
42 }
43 else
44 {
45 SpMat<eT> tmp;
46
47 spglue_plus::apply_noalias(tmp, pa, pb);
48
49 out.steal_mem(tmp);
50 }
51 }
52
53
54
55 template<typename eT, typename T1, typename T2>
56 arma_hot
57 inline
58 void
apply_noalias(SpMat<eT> & out,const SpProxy<T1> & pa,const SpProxy<T2> & pb)59 spglue_plus::apply_noalias(SpMat<eT>& out, const SpProxy<T1>& pa, const SpProxy<T2>& pb)
60 {
61 arma_extra_debug_sigprint();
62
63 arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "addition");
64
65 if(pa.get_n_nonzero() == 0) { out = pb.Q; return; }
66 if(pb.get_n_nonzero() == 0) { out = pa.Q; return; }
67
68 const uword max_n_nonzero = pa.get_n_nonzero() + pb.get_n_nonzero();
69
70 // Resize memory to upper bound
71 out.reserve(pa.get_n_rows(), pa.get_n_cols(), max_n_nonzero);
72
73 // Now iterate across both matrices.
74 typename SpProxy<T1>::const_iterator_type x_it = pa.begin();
75 typename SpProxy<T1>::const_iterator_type x_end = pa.end();
76
77 typename SpProxy<T2>::const_iterator_type y_it = pb.begin();
78 typename SpProxy<T2>::const_iterator_type y_end = pb.end();
79
80 uword count = 0;
81
82 while( (x_it != x_end) || (y_it != y_end) )
83 {
84 eT out_val;
85
86 const uword x_it_col = x_it.col();
87 const uword x_it_row = x_it.row();
88
89 const uword y_it_col = y_it.col();
90 const uword y_it_row = y_it.row();
91
92 bool use_y_loc = false;
93
94 if(x_it == y_it)
95 {
96 out_val = (*x_it) + (*y_it);
97
98 ++x_it;
99 ++y_it;
100 }
101 else
102 {
103 if((x_it_col < y_it_col) || ((x_it_col == y_it_col) && (x_it_row < y_it_row))) // if y is closer to the end
104 {
105 out_val = (*x_it);
106
107 ++x_it;
108 }
109 else
110 {
111 out_val = (*y_it);
112
113 ++y_it;
114
115 use_y_loc = true;
116 }
117 }
118
119 if(out_val != eT(0))
120 {
121 access::rw(out.values[count]) = out_val;
122
123 const uword out_row = (use_y_loc == false) ? x_it_row : y_it_row;
124 const uword out_col = (use_y_loc == false) ? x_it_col : y_it_col;
125
126 access::rw(out.row_indices[count]) = out_row;
127 access::rw(out.col_ptrs[out_col + 1])++;
128 ++count;
129 }
130
131 arma_check( (count > max_n_nonzero), "internal error: spglue_plus::apply_noalias(): count > max_n_nonzero" );
132 }
133
134 const uword out_n_cols = out.n_cols;
135
136 uword* col_ptrs = access::rwp(out.col_ptrs);
137
138 // Fix column pointers to be cumulative.
139 for(uword c = 1; c <= out_n_cols; ++c)
140 {
141 col_ptrs[c] += col_ptrs[c - 1];
142 }
143
144 if(count < max_n_nonzero)
145 {
146 if(count <= (max_n_nonzero/2))
147 {
148 out.mem_resize(count);
149 }
150 else
151 {
152 // quick resize without reallocating memory and copying data
153 access::rw( out.n_nonzero) = count;
154 access::rw( out.values[count]) = eT(0);
155 access::rw(out.row_indices[count]) = uword(0);
156 }
157 }
158 }
159
160
161
162 template<typename eT>
163 arma_hot
164 inline
165 void
apply_noalias(SpMat<eT> & out,const SpMat<eT> & A,const SpMat<eT> & B)166 spglue_plus::apply_noalias(SpMat<eT>& out, const SpMat<eT>& A, const SpMat<eT>& B)
167 {
168 arma_extra_debug_sigprint();
169
170 const SpProxy< SpMat<eT> > pa(A);
171 const SpProxy< SpMat<eT> > pb(B);
172
173 spglue_plus::apply_noalias(out, pa, pb);
174 }
175
176
177
178 //
179
180
181
182 template<typename T1, typename T2>
183 inline
184 void
apply(SpMat<typename eT_promoter<T1,T2>::eT> & out,const mtSpGlue<typename eT_promoter<T1,T2>::eT,T1,T2,spglue_plus_mixed> & expr)185 spglue_plus_mixed::apply(SpMat<typename eT_promoter<T1,T2>::eT>& out, const mtSpGlue<typename eT_promoter<T1,T2>::eT, T1, T2, spglue_plus_mixed>& expr)
186 {
187 arma_extra_debug_sigprint();
188
189 typedef typename T1::elem_type eT1;
190 typedef typename T2::elem_type eT2;
191
192 typedef typename promote_type<eT1,eT2>::result out_eT;
193
194 promote_type<eT1,eT2>::check();
195
196 if( (is_same_type<eT1,out_eT>::no) && (is_same_type<eT2,out_eT>::yes) )
197 {
198 // upgrade T1
199
200 const unwrap_spmat<T1> UA(expr.A);
201 const unwrap_spmat<T2> UB(expr.B);
202
203 const SpMat<eT1>& A = UA.M;
204 const SpMat<eT2>& B = UB.M;
205
206 SpMat<out_eT> AA(arma_layout_indicator(), A);
207
208 for(uword i=0; i < A.n_nonzero; ++i) { access::rw(AA.values[i]) = out_eT(A.values[i]); }
209
210 const SpMat<out_eT>& BB = reinterpret_cast< const SpMat<out_eT>& >(B);
211
212 out = AA + BB;
213 }
214 else
215 if( (is_same_type<eT1,out_eT>::yes) && (is_same_type<eT2,out_eT>::no) )
216 {
217 // upgrade T2
218
219 const unwrap_spmat<T1> UA(expr.A);
220 const unwrap_spmat<T2> UB(expr.B);
221
222 const SpMat<eT1>& A = UA.M;
223 const SpMat<eT2>& B = UB.M;
224
225 const SpMat<out_eT>& AA = reinterpret_cast< const SpMat<out_eT>& >(A);
226
227 SpMat<out_eT> BB(arma_layout_indicator(), B);
228
229 for(uword i=0; i < B.n_nonzero; ++i) { access::rw(BB.values[i]) = out_eT(B.values[i]); }
230
231 out = AA + BB;
232 }
233 else
234 {
235 // upgrade T1 and T2
236
237 const unwrap_spmat<T1> UA(expr.A);
238 const unwrap_spmat<T2> UB(expr.B);
239
240 const SpMat<eT1>& A = UA.M;
241 const SpMat<eT2>& B = UB.M;
242
243 SpMat<out_eT> AA(arma_layout_indicator(), A);
244 SpMat<out_eT> BB(arma_layout_indicator(), B);
245
246 for(uword i=0; i < A.n_nonzero; ++i) { access::rw(AA.values[i]) = out_eT(A.values[i]); }
247 for(uword i=0; i < B.n_nonzero; ++i) { access::rw(BB.values[i]) = out_eT(B.values[i]); }
248
249 out = AA + BB;
250 }
251 }
252
253
254
255 template<typename T1, typename T2>
256 inline
257 void
dense_plus_sparse(Mat<typename promote_type<typename T1::elem_type,typename T2::elem_type>::result> & out,const T1 & X,const T2 & Y)258 spglue_plus_mixed::dense_plus_sparse(Mat< typename promote_type<typename T1::elem_type, typename T2::elem_type >::result>& out, const T1& X, const T2& Y)
259 {
260 arma_extra_debug_sigprint();
261
262 typedef typename T1::elem_type eT1;
263 typedef typename T2::elem_type eT2;
264
265 typedef typename promote_type<eT1,eT2>::result out_eT;
266
267 promote_type<eT1,eT2>::check();
268
269 if(is_same_type<eT1,out_eT>::no)
270 {
271 out = conv_to< Mat<out_eT> >::from(X);
272 }
273 else
274 {
275 const quasi_unwrap<T1> UA(X);
276
277 const Mat<eT1>& A = UA.M;
278
279 out = reinterpret_cast< const Mat<out_eT>& >(A);
280 }
281
282 const SpProxy<T2> pb(Y);
283
284 arma_debug_assert_same_size( out.n_rows, out.n_cols, pb.get_n_rows(), pb.get_n_cols(), "addition" );
285
286 typename SpProxy<T2>::const_iterator_type it = pb.begin();
287 typename SpProxy<T2>::const_iterator_type it_end = pb.end();
288
289 while(it != it_end)
290 {
291 out.at(it.row(), it.col()) += out_eT(*it);
292 ++it;
293 }
294 }
295
296
297
298 //! @}
299