vnl_sparse_matrix.txx

InsightToolkit-1.4.0(有大量的优化算法程序)

//------------------------------------------------------------
//: Add rhs to this.
template <class T>
void vnl_sparse_matrix<T>::add(const vnl_sparse_matrix<T>& rhs,
                               vnl_sparse_matrix<T>& result) const
{
  assert((rhs.rows() == rows()) && (rhs.columns() == columns()));

  // Clear result matrix.
  result.elements.clear();

  // Now give the result matrix enough rows.
  result.elements.resize(rows());
  result.rs_ = rows();
  result.cs_ = columns();

  // Now, iterate over non-zero rows of this.
  unsigned int row_id = 0;
  for (typename vcl_vector<row>::const_iterator row_iter = elements.begin();
       row_iter != elements.end();
       ++row_iter, ++row_id)
    {
      // Get the row from this matrix (lhs).
      row const & this_row = *row_iter;

      // Get the new row in the result matrix.
      row& result_row = result.elements[row_id];

      // Store this into result row.
      result_row = this_row;

      // If rhs row is empty, we are done.
      if (rhs.empty_row(row_id))
        continue;

      // Get the rhs row.
      row const& rhs_row = rhs.elements[row_id];

      // Iterate over the rhs row.
      for (typename row::const_iterator col_iter = rhs_row.begin();
           col_iter != rhs_row.end();
           ++col_iter)
        {
          // Get the element from the row.
          vnl_sparse_matrix_pair<T> const& entry = *col_iter;
          unsigned const col_id = entry.first;

          // So we are at (row_id,col_id) in rhs matrix.
          result(row_id,col_id) += entry.second;
        }
    }
}

//------------------------------------------------------------
//: Subtract rhs from this.
template <class T>
void vnl_sparse_matrix<T>::subtract(const vnl_sparse_matrix<T>& rhs,
                                    vnl_sparse_matrix<T>& result) const
{
  assert((rhs.rows() == rows()) && (rhs.columns() == columns()));

  // Clear result matrix.
  result.elements.clear();

  // Now give the result matrix enough rows.
  result.elements.resize(rows());
  result.rs_ = rows();
  result.cs_ = columns();

  // Now, iterate over non-zero rows of this.
  unsigned int row_id = 0;
  for (typename vcl_vector<row>::const_iterator row_iter = elements.begin();
       row_iter != elements.end();
       ++row_iter, ++row_id)
    {
      // Get the row from this matrix (lhs).
      row const& this_row = *row_iter;

      // Get the new row in the result matrix.
      row& result_row = result.elements[row_id];

      // Store this into result row.
      result_row = this_row;

      // If rhs row is empty, we are done.
      if (rhs.empty_row(row_id))
        continue;

      // Get the rhs row.
      row const& rhs_row = rhs.elements[row_id];

      // Iterate over the rhs row.
      for (typename row::const_iterator col_iter = rhs_row.begin();
           col_iter != rhs_row.end();
           ++col_iter)
        {
          // Get the element from the row.
          vnl_sparse_matrix_pair<T> const& entry = *col_iter;
          unsigned const col_id = entry.first;

          // So we are at (row_id,col_id) in rhs matrix.
          result(row_id,col_id) -= entry.second;
        }
    }
}

//------------------------------------------------------------
//: Get a reference to an entry in the matrix.
template <class T>
T& vnl_sparse_matrix<T>::operator()(unsigned int r, unsigned int c)
{
  assert((r < rows()) && (c < columns()));
  row& rw = elements[r];
  typename row::iterator ri;
  for (ri = rw.begin(); (ri != rw.end()) && ((*ri).first < c); ++ri);

  if ((ri == rw.end()) || ((*ri).first != c)) {
    // Add new column to the row.
    ri = rw.insert(ri, vnl_sparse_matrix_pair<T>(c,T(0)));
  }

  return (*ri).second;
}

template <class T>
void vnl_sparse_matrix<T>::diag_AtA(vnl_vector<T> & result) const
{
  result.resize( columns() );
  result.fill(T(0));

  typename vcl_vector<row>::const_iterator row_iter = elements.begin();
  for ( ; row_iter != elements.end(); ++row_iter) {
    row const& this_row = *row_iter;
    typename row::const_iterator col_iter = this_row.begin();
    for ( ; col_iter != this_row.end(); ++col_iter) {
      vnl_sparse_matrix_pair<T> const& entry = *col_iter;
      unsigned const col_id = entry.first;
      result[col_id] += entry.second * entry.second;
    }
  }
}

//------------------------------------------------------------
//: Set row in the matrix.

template <class T>
void vnl_sparse_matrix<T>::set_row(unsigned int r,
                                   vcl_vector<int> const& cols,
                                   vcl_vector<T> const& vals)
{
  assert (r < rows());
  assert (cols.size() == vals.size());

  row& rw = elements[r];
  if (rw.size() != cols.size()) rw = row(cols.size());
  for (unsigned int i=0; i < cols.size(); ++i)
    rw[i] = vnl_sparse_matrix_pair<T>(cols[i], vals[i]);
  vcl_sort(rw.begin(), rw.end(), vnl_sparse_matrix_pair<T>::less());
}

template <class T>
vnl_sparse_matrix<T>& vnl_sparse_matrix<T>::vcat(vnl_sparse_matrix<T> const& A)
{
  if (rs_ == 0) {
    rs_ = A.rs_;
    cs_ = A.cs_;
    elements = A.elements;
  }
  else {
    assert(cs_ == A.cs_);
    rs_ += A.rs_;
    elements.insert(elements.end(), A.elements.begin(), A.elements.end());
  }
  return *this;
}


//------------------------------------------------------------
//: This is occasionally useful.  Sums a row of the matrix efficiently.
template <class T>
T vnl_sparse_matrix<T>::sum_row(unsigned int r)
{
  assert(r < rows());
  row & rw = elements[r];
  T sum = T(0);
  for (typename row::iterator ri = rw.begin(); ri != rw.end(); ++ri)
    sum += (*ri).second;

  return sum;
}

template <class T>
void vnl_sparse_matrix<T>::scale_row(unsigned int r, T scale)
{
  assert(r < rows());
  row& rw = elements[r];
  for (typename row::iterator ri = rw.begin(); ri != rw.end(); ++ri)
    (*ri).second *= scale;
}

//------------------------------------------------------------
//: Resizes the matrix so that it has r rows and c columns.
//    Currently not implemented.
//
template <class T>
void vnl_sparse_matrix<T>::resize( int /*r*/, int /*c*/)
{
  vcl_cerr << "Warning: vnl_sparse_matrix::resize not implemented.\n";
  vcl_abort();
}

//------------------------------------------------------------
//: Resets the internal iterator
template <class T>
void vnl_sparse_matrix<T>::reset()
{
  itr_isreset = true;
  itr_row = 0;
}

//------------------------------------------------------------
//: Moves the internal iterator to next non-zero entry in matrix.
// Returns true if there is another value, false otherwise. Use
// in combination with methods reset, getrow, getcolumn, and value.
//
template <class T>
bool vnl_sparse_matrix<T>::next()
{
  if ( itr_row >= rows() )
    return false;

  if ( itr_isreset ) {
    // itr_cur is not pointing to a entry
    itr_row = 0;
    itr_isreset = false;
  } else {
    // itr_cur is pointing to an entry.
    // Try to move to next entry in current row.
    itr_cur++;
    if ( itr_cur != elements[itr_row].end() )
      return true;  // found next entry in current row
    else
      itr_row++;
  }

  // search for next entry starting at row itr_row
  while ( itr_row < rows() ) {
    itr_cur = elements[itr_row].begin();
    if ( itr_cur != elements[itr_row].end() )
      return true;
    else
      itr_row++;
  }

  return itr_row < rows();
}

//------------------------------------------------------------
//: Returns the row of the entry pointed to by internal iterator.
//
template <class T>
int vnl_sparse_matrix<T>::getrow()
{
  return itr_row;
}

//------------------------------------------------------------
//: Returns the column of the entry pointed to by internal iterator.
//
template <class T>
int vnl_sparse_matrix<T>::getcolumn()
{
  return (*itr_cur).first;
}

//------------------------------------------------------------
//: Returns the value pointed to by the internal iterator.
//
template <class T>
T vnl_sparse_matrix<T>::value()
{
  return (*itr_cur).second;
}

#define VNL_SPARSE_MATRIX_INSTANTIATE(T) \
template class vnl_sparse_matrix<T >

#endif // vnl_sparse_matrix_txx_