matpackII.h

Go to the documentation of this file.

/* Copyright (C) 2001 Stefan Buehler <sbuehler@uni-bremen.de>
 
   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.
 
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */
 
#include <vector>
 
/*
  Notes:
 
  There are two different ways to index: 
  S.rw(3,4) = 1;                // Read and write
  cout << S.ro(3,4);            // Read only
 
  This distinction is necessary, because rw() creates elements if they
  don't already exist.
 
 */
 
 
class SparseView {
public:
  // Member functions:
  Index nrows() const;
  Index ncols() const;
  Index nnz()   const;
 
  // Index Operators:
  Numeric& rw(Index r, Index c);
  Numeric  ro(Index r, Index c) const;
 
  // Friends:
  friend std::ostream& operator<<(std::ostream& os, const SparseView& v);
 
protected:
  // Constructors:
  SparseView();
  SparseView(std::vector<Numeric> *data,
                  std::vector<Index>   *rowind,
                  std::vector<Index>   *colptr,
                  const Range&         rr,
                  const Range&         cr                 );
 
  // Data members:
  std::vector<Numeric> *mdata;
  std::vector<Index> *mrowind;
  std::vector<Index> *mcolptr;  
  Index mnr;
};
 
class Sparse : public SparseView {
public:
  // Constructors:
  Sparse();
  Sparse(Index r, Index c);
  Sparse(const Sparse& m);
 
  // Destructor:
  ~Sparse();
};
 
 
// Functions for SparseView:
// -------------------------------
 
inline Index SparseView::nrows() const
{
  return mrr.mextent;
}
 
inline Index SparseView::ncols() const
{
  return mcr.mextent;
}
 
inline Index SparseView::nnz() const
{
  return *(mcolptr->end()-1);
}
 
inline Numeric& SparseView::rw(Index r, Index c)
{
  // Check if indices are valid:
  assert( 0<=r );
  assert( 0<=c );
  assert( r<mrr.mextent );
  assert( c<mcr.mextent );
 
  // Convert to true indices. We use r and c, which are local
  // variables because we used call be value.
  r = mrr.mstart + r*mrr.mstride;
  c = mcr.mstart + c*mcr.mstride;
  
  // Get index of first data element of this column:
  Index i = (*mcolptr)[c];
 
  // Get index of first data element of next column. (This always works,
  // because mcolptr has one extra element pointing behind the last
  // column.)
  const Index end = (*mcolptr)[c+1];
 
  // See if we find an element with the right row index in this
  // range. We assume that the elements are sorted by ascending row
  // index. We use the index i as the loop counter, which we have
  // initialized above. 
  for ( ; i<end; ++i )
    {
      Index rowi = (*mrowind)[i]; 
      if ( r <  rowi )
        break;
      if ( r == rowi )
        return (*mdata)[i];
    }
 
  // If we are here, then the requested data element does not yet
  // exist. 
 
  // We have to adjust the array of column pointers. The values
  // in all columns above the current one have to be increased by 1. 
  for ( std::vector<Index>::iterator j = mcolptr->begin() + c + 1;
        j < mcolptr->end();
        ++j )
    ++(*j);
 
  // We have to insert the new element in *mrowind and *mdata just one
  // position before the index i. We can use std::insert to achieve
  // this. Because they return an iterator to the newly inserted
  // element, we can return directly from the second call.
  mrowind->insert( mrowind->begin()+i, r );
  return *( mdata->insert( mdata->begin()+i, 0 ) );
}
 
inline Numeric SparseView::ro(Index r, Index c) const
{
  // Check if indices are valid:
  assert( 0<=r );
  assert( 0<=c );
  assert( r<mrr.mextent );
  assert( c<mcr.mextent );
 
  // Convert to true indices:
  r = mrr.mstart + r*mrr.mstride;
  c = mcr.mstart + c*mcr.mstride;
  
  // Get index of first data element of this column:
  Index begin = (*mcolptr)[c];
 
  // Get index of first data element of next column. (This always works,
  // because mcolptr has one extra element pointing behind the last
  // column.)
  const Index end = (*mcolptr)[c+1];
 
  // See if we find an element with the right row index in this
  // range. We assume that the elements are sorted by ascending row
  // index. 
  for ( Index i=begin; i<end; ++i )
    {
      Index rowi = (*mrowind)[i]; 
      if ( r <  rowi )
        return 0;
      if ( r == rowi )
        return (*mdata)[i];
    }
  return 0;
}
 
inline SparseView::SparseView() :
  mdata(NULL),
  mrowind(NULL),
  mcolptr(NULL),
  mrr(0,0),
  mcr(0,0)
{
  // Nothing to do here
}
 
inline SparseView::SparseView(std::vector<Numeric> *data,   
                                        std::vector<Index>   *rowind, 
                                        std::vector<Index>   *colptr, 
                                        const Range&         rr,      
                                        const Range&         cr       ) :
  mdata(data),
  mrowind(rowind),
  mcolptr(colptr),
  mrr(rr),
  mcr(cr)
{
  // Nothing to do here.
}
 
 
 
// Functions for Sparse:
// ---------------------------
 
inline Sparse::Sparse() 
{
  // Nothing to do here
}
 
inline Sparse::Sparse(Index r, Index c) :
  SparseView( new std::vector<Numeric>,
                   new std::vector<Index>,
                   new std::vector<Index>(c+1,0),
                   Range(0,r),
                   Range(0,c) )
{
  // Nothing to do here.
}
 
inline Sparse::Sparse(const Sparse& m) :
  SparseView( new std::vector<Numeric>(*m.mdata),
                   new std::vector<Index>(*m.mrowind),
                   new std::vector<Index>(*m.mcolptr),
                   m.mrr,
                   m.mcr )
{
  // Nothing to do here. 
}
 
 
inline Sparse::~Sparse()
{
  delete mdata;
  delete mrowind;
  delete mcolptr;
}
 
// Output operator for SparseView:
 
inline std::ostream& operator<<(std::ostream& os, const SparseView& v)
{
  for (size_t c=0; c<v.mcolptr->size()-1; ++c)
    {
      // Get index of first data element of this column:
      Index begin = (*v.mcolptr)[c];
 
      // Get index of first data element of next column. (This always works,
      // because mcolptr has one extra element pointing behind the last
      // column.)
      const Index end = (*v.mcolptr)[c+1];
 
      // Loop through the elements in this column:
      for ( Index i=begin; i<end; ++i )
        {
          Index r = (*v.mrowind)[i]; 
 
          // Now we know the true row r and column c. Convert to aparent r
          // and c using the Ranges mrr and mcr,
 
          Index ra = (r-v.mrr.mstart)/v.mrr.mstride;
          Index ca = (c-v.mcr.mstart)/v.mcr.mstride;
 
          // Are the aparent row ra and column ca inside the active range?
          if ( 0 <= ra &&
               ra < v.mrr.mextent &&
               0 <= ca &&
               ca < v.mcr.mextent     )
            {
              // Yes, they are! Let's output this element.
              os << setw(3) << ra << " "
                 << setw(3) << ca << " "
                 << setw(3) << (*v.mdata)[i] << "\n";
            }
        }
    }
 
  return os;
}