compressed2d.h

MTL C++ Numeric Library

  

  //: Return an iterator pointing to the first 1D container
  inline iterator begin() {
    return iterator(values, indices, starts, 0);
  }
  //: Return an iterator pointing past the end of the 2D container
  inline iterator end() {
    return iterator(values, indices, starts, dim.first());
  }

  //: Return a const iterator pointing to the first 1D container
  inline const_iterator begin() const {
    return const_iterator(values, indices, starts, 0);
  }
  //: Return a const iterator pointing past the end of the 2D container
  const_iterator end() const {
    return const_iterator(values, indices, starts, dim.first());
  }

  /* reverse iterators */

  //: Return a reverse iterator pointing to the last 1D container
  inline reverse_iterator rbegin() {
    return reverse_iterator(end());
  }
  //: Return a reverse iterator pointing past the start of the 2D container
  inline reverse_iterator rend() {
    return reverse_iterator(begin());
  }

  //: Return a const reverse iterator pointing to the last 1D container
  inline const_reverse_iterator rbegin() const {
    return const_reverse_iterator(end());
  }
  //: Return a const reverse iterator pointing past the start of the 2D container
  inline const_reverse_iterator rend() const {
    return const_reverse_iterator(begin());
  }

  /* Element Access Methods */
  
  //: Return a reference to the (i,j) element, where (i,j) is in the 2D coordinate system
  inline typename reference::reference 
  operator()(size_type i, size_type j) {
    return MajorVectorRef(values, indices, starts, i)[j];
  }
  //: Return a const reference to the (i,j) element, where (i,j) is in the 2D coordinate system
  inline typename const_reference:: const_reference
  operator()(size_type i, size_type j) const {
    return ConstMajorVectorRef(values, indices, starts, i)[j];
  }
  

  /* Size Methods */
  
  //: The dimension of the 2D container
  inline size_type major() const { return dim.first(); }
  //: The dimension of the 1D containers
  inline size_type minor() const { return dim.second(); } 
  //: The number of non-zeros 
  inline size_type nnz() const { return values->size(); }
  
  /* Vector Access Methods */

  //: Return the ith 1D container
  inline value_type operator[](size_type i) const {
    return MajorVector(values, indices, starts, i);
  }

  /* return the raw data */

  //: Return a pointer to the values array
  inline TT* get_val() { return values->data(); }
  //:
  inline const TT* get_val() const { return values->data(); }

  //: Return a pointer to the indices array
  inline size_type* get_ind() { return indices->data(); }
  //:
  inline const size_type* get_ind() const { return indices->data(); }

  //: Return a pointer to the pointers array
  inline size_type* get_ptr() { return starts->data(); }
  //:
  inline const size_type* get_ptr() const { return starts->data(); }

  //: A fast specialization for copying from a sparse matrix
  template <class SparseMat>
  void fast_copy(const SparseMat& x, sparse_tag) {
    size_type nnz = x.nnz();
    typename SparseMat::const_iterator i;
    typename SparseMat::OneD::const_iterator j, jend;

    values->resize(nnz);
    indices->resize(nnz);
    starts->resize(x.major() + 1);

    size_type curr = 0;

    for (i = x.begin(); i != x.end(); ++i) {
      (*starts)[i.index()] = curr - IND_OFFSET; /* F to C */
      jend = (*i).end();
      for (j = (*i).begin(); j != jend; ++j) {
        (*values)[curr] = *j;
        (*indices)[curr++] = j.index() - IND_OFFSET; /* F to C */
      }
    }
    (*starts)[i.index()] = curr - IND_OFFSET; /* F to C */
  }
  //: A fast specialization for copying from a dense matrix
  template <class DenseMat>
  void fast_copy(const DenseMat x, dense_tag) {
    typename DenseMat::const_iterator i;
    typename DenseMat::OneD::const_iterator j;

    size_type nnz = 0;

    /* count non-zeros in each row */
    for (i = x.begin(); i != x.end(); ++i)
      for (j = (*i).begin(); j != (*i).end(); ++j)
        if (*j != TT(0))
          ++nnz;

    values->resize(nnz);
    indices->resize(nnz);
    starts->resize(x.major() + 1);

    size_type curr = 0;

    for (i = x.begin(); i != x.end(); ++i) {
      (*starts)[i.index()] = curr - IND_OFFSET; /* F to C */
      for (j = (*i).begin(); j != (*i).end(); ++j) {
        if (*j != TT(0)) {
          (*values)[curr] = *j;
          (*indices)[curr++] = j.index() - IND_OFFSET; /* F to C */
        }
      }
    }
    (*starts)[i.index()] = curr - IND_OFFSET; /* F to C */
  }
  template <class TwoD__>
  void fast_copy(const TwoD__& x) {
    typedef typename TwoD__::sparsity Sparsity;
    fast_copy(x, Sparsity());
  }

  void print() const {
    cout << "values ";
    print_vector(*values);
    cout << "indices ";
    print_vector(*indices);
    cout << "starts ";
    print_vector(*starts);
  }

protected:
  dim_type dim;
  ValPtr values;
  IndPtr indices;
  IndPtr starts;  /* starting point for each minor vector */
};

//: Sparse matrix storage format with internal storage
//
// This version of the compressed matrix format keeps track 
// of its own memory using reference counting.
//
//!component: type
//!category: container
//!definition: compressed2D.h
//!models: TwoDContainerRef
//!tparam: T - The element type
//!tparam: SizeType - The type for the stored indices
//!tparam: IND_OFFSET - To handle indexing from 0 or 1
template <class T, class SizeType, int IND_OFFSET>
class compressed2D : public generic_comp2D< dense1D<T>, dense1D<T>*,
                     dense1D<SizeType>,
                     dense1D<SizeType>* , IND_OFFSET>
{
  //: The base class
  typedef generic_comp2D< dense1D<T>, dense1D<T>*,
                     dense1D<SizeType>,
                     dense1D<SizeType>* , IND_OFFSET> Base;
  //: The type of this class
  typedef compressed2D<T, SizeType,IND_OFFSET> self;
public:
  typedef typename Base::dim_type dim_type;

  //: The integral type for dimensions, indices, etc.
  typedef SizeType size_type;
  //: A variant of a pair type to describe the band width
  typedef dimension<int> band_type;

  //: Default Constructor
  inline compressed2D() : Base(dim_type(0,0), &vals, &inds, &ptrs) { }

  //: Normal Constructor
  inline compressed2D(dim_type d) 
    : Base(d, &vals, &inds, &ptrs),
      ptrs(d.first() + 1, -IND_OFFSET) /* F to C */
  {
    vals.reserve(dim.first() * 5);
    inds.reserve(dim.first() * 5);
  }

  inline compressed2D(dim_type d, size_type nnz)
    : Base(d, &vals, &inds, &ptrs),
      ptrs(d.first() + 1, -IND_OFFSET) /* F to C */
  {
    vals.reserve(nnz);
    inds.reserve(nnz);
  }

  //: Banded Constructor
  inline compressed2D(dim_type d, band_type) 
    : Base(d, &vals, &inds, &ptrs), 
      ptrs(d.first() + 1, -IND_OFFSET) /* F to C */
  {
    vals.reserve(dim.first() * 5);
    inds.reserve(dim.first() * 5);
  }

  //: Copy Constructor
  inline compressed2D(const self& x)
    : Base(x.dim, &vals, &inds, &ptrs), 
      vals(x.vals), ptrs(x.ptrs), inds(x.inds)
  { }

  // for banded view
  inline compressed2D(const self& x, band_type, banded_tag)
    : Base(x.dim, &vals, &inds, &ptrs), 
      vals(x.vals), ptrs(x.ptrs), inds(x.inds) 
  { }

  inline self& operator=(const self& x) {
    vals = x.vals; ptrs = x.ptrs; inds = x.inds;
    /* fill out inhereted part */
    dim = x.dim;
    values = &vals; indices = &inds; starts = &ptrs;
    return *this;
  }

  //: Matrix Stream Constructor
  template <class MatrixStream, class Orien>
  inline compressed2D(MatrixStream& s, Orien) /* F to C */
    : Base(Orien::map(dim_type(s.nrows(),s.ncols())), &vals, &inds, &ptrs),
      ptrs(1 + Orien::map(dim_type(s.nrows(),s.ncols())).first(), -IND_OFFSET)
  { 
    vals.reserve(s.nnz());
    inds.reserve(s.nnz());
  }

  //: Banded Matrix Stream Constructor
  template <class MatrixStream, class Orien>
  inline compressed2D(MatrixStream& s, Orien, band_type) /* F to C */
    : Base(Orien::map(dim_type(s.nrows(),s.ncols())), &vals, &inds, &ptrs),
      ptrs(1 + Orien::map(dim_type(s.nrows(),s.ncols())).first(), -IND_OFFSET)
  { 
    vals.reserve(s.nnz());
    inds.reserve(s.nnz());
  }
#if 0
  // deprecated
  template <class SubMatrix>
  struct partitioned {
    typedef compressed2D<SubMatrix, size_type,IND_OFFSET> type;
    typedef gen_compressed2D<SubMatrix, size_type,IND_OFFSET> generator;
  };
#endif
  inline size_type capacity() const { return inds.capacity(); }

protected:
  dense1D<T> vals;
  dense1D<size_type> ptrs;
  dense1D<size_type> inds;
};


//: Sparse matrix storage format with external storage
//
// This version of the compressed matrix format uses memory
// provided by the user for the matrix storage.
//
//!component: type
//!category: container
//!definition: compressed2D.h
//!models: TwoDContainerRef
//!tparam: T - The element type
//!tparam: SizeType - The type for the stored indices
//!tparam: IND_OFFSET - To handle indexing from 0 or 1
template <class T, class SizeType, int IND_OFFSET>
class ext_comp2D : public generic_comp2D< external_vec<T,0,SizeType>,
                                          external_vec<T,0,SizeType>*,
                                          external_vec<SizeType,0,SizeType>, 
                                          external_vec<SizeType,0,SizeType>*,
                             IND_OFFSET >
{
  //: The type of this class
  typedef ext_comp2D<T, SizeType,IND_OFFSET> self;
  //: The base class
  typedef generic_comp2D< external_vec<T,0,SizeType>, 
                          external_vec<T,0,SizeType>*,
                          external_vec<SizeType,0,SizeType>,
                          external_vec<SizeType,0,SizeType>*,
                          IND_OFFSET > Base;
public:
  typedef typename Base::dim_type dim_type;

  external_vec<T,0,SizeType> vals;
  typedef SizeType size_type;
  external_vec<size_type,0,SizeType> inds;
  external_vec<size_type,0,SizeType> ptrs;

  typedef typename external_vec<size_type,0,SizeType>::size_type rep_size_t;

  inline ext_comp2D() { }
  //: Preallocated Memory Constructor
  inline ext_comp2D(dim_type d, size_type nz,
                    T* val, size_type* ptr, size_type* ind)
    : Base(d, &vals, &inds, &ptrs),
      vals(val, nz), 
      inds(ind, nz),
      ptrs(ptr, rep_size_t(d.first() + 1))
  { }

  //: Copy Constructor
  inline ext_comp2D(const self& x)
    : Base(dim_type(x.major(), x.minor()), &vals, &inds, &ptrs),
      vals(x.vals), 
      inds(x.inds), 
      ptrs(x.ptrs)
  { }

  inline self& operator=(const self& x) {
    vals = x.vals; inds = x.inds; ptrs = x.ptrs;
    /* fill out inhereted part */
    dim = x.dim;
    values = &vals; indices = &inds; starts = &ptrs;
    return *this;
  }

  //: Banded View (including Symm and Tri) Constructor
  typedef dimension<int> band_type;
  inline ext_comp2D(const self& x, band_type, banded_tag)
    : Base(x.dim, &vals, &inds, &ptrs),
      vals(x.vals), 
      inds(x.inds), 
      ptrs(x.ptrs)
  { }
#if 0
  // deprecated
  template <class SubMatrix>
  struct partitioned {
    typedef ext_comp2D<SubMatrix, size_type,IND_OFFSET> type;
    typedef gen_ext_comp2D<SubMatrix, size_type,IND_OFFSET> generator;
  };
#endif
};

//: blah
//!noindex:
template <class T, class SizeType, int INDEX, int M, int N>
struct gen_compressed2D {
  typedef gen_compressed2D<T,SizeType,INDEX,M,N> submatrix_type; /* bogus */
  typedef gen_dense2D<T,gen_rect_offset<M,N>,N,M> transpose_type; /* bogus */
  typedef gen_dense2D<T,gen_rect_offset<M,N>,M,N> banded_view_type; /* bogus */

  typedef compressed2D<T, SizeType,INDEX> type;
};



//: blah
//!noindex:
template <class T, class SizeType, int INDEX, int M, int N>
struct gen_ext_comp2D {
  typedef gen_ext_comp2D<T, SizeType,INDEX,M,N> submatrix_type; /* bogus */
  typedef gen_dense2D<T,gen_rect_offset<N,M>,N,M> transpose_type; /* bogus */
  typedef gen_dense2D<T,gen_rect_offset<M,N>,M,N> banded_view_type; /* bogus */
  typedef ext_comp2D<T, SizeType,INDEX> type;
};




} /* namespace mtl */

#endif




#if 0
  template <class Sparse2D>
  inline compressed2D(const nz_structure<Sparse2D>& x)
    : major_(x.major()), minor_(x.minor()),
      starts(1 + x.major())
  {
    int nnz = x.nnz();
    values->resize(nnz);
    indices->resize(nnz);

    int curr = 0;
    int i;
    for (i = 0; i < major_; ++i) {
      starts[i] = curr;
      
      typename Sparse2D::MajorVector::const_iterator j = x[i].begin();
      typename Sparse2D::MajorVector::const_iterator jend = x[i].end();
      while (not_at(j, jend)) {
        indices[curr++] = j.index();
        ++j;
      }
    }
    starts[i] = curr;
  }
#endif

#if 0
    inline TT get(int i) const {
      int first = (*starts)[major];
      int last = (*starts)[major + 1];
      const_index_iterator indices_begin = 
        ((const indices_t*)indices)->begin();
      const_index_iterator evi = find(indices_begin + first,
                                      indices_begin + last, i);
      return evi == (indices_begin + last) ?
        0.0 : (*values)[evi - indices_begin];
    }
    
    inline void set(int i, TT v) {
      int first = (*starts)[major];
      int last = (*starts)[major + 1];
      
      if (first == last) {
        /* this row had no entries */
        (*starts)[i] = indices->insert(indices->begin() + first, i)
          - indices->begin();
        values->insert(values->begin() + first, v);
        increment_starts(major + 1);
      }
      else {
        index_iterator evi = lower_bound(indices->begin() + first,
                                         indices->beign() + last, i);
        if (evi == (indices->begin() + last)) {
          evi = indices->insert(evi, i);
          int n = evi - indices->begin();
          values->insert(values->begin() + n, v);       
          increment_starts(major + 1);
        }
        else if (*evi == j) {
          int n = evi - indices->begin();
          (*values)[n] = v;
        }
        else {
          evi = indices->insert(evi, i);
          int n = evi - indices->begin();
          values->insert(values->begin() + n, v);       
          increment_starts(major + 1);
        }
      }
    }
#endif