C++ SpProxy::get_n_nonzero方法代码示例

inline
arma_warn_unused
typename
enable_if2
  <(is_arma_sparse_type<T1>::value) && (is_arma_sparse_type<T2>::value) && (is_same_type<typename T1::elem_type, typename T2::elem_type>::value),
   typename T1::elem_type
  >::result
dot
  (
  const SpBase<typename T1::elem_type, T1>& x,
  const SpBase<typename T2::elem_type, T2>& y
  )
  {
  arma_extra_debug_sigprint();

  const SpProxy<T1> pa(x.get_ref());
  const SpProxy<T2> pb(y.get_ref());

  arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "dot()");

  typedef typename T1::elem_type eT;

  if((&(x.get_ref()) == &(y.get_ref())) && (SpProxy<T1>::must_use_iterator == false))
    {
    // We can do it directly!
    return op_dot::direct_dot_arma(pa.get_n_nonzero(), pa.get_values(), pa.get_values());
    }
  else
    {
    // Iterate over both objects and see when they are the same
    eT result = eT(0);

    typename SpProxy<T1>::const_iterator_type a_it = pa.begin();
    typename SpProxy<T2>::const_iterator_type b_it = pb.begin();

    while((a_it.pos() < pa.get_n_nonzero()) && (b_it.pos() < pb.get_n_nonzero()))
      {
      if(a_it == b_it)
        {
        result += (*a_it) * (*b_it);

        ++a_it;
        ++b_it;
        }
      else if((a_it.col() < b_it.col()) || ((a_it.col() == b_it.col()) && (a_it.row() < b_it.row())))
        {
        // a_it is "behind"
        ++a_it;
        }
      else
        {
        // b_it is "behind"
        ++b_it;
        }
      }

    return result;
    }
  }

inline
void
op_nonzeros::apply_noalias(Mat<typename T1::elem_type>& out, const SpBase<typename T1::elem_type,T1>& X)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
  const SpProxy<T1> P(X.get_ref());
  
  const uword N = P.get_n_nonzero();
  
  out.set_size(N,1);
  
  if(N > 0)
    {
    if(is_SpMat<typename SpProxy<T1>::stored_type>::value)
      {
      const unwrap_spmat<typename SpProxy<T1>::stored_type> U(P.Q);
      
      arrayops::copy(out.memptr(), U.M.values, N);
      }
    else
      {
      eT* out_mem = out.memptr();
      
      typename SpProxy<T1>::const_iterator_type it = P.begin();
      
      for(uword i=0; i<N; ++i)  { out_mem[i] = (*it); ++it; }
      }
    }
  }

inline
arma_warn_unused
typename
enable_if2
  <(is_arma_type<T1>::value) && (is_arma_sparse_type<T2>::value) && (is_same_type<typename T1::elem_type, typename T2::elem_type>::value),
   typename T1::elem_type
  >::result
dot
  (
  const Base<typename T1::elem_type, T1>& x,
  const SpBase<typename T2::elem_type, T2>& y
  )
  {
  arma_extra_debug_sigprint();

  const Proxy<T1> pa(x.get_ref());
  const SpProxy<T2> pb(y.get_ref());

  arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "dot()");

  typedef typename T1::elem_type eT;

  eT result = eT(0);

  typename SpProxy<T2>::const_iterator_type it = pb.begin();

  // prefer_at_accessor won't save us operations
  while(it.pos() < pb.get_n_nonzero())
    {
    result += (*it) * pa.at(it.row(), it.col());
    ++it;
    }

  return result;
  }

arma_warn_unused
inline
Col<uword>
find(const SpBase<typename T1::elem_type,T1>& X, const uword k = 0)
  {
  arma_extra_debug_sigprint();
  
  const SpProxy<T1> P(X.get_ref());
  
  const uword n_rows = P.get_n_rows();
  const uword n_nz   = P.get_n_nonzero();
  
  Mat<uword> tmp(n_nz,1);
  
  uword* tmp_mem = tmp.memptr();
  
  typename SpProxy<T1>::const_iterator_type it = P.begin();
  
  for(uword i=0; i<n_nz; ++i)
    {
    const uword index = it.row() + it.col()*n_rows;
    
    tmp_mem[i] = index;
    
    ++it;
    }
  
  Col<uword> out;
  
  const uword count = (k == 0) ? uword(n_nz) : uword( (std::min)(n_nz, k) );
  
  out.steal_mem_col(tmp, count);
  
  return out;
  }

arma_hot
inline
void
spop_htrans::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1,spop_htrans>& in, const typename arma_cx_only<typename T1::elem_type>::result* junk)
{
    arma_extra_debug_sigprint();
    arma_ignore(junk);

    typedef typename   T1::elem_type  eT;
    typedef typename umat::elem_type ueT;

    const SpProxy<T1> p(in.m);

    const uword N = p.get_n_nonzero();

    if(N == uword(0))
    {
        out.set_size(p.get_n_cols(), p.get_n_rows());
        return;
    }

    umat locs(2, N);

    Col<eT> vals(N);

    eT* vals_ptr = vals.memptr();

    typename SpProxy<T1>::const_iterator_type it = p.begin();

    for(uword count = 0; count < N; ++count)
    {
        ueT* locs_ptr = locs.colptr(count);

        locs_ptr[0] = it.col();
        locs_ptr[1] = it.row();

        vals_ptr[count] = std::conj(*it);

        ++it;
    }

    SpMat<eT> tmp(locs, vals, p.get_n_cols(), p.get_n_rows());

    out.steal_mem(tmp);
}

inline
typename
enable_if2
  <
  (is_arma_type<T1>::value && is_arma_sparse_type<T2>::value && is_same_type<typename T1::elem_type, typename T2::elem_type>::value),
  Mat<typename T1::elem_type>
  >::result
operator/
  (
  const   Base<typename T1::elem_type, T1>& x,
  const SpBase<typename T2::elem_type, T2>& y
  )
  {
  arma_extra_debug_sigprint();

  const   Proxy<T1> pa(x.get_ref());
  const SpProxy<T2> pb(y.get_ref());

  arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "element-wise division");

  Mat<typename T1::elem_type> result(pa.get_n_rows(), pa.get_n_cols());

  result.fill(Datum<typename T1::elem_type>::inf);

  // Now divide each element
  typename SpProxy<T2>::const_iterator_type it = pb.begin();

  while(it.pos() < pb.get_n_nonzero())
    {
    if(Proxy<T1>::prefer_at_accessor == false)
      {
      const uword index = (it.col() * result.n_rows) + it.row();
      result[index] = pa[index] / (*it);
      }
    else
      {
      result.at(it.row(), it.col()) = pa.at(it.row(), it.col()) / (*it);
      }

    ++it;
    }

  return result;
  }

inline
void
spop_mean::apply_noalias_fast
  (
        SpMat<typename T1::elem_type>& out,
  const SpProxy<T1>&                   p,
  const uword                          dim
  )
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  typedef typename T1::pod_type   T;
  
  const uword p_n_rows = p.get_n_rows();
  const uword p_n_cols = p.get_n_cols();
  
  if( (p_n_rows == 0) || (p_n_cols == 0) || (p.get_n_nonzero() == 0) )
    {
    if(dim == 0)  { out.zeros((p_n_rows > 0) ? 1 : 0, p_n_cols); }
    if(dim == 1)  { out.zeros(p_n_rows, (p_n_cols > 0) ? 1 : 0); }
    
    return;
    }
  
  if(dim == 0) // find the mean in each column
    {
    Row<eT> acc(p_n_cols, fill::zeros);
    
    if(SpProxy<T1>::must_use_iterator)
      {
      typename SpProxy<T1>::const_iterator_type it     = p.begin();
      typename SpProxy<T1>::const_iterator_type it_end = p.end();
      
      while(it != it_end)  { acc[it.col()] += (*it);  ++it; }
      
      acc /= T(p_n_rows);
      }
    else
      {
      for(uword col = 0; col < p_n_cols; ++col)
        {
        acc[col] = arrayops::accumulate
          (
          &p.get_values()[p.get_col_ptrs()[col]],
          p.get_col_ptrs()[col + 1] - p.get_col_ptrs()[col]
          ) / T(p_n_rows);
        }
      }
    
    out = acc;
    }
  else
  if(dim == 1)  // find the mean in each row
    {
    Col<eT> acc(p_n_rows, fill::zeros);
    
    typename SpProxy<T1>::const_iterator_type it     = p.begin();
    typename SpProxy<T1>::const_iterator_type it_end = p.end();
    
    while(it != it_end)  { acc[it.row()] += (*it);  ++it; }
    
    acc /= T(p_n_cols);
    
    out = acc;
    }
  
  if(out.is_finite() == false)
    {
    spop_mean::apply_noalias_slow(out, p, dim);
    }
  }

inline
void
spop_mean::apply_noalias_slow
  (
        SpMat<typename T1::elem_type>& out,
  const SpProxy<T1>&                   p,
  const uword                          dim
  )
  {
  arma_extra_debug_sigprint();

  typedef typename T1::elem_type eT;
  
  const uword p_n_rows = p.get_n_rows();
  const uword p_n_cols = p.get_n_cols();

  if(dim == 0)  // find the mean in each column
    {
    arma_extra_debug_print("spop_mean::apply_noalias(): dim = 0");
    
    out.set_size((p_n_rows > 0) ? 1 : 0, p_n_cols);
    
    if( (p_n_rows == 0) || (p.get_n_nonzero() == 0) )  { return; }
    
    for(uword col = 0; col < p_n_cols; ++col)
      {
      // Do we have to use an iterator or can we use memory directly?
      if(SpProxy<T1>::must_use_iterator)
        {
        typename SpProxy<T1>::const_iterator_type it  = p.begin_col(col);
        typename SpProxy<T1>::const_iterator_type end = p.begin_col(col + 1);
        
        const uword n_zero = p_n_rows - (end.pos() - it.pos());
        
        out.at(0,col) = spop_mean::iterator_mean(it, end, n_zero, eT(0));
        }
      else
        {
        out.at(0,col) = spop_mean::direct_mean
          (
          &p.get_values()[p.get_col_ptrs()[col]],
          p.get_col_ptrs()[col + 1] - p.get_col_ptrs()[col],
          p_n_rows
          );
        }
      }
    }
  else
  if(dim == 1)  // find the mean in each row
    {
    arma_extra_debug_print("spop_mean::apply_noalias(): dim = 1");
    
    out.set_size(p_n_rows, (p_n_cols > 0) ? 1 : 0);
    
    if( (p_n_cols == 0) || (p.get_n_nonzero() == 0) )  { return; }
    
    for(uword row = 0; row < p_n_rows; ++row)
      {
      // We must use an iterator regardless of how it is stored.
      typename SpProxy<T1>::const_row_iterator_type it  = p.begin_row(row);
      typename SpProxy<T1>::const_row_iterator_type end = p.end_row(row);
      
      const uword n_zero = p_n_cols - (end.pos() - it.pos());
      
      out.at(row,0) = spop_mean::iterator_mean(it, end, n_zero, eT(0));
      }
    }
  }

inline
typename
enable_if2
  <
  (is_arma_sparse_type<T1>::value && is_arma_sparse_type<T2>::value && is_same_type<typename T1::elem_type, typename T2::elem_type>::value),
  SpMat<typename T1::elem_type>
  >::result
operator%
  (
  const SpBase<typename T1::elem_type, T1>& x,
  const SpBase<typename T2::elem_type, T2>& y
  )
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;

  const SpProxy<T1> pa(x.get_ref());
  const SpProxy<T2> pb(y.get_ref());

  arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "element-wise multiplication");

  SpMat<typename T1::elem_type> result(pa.get_n_rows(), pa.get_n_cols());
  
  if( (pa.get_n_nonzero() != 0) && (pb.get_n_nonzero() != 0) )
    {
    // Resize memory to correct size.
    result.mem_resize(n_unique(x, y, op_n_unique_mul()));
    
    // Now iterate across both matrices.
    typename SpProxy<T1>::const_iterator_type x_it = pa.begin();
    typename SpProxy<T2>::const_iterator_type y_it = pb.begin();
    
    typename SpProxy<T1>::const_iterator_type x_end = pa.end();
    typename SpProxy<T2>::const_iterator_type y_end = pb.end();
    
    uword cur_val = 0;
    while((x_it != x_end) || (y_it != y_end))
      {
      if(x_it == y_it)
        {
        const eT val = (*x_it) * (*y_it);
        
        if (val != eT(0))
          {
          access::rw(result.values[cur_val]) = val;
          access::rw(result.row_indices[cur_val]) = x_it.row();
          ++access::rw(result.col_ptrs[x_it.col() + 1]);
          ++cur_val;
          }
        
        ++x_it;
        ++y_it;
        }
      else
        {
        const uword x_it_row = x_it.row();
        const uword x_it_col = x_it.col();
        
        const uword y_it_row = y_it.row();
        const uword y_it_col = y_it.col();
        
        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
          {
          ++x_it;
          }
        else
          {
          ++y_it;
          }
        }
      }
    
    // Fix column pointers to be cumulative.
    for(uword c = 1; c <= result.n_cols; ++c)
      {
      access::rw(result.col_ptrs[c]) += result.col_ptrs[c - 1];
      }
    }
  
  return result;
  }

inline
void
spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpProxy<T1>& p)
  {
  arma_extra_debug_sigprint();
  
  const uword n_rows = p.get_n_rows();
  const uword n_cols = p.get_n_cols();
  
  const bool p_is_vec = (n_rows == 1) || (n_cols == 1);
  
  if(p_is_vec)    // generate a diagonal matrix out of a vector
    {
    const uword N = (n_rows == 1) ? n_cols : n_rows;
    
    out.zeros(N, N);
    
    if(p.get_n_nonzero() == 0)  { return; }
    
    typename SpProxy<T1>::const_iterator_type it     = p.begin();
    typename SpProxy<T1>::const_iterator_type it_end = p.end();
      
    if(n_cols == 1)
      {
      while(it != it_end)
        {
        const uword row = it.row();
        
        out.at(row,row) = (*it);
        
        ++it;
        }
      }
    else
    if(n_rows == 1)
      {
      while(it != it_end)
        {
        const uword col = it.col();
        
        out.at(col,col) = (*it);
        
        ++it;
        }
      }
    }
  else   // generate a diagonal matrix out of a matrix
    {
    arma_debug_check( (n_rows != n_cols), "diagmat(): given matrix is not square" );
    
    out.zeros(n_rows, n_rows);
    
    if(p.get_n_nonzero() == 0)  { return; }
    
    typename SpProxy<T1>::const_iterator_type it     = p.begin();
    typename SpProxy<T1>::const_iterator_type it_end = p.end();
      
    while(it != it_end)
      {
      const uword row = it.row();
      const uword col = it.col();
      
      if(row == col)
        {
        out.at(row,row) = (*it);
        }
      
      ++it;
      }
    }
  }

inline
const SpSubview<eT>&
SpSubview<eT>::operator_equ_common(const SpBase<eT, T1>& in)
  {
  arma_extra_debug_sigprint();
  
  // algorithm:
  // instead of directly inserting values into the matrix underlying the subview,
  // create a new matrix by merging the underlying matrix with the input object,
  // and then replacing the underlying matrix with the created matrix.
  // 
  // the merging process requires pretending that the input object
  // has the same size as the underlying matrix.
  // while iterating through the elements of the input object,
  // this requires adjusting the row and column locations of each element,
  // as well as providing fake zero elements.
  // in effect there is a proxy for a proxy.
  
  
  const SpProxy< SpMat<eT> > pa((*this).m   );
  const SpProxy< T1        > pb(in.get_ref());
  
  arma_debug_assert_same_size(n_rows, n_cols, pb.get_n_rows(), pb.get_n_cols(), "insertion into sparse submatrix");
  
  const uword pa_start_row = (*this).aux_row1;
  const uword pa_start_col = (*this).aux_col1;
  
  const uword pa_end_row = pa_start_row + (*this).n_rows - 1;
  const uword pa_end_col = pa_start_col + (*this).n_cols - 1;
  
  const uword pa_n_rows = pa.get_n_rows();
  
  SpMat<eT> out(pa.get_n_rows(), pa.get_n_cols());
  
  const uword alt_count = pa.get_n_nonzero() - (*this).n_nonzero + pb.get_n_nonzero();
  
  // Resize memory to correct size.
  out.mem_resize(alt_count);
  
  typename SpProxy< SpMat<eT> >::const_iterator_type x_it  = pa.begin();
  typename SpProxy< SpMat<eT> >::const_iterator_type x_end = pa.end();
  
  typename SpProxy<T1>::const_iterator_type y_it  = pb.begin();
  typename SpProxy<T1>::const_iterator_type y_end = pb.end();
  
  bool x_it_ok = (x_it != x_end);
  bool y_it_ok = (y_it != y_end);
  
  uword x_it_row = (x_it_ok) ? x_it.row() : 0;
  uword x_it_col = (x_it_ok) ? x_it.col() : 0;
  
  uword y_it_row = (y_it_ok) ? y_it.row() + pa_start_row : 0;
  uword y_it_col = (y_it_ok) ? y_it.col() + pa_start_col : 0;
    
  uword cur_val = 0;
  while(x_it_ok || y_it_ok)
    {
    const bool x_inside_box = (x_it_row >= pa_start_row) && (x_it_row <= pa_end_row) && (x_it_col >= pa_start_col) && (x_it_col <= pa_end_col);
    const bool y_inside_box = (y_it_row >= pa_start_row) && (y_it_row <= pa_end_row) && (y_it_col >= pa_start_col) && (y_it_col <= pa_end_col);
    
    const eT x_val = x_inside_box ? eT(0) : ( x_it_ok ? (*x_it) : eT(0) );
    
    const eT y_val = y_inside_box ? ( y_it_ok ? (*y_it) : eT(0) ) : eT(0);
    
    if( (x_it_row == y_it_row) && (x_it_col == y_it_col) )
      {
      if( (x_val != eT(0)) || (y_val != eT(0)) )  
        {
        access::rw(out.values[cur_val]) = (x_val != eT(0)) ? x_val : y_val;
        access::rw(out.row_indices[cur_val]) = x_it_row;
        ++access::rw(out.col_ptrs[x_it_col + 1]);
        ++cur_val;
        }
      
      if(x_it_ok)
        {
        ++x_it;
        
        if(x_it == x_end)  { x_it_ok = false; }
        }
      
      if(x_it_ok)
        {
        x_it_row = x_it.row();
        x_it_col = x_it.col();
        }
      else
        {
        x_it_row++;
        
        if(x_it_row >= pa_n_rows)  { x_it_row = 0; x_it_col++; }
        }
      
      if(y_it_ok)
        {
        ++y_it;
        
        if(y_it == y_end)  { y_it_ok = false; }
        }
      
//.........这里部分代码省略.........

arma_hot
inline
void
spglue_minus::apply_noalias(SpMat<eT>& out, const SpProxy<T1>& pa, const SpProxy<T2>& pb)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "subtraction");
  
  if(pa.get_n_nonzero() == 0)  { out = pb.Q; out *= eT(-1); return; }
  if(pb.get_n_nonzero() == 0)  { out = pa.Q;                return; }
  
  const uword max_n_nonzero = spglue_elem_helper::max_n_nonzero_plus(pa, pb);
  
  // Resize memory to upper bound
  out.reserve(pa.get_n_rows(), pa.get_n_cols(), max_n_nonzero);
  
  // Now iterate across both matrices.
  typename SpProxy<T1>::const_iterator_type x_it  = pa.begin();
  typename SpProxy<T1>::const_iterator_type x_end = pa.end();
  
  typename SpProxy<T2>::const_iterator_type y_it  = pb.begin();
  typename SpProxy<T2>::const_iterator_type y_end = pb.end();
  
  uword count = 0;
  
  while( (x_it != x_end) || (y_it != y_end) )
    {
    eT out_val;
    
    const uword x_it_row = x_it.row();
    const uword x_it_col = x_it.col();
    
    const uword y_it_row = y_it.row();
    const uword y_it_col = y_it.col();
    
    bool use_y_loc = false;
    
    if(x_it == y_it)
      {
      out_val = (*x_it) - (*y_it);
      
      ++x_it;
      ++y_it;
      }
    else
      {
      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
        {
        out_val = (*x_it);
        
        ++x_it;
        }
      else
        {
        out_val = -(*y_it);  // take the negative
        
        ++y_it;
        
        use_y_loc = true;
        }
      }
    
    if(out_val != eT(0))
      {
      access::rw(out.values[count]) = out_val;
      
      const uword out_row = (use_y_loc == false) ? x_it_row : y_it_row;
      const uword out_col = (use_y_loc == false) ? x_it_col : y_it_col;
      
      access::rw(out.row_indices[count]) = out_row;
      access::rw(out.col_ptrs[out_col + 1])++;
      ++count;
      }
    }
  
  const uword out_n_cols = out.n_cols;
  
  uword* col_ptrs = access::rwp(out.col_ptrs);
  
  // Fix column pointers to be cumulative.
  for(uword c = 1; c <= out_n_cols; ++c)
    {
    col_ptrs[c] += col_ptrs[c - 1];
    }
  
  if(count < max_n_nonzero)
    {
    if(count <= (max_n_nonzero/2))
      {
      out.mem_resize(count);
      }
    else
      {
      // quick resize without reallocating memory and copying data
      access::rw(         out.n_nonzero) = count;
      access::rw(     out.values[count]) = eT(0);
      access::rw(out.row_indices[count]) = uword(0);
      }
    }
//.........这里部分代码省略.........

inline
typename
enable_if2
<
(is_arma_type<T1>::value && is_arma_sparse_type<T2>::value && is_same_type<typename T1::elem_type, typename T2::elem_type>::value),
SpMat<typename T1::elem_type>
>::result
operator%
(
    const   Base<typename T1::elem_type, T1>& x,
    const SpBase<typename T2::elem_type, T2>& y
)
{
    arma_extra_debug_sigprint();

    const   Proxy<T1> pa(x.get_ref());
    const SpProxy<T2> pb(y.get_ref());

    arma_debug_assert_same_size(pa.get_n_rows(), pa.get_n_cols(), pb.get_n_rows(), pb.get_n_cols(), "element-wise multiplication");

    SpMat<typename T1::elem_type> result(pa.get_n_rows(), pa.get_n_cols());

    if(Proxy<T1>::prefer_at_accessor == false)
    {
        // use direct operator[] access
        // count new size
        uword new_n_nonzero = 0;
        typename SpProxy<T2>::const_iterator_type it = pb.begin();
        while(it.pos() < pb.get_n_nonzero())
        {
            if(((*it) * pa[(it.col() * pa.get_n_rows()) + it.row()]) != 0)
            {
                ++new_n_nonzero;
            }

            ++it;
        }

        // Resize memory accordingly.
        result.mem_resize(new_n_nonzero);

        uword cur_val = 0;
        typename SpProxy<T2>::const_iterator_type it2 = pb.begin();
        while(it2.pos() < pb.get_n_nonzero())
        {
            const typename T1::elem_type val = (*it2) * pa[(it2.col() * pa.get_n_rows()) + it2.row()];
            if(val != 0)
            {
                access::rw(result.values[cur_val]) = val;
                access::rw(result.row_indices[cur_val]) = it2.row();
                ++access::rw(result.col_ptrs[it2.col() + 1]);
                ++cur_val;
            }

            ++it2;
        }
    }
    else
    {
        // use at() access
        // count new size
        uword new_n_nonzero = 0;
        typename SpProxy<T2>::const_iterator_type it = pb.begin();
        while(it.pos() < pb.get_n_nonzero())
        {
            if(((*it) * pa.at(it.row(), it.col())) != 0)
            {
                ++new_n_nonzero;
            }

            ++it;
        }

        // Resize memory accordingly.
        result.mem_resize(new_n_nonzero);

        uword cur_val = 0;
        typename SpProxy<T2>::const_iterator_type it2 = pb.begin();
        while(it2.pos() < pb.get_n_nonzero())
        {
            const typename T1::elem_type val = (*it2) * pa.at(it2.row(), it2.col());
            if(val != 0)
            {
                access::rw(result.values[cur_val]) = val;
                access::rw(result.row_indices[cur_val]) = it2.row();
                ++access::rw(result.col_ptrs[it2.col() + 1]);
                ++cur_val;
            }

            ++it2;
        }
    }

    // Fix column pointers.
    for(uword c = 1; c <= result.n_cols; ++c)
    {
        access::rw(result.col_ptrs[c]) += result.col_ptrs[c - 1];
    }

    return result;
//.........这里部分代码省略.........

arma_hot
inline
void
spop_sum::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1,spop_sum>& in)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
  const uword dim = in.aux_uword_a;
  arma_debug_check( (dim > 1), "sum(): parameter 'dim' must be 0 or 1" );
  
  const SpProxy<T1> p(in.m);
  
  const uword p_n_rows = p.get_n_rows();
  const uword p_n_cols = p.get_n_cols();
  
  if(p.get_n_nonzero() == 0)
    {
    if(dim == 0)  { out.zeros(1,p_n_cols); }
    if(dim == 1)  { out.zeros(p_n_rows,1); }
    
    return;
    }
  
  if(dim == 0) // find the sum in each column
    {
    Row<eT> acc(p_n_cols, fill::zeros);
    
    if(SpProxy<T1>::must_use_iterator)
      {
      typename SpProxy<T1>::const_iterator_type it     = p.begin();
      typename SpProxy<T1>::const_iterator_type it_end = p.end();
      
      while(it != it_end)  { acc[it.col()] += (*it);  ++it; }
      }
    else
      {
      for(uword col = 0; col < p_n_cols; ++col)
        {
        acc[col] = arrayops::accumulate
          (
          &p.get_values()[p.get_col_ptrs()[col]],
          p.get_col_ptrs()[col + 1] - p.get_col_ptrs()[col]
          );
        }
      }
    
    out = acc;
    }
  else
  if(dim == 1)  // find the sum in each row
    {
    Col<eT> acc(p_n_rows, fill::zeros);
    
    typename SpProxy<T1>::const_iterator_type it     = p.begin();
    typename SpProxy<T1>::const_iterator_type it_end = p.end();
    
    while(it != it_end)  { acc[it.row()] += (*it);  ++it; }
    
    out = acc;
    }
  }

inline
void
spop_var::apply_noalias
  (
        SpMat<typename T1::pod_type>& out,
  const SpProxy<T1>&                  p,
  const uword                         norm_type,
  const uword                         dim
  )
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type  in_eT;
  //typedef typename T1::pod_type  out_eT;
  
  const uword p_n_rows = p.get_n_rows();
  const uword p_n_cols = p.get_n_cols();
  
  // TODO: this is slow; rewrite based on the approach used by sparse mean()
  
  if(dim == 0)  // find variance in each column
    {
    arma_extra_debug_print("spop_var::apply_noalias(): dim = 0");
    
    out.set_size((p_n_rows > 0) ? 1 : 0, p_n_cols);
    
    if( (p_n_rows == 0) || (p.get_n_nonzero() == 0) )  { return; }
    
    for(uword col = 0; col < p_n_cols; ++col)
      {
      if(SpProxy<T1>::must_use_iterator)
        {
        // We must use an iterator; we can't access memory directly.
        typename SpProxy<T1>::const_iterator_type it  = p.begin_col(col);
        typename SpProxy<T1>::const_iterator_type end = p.begin_col(col + 1);
        
        const uword n_zero = p_n_rows - (end.pos() - it.pos());
        
        // in_eT is used just to get the specialization right (complex / noncomplex)
        out.at(0, col) = spop_var::iterator_var(it, end, n_zero, norm_type, in_eT(0));
        }
      else
        {
        // We can use direct memory access to calculate the variance.
        out.at(0, col) = spop_var::direct_var
          (
          &p.get_values()[p.get_col_ptrs()[col]],
          p.get_col_ptrs()[col + 1] - p.get_col_ptrs()[col],
          p_n_rows,
          norm_type
          );
        }
      }
    }
  else
  if(dim == 1)  // find variance in each row
    {
    arma_extra_debug_print("spop_var::apply_noalias(): dim = 1");
    
    out.set_size(p_n_rows, (p_n_cols > 0) ? 1 : 0);
    
    if( (p_n_cols == 0) || (p.get_n_nonzero() == 0) )  { return; }
    
    for(uword row = 0; row < p_n_rows; ++row)
      {
      // We have to use an iterator here regardless of whether or not we can
      // directly access memory.
      typename SpProxy<T1>::const_row_iterator_type it  = p.begin_row(row);
      typename SpProxy<T1>::const_row_iterator_type end = p.end_row(row);
      
      const uword n_zero = p_n_cols - (end.pos() - it.pos());
      
      out.at(row, 0) = spop_var::iterator_var(it, end, n_zero, norm_type, in_eT(0));
      }
    }
  }