Chris@49: // Copyright (C) 2012 Ryan Curtin
Chris@49: // Copyright (C) 2012 Conrad Sanderson
Chris@49: // 
Chris@49: // This Source Code Form is subject to the terms of the Mozilla Public
Chris@49: // License, v. 2.0. If a copy of the MPL was not distributed with this
Chris@49: // file, You can obtain one at http://mozilla.org/MPL/2.0/.
Chris@49: 
Chris@49: 
Chris@49: //! \addtogroup spglue_minus
Chris@49: //! @{
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1, typename T2>
Chris@49: inline
Chris@49: void
Chris@49: spglue_minus::apply(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_minus>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   
Chris@49:   const SpProxy<T1> pa(X.A);
Chris@49:   const SpProxy<T2> pb(X.B);
Chris@49:   
Chris@49:   const bool is_alias = pa.is_alias(out) || pb.is_alias(out);
Chris@49:   
Chris@49:   if(is_alias == false)
Chris@49:     {
Chris@49:     spglue_minus::apply_noalias(out, pa, pb);
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     SpMat<eT> tmp;
Chris@49:     spglue_minus::apply_noalias(tmp, pa, pb);
Chris@49:     
Chris@49:     out.steal_mem(tmp);
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT, typename T1, typename T2>
Chris@49: arma_hot
Chris@49: inline
Chris@49: void
Chris@49: spglue_minus::apply_noalias(SpMat<eT>& result, const SpProxy<T1>& pa, const SpProxy<T2>& pb)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "subtraction");
Chris@49:   
Chris@49:   if( (pa.get_n_nonzero() != 0) && (pb.get_n_nonzero() != 0) )
Chris@49:     {
Chris@49:     result.set_size(pa.get_n_rows(), pa.get_n_cols());
Chris@49: 
Chris@49:     // Resize memory to correct size.
Chris@49:     result.mem_resize(n_unique(pa, pb, op_n_unique_sub()));
Chris@49: 
Chris@49:     // Now iterate across both matrices.
Chris@49:     typename SpProxy<T1>::const_iterator_type x_it = pa.begin();
Chris@49:     typename SpProxy<T2>::const_iterator_type y_it = pb.begin();
Chris@49: 
Chris@49:     typename SpProxy<T1>::const_iterator_type x_end = pa.end();
Chris@49:     typename SpProxy<T2>::const_iterator_type y_end = pb.end();
Chris@49:     
Chris@49:     uword cur_val = 0;
Chris@49:     while((x_it != x_end) || (y_it != y_end))
Chris@49:       {
Chris@49:       if(x_it == y_it)
Chris@49:         {
Chris@49:         const eT val = (*x_it) - (*y_it);
Chris@49:         
Chris@49:         if (val != eT(0))
Chris@49:           {
Chris@49:           access::rw(result.values[cur_val]) = val;
Chris@49:           access::rw(result.row_indices[cur_val]) = x_it.row();
Chris@49:           ++access::rw(result.col_ptrs[x_it.col() + 1]);
Chris@49:           ++cur_val;
Chris@49:           }
Chris@49: 
Chris@49:         ++x_it;
Chris@49:         ++y_it;
Chris@49:         }
Chris@49:       else
Chris@49:         {
Chris@49:         const uword x_it_row = x_it.row();
Chris@49:         const uword x_it_col = x_it.col();
Chris@49:         
Chris@49:         const uword y_it_row = y_it.row();
Chris@49:         const uword y_it_col = y_it.col();
Chris@49:         
Chris@49:         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
Chris@49:           {
Chris@49:           access::rw(result.values[cur_val]) = (*x_it);
Chris@49:           access::rw(result.row_indices[cur_val]) = x_it_row;
Chris@49:           ++access::rw(result.col_ptrs[x_it_col + 1]);
Chris@49:           ++cur_val;
Chris@49:           ++x_it;
Chris@49:           }
Chris@49:         else
Chris@49:           {
Chris@49:           access::rw(result.values[cur_val]) = -(*y_it);
Chris@49:           access::rw(result.row_indices[cur_val]) = y_it_row;
Chris@49:           ++access::rw(result.col_ptrs[y_it_col + 1]);
Chris@49:           ++cur_val;
Chris@49:           ++y_it;
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49: 
Chris@49:     // Fix column pointers to be cumulative.
Chris@49:     for(uword c = 1; c <= result.n_cols; ++c)
Chris@49:       {
Chris@49:       access::rw(result.col_ptrs[c]) += result.col_ptrs[c - 1];
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     if(pa.get_n_nonzero() == 0)
Chris@49:       {
Chris@49:       result = pb.Q;
Chris@49:       result *= eT(-1);
Chris@49:       
Chris@49:       return;
Chris@49:       }
Chris@49:     
Chris@49:     if(pb.get_n_nonzero() == 0)
Chris@49:       {
Chris@49:       result = pa.Q;
Chris@49:       return;
Chris@49:       }
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //
Chris@49: //
Chris@49: // spglue_minus2: scalar*(A - B)
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1, typename T2>
Chris@49: inline
Chris@49: void
Chris@49: spglue_minus2::apply(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_minus2>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   
Chris@49:   const SpProxy<T1> pa(X.A);
Chris@49:   const SpProxy<T2> pb(X.B);
Chris@49:   
Chris@49:   const bool is_alias = pa.is_alias(out) || pb.is_alias(out);
Chris@49:   
Chris@49:   if(is_alias == false)
Chris@49:     {
Chris@49:     spglue_minus::apply_noalias(out, pa, pb);
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     SpMat<eT> tmp;
Chris@49:     spglue_minus::apply_noalias(tmp, pa, pb);
Chris@49:     
Chris@49:     out.steal_mem(tmp);
Chris@49:     }
Chris@49:   
Chris@49:   out *= X.aux;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}