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