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_minus
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_minus> & X)28 spglue_minus::apply(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_minus>& 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_minus::apply_noalias(out, pa, pb);
42 }
43 else
44 {
45 SpMat<eT> tmp;
46
47 spglue_minus::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_minus::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(), "subtraction");
64
65 if(pa.get_n_nonzero() == 0) { out = pb.Q; out *= eT(-1); 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_row = x_it.row();
87 const uword x_it_col = x_it.col();
88
89 const uword y_it_row = y_it.row();
90 const uword y_it_col = y_it.col();
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); // take the negative
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_minus::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_minus::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_minus::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_minus_mixed> & expr)185 spglue_minus_mixed::apply(SpMat<typename eT_promoter<T1,T2>::eT>& out, const mtSpGlue<typename eT_promoter<T1,T2>::eT, T1, T2, spglue_minus_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
sparse_minus_dense(Mat<typename promote_type<typename T1::elem_type,typename T2::elem_type>::result> & out,const T1 & X,const T2 & Y)258 spglue_minus_mixed::sparse_minus_dense(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 const quasi_unwrap<T2> UB(Y);
270 const Mat<eT2>& B = UB.M;
271
272 const uword B_n_elem = B.n_elem;
273 const eT2* B_mem = B.memptr();
274
275 out.set_size(B.n_rows, B.n_cols);
276
277 out_eT* out_mem = out.memptr();
278
279 for(uword i=0; i<B_n_elem; ++i)
280 {
281 out_mem[i] = out_eT(-B_mem[i]);
282 }
283
284 const SpProxy<T1> pa(X);
285
286 arma_debug_assert_same_size( pa.get_n_rows(), pa.get_n_cols(), out.n_rows, out.n_cols, "subtraction" );
287
288 typename SpProxy<T1>::const_iterator_type it = pa.begin();
289 typename SpProxy<T1>::const_iterator_type it_end = pa.end();
290
291 while(it != it_end)
292 {
293 out.at(it.row(), it.col()) += out_eT(*it);
294 ++it;
295 }
296 }
297
298
299
300 template<typename T1, typename T2>
301 inline
302 void
dense_minus_sparse(Mat<typename promote_type<typename T1::elem_type,typename T2::elem_type>::result> & out,const T1 & X,const T2 & Y)303 spglue_minus_mixed::dense_minus_sparse(Mat< typename promote_type<typename T1::elem_type, typename T2::elem_type >::result>& out, const T1& X, const T2& Y)
304 {
305 arma_extra_debug_sigprint();
306
307 typedef typename T1::elem_type eT1;
308 typedef typename T2::elem_type eT2;
309
310 typedef typename promote_type<eT1,eT2>::result out_eT;
311
312 promote_type<eT1,eT2>::check();
313
314 if(is_same_type<eT1,out_eT>::no)
315 {
316 out = conv_to< Mat<out_eT> >::from(X);
317 }
318 else
319 {
320 const quasi_unwrap<T1> UA(X);
321
322 const Mat<eT1>& A = UA.M;
323
324 out = reinterpret_cast< const Mat<out_eT>& >(A);
325 }
326
327 const SpProxy<T2> pb(Y);
328
329 arma_debug_assert_same_size( out.n_rows, out.n_cols, pb.get_n_rows(), pb.get_n_cols(), "subtraction" );
330
331 typename SpProxy<T2>::const_iterator_type it = pb.begin();
332 typename SpProxy<T2>::const_iterator_type it_end = pb.end();
333
334 while(it != it_end)
335 {
336 out.at(it.row(), it.col()) -= out_eT(*it);
337 ++it;
338 }
339 }
340
341
342
343 //! @}
344