compressed2d.h

MTL C++ Numeric Library

// -*- c++ -*-
//
// Copyright 1997, 1998, 1999 University of Notre Dame.
// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee
//
// This file is part of the Matrix Template Library
//
// You should have received a copy of the License Agreement for the
// Matrix Template Library along with the software;  see the
// file LICENSE.  If not, contact Office of Research, University of Notre
// Dame, Notre Dame, IN  46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//
//===========================================================================

#ifndef MTL_COMPRESSED_CONTIG2D_H
#define MTL_COMPRESSED_CONTIG2D_H

#include "mtl/mtl_iterator.h"
#include <algorithm>

#include "mtl/entry.h"
#include "mtl/compressed_iter.h"
#include "mtl/reverse_iter.h"
#include "mtl/linalg_vec.h"
#include "mtl/matrix_traits.h"
#include "mtl/dense2D.h"
#include "mtl/dimension.h"
#include "mtl/dense1D.h"

#include "mtl/utils.h" /* debug */

extern "C" void foobar(int);

namespace mtl {

using std::less;
using std::lower_bound;

template <class T, class SizeType, int IndexOff, int M, int N>
struct gen_compressed2D;
template <class T, class SizeType, int IndexOff, int M, int N>
struct gen_ext_comp2D;

//: Compressed 2-D Container
//
//   The compressed2D container uses the compressed row/column storage
//   format. There is an array of matrix elements, an array of indices,
//   and an array of pointers to the row or column vector starts.
//   A compressed2D matrix can be created from scratch
//   with the constructor:
//   <codeblock>
//   compressed2D(size_type m, size_type n)
//   </codeblock>
//   One can also use preexsisting arrays with the constructor:
//   <codeblock>
//   compressed2D(size_type m, size_type n, size_type nnz,
//                T* val, size_type* ptr, size_type* ind)
//   </codeblock>
//
//   The stored indices (in the ptr and int arrays) are indexed
//   from 1 ala LAPACK and BLAS conventions (Fortran style).
//
//!category: containers
//!component: type
//!definition: compressed2D.h
//!models: TwoDContainerRef
//!example: sparse_matrix.h
//!tparam: ValsType - The container type to use for the values array
//!tparam: ValPtr - A pointer to the ValsType
//!tparam: IndType - The container type to use for the indices and pointer arrays
//!tparam: IndPtr - A pointer to IndType
//!tparam: IND_OFFSET - To handle indexing from 0 or 1
template <class ValsType, class ValPtr, class IndType, class IndPtr, int IND_OFFSET>
class generic_comp2D {
public:
  typedef generic_comp2D self;
  typedef ValsType values_t;
  typedef typename ValsType::value_type TT;
  typedef typename values_t::iterator        value_iterator;
  typedef typename values_t::const_iterator  const_value_iterator;
  typedef IndType indices_t;
  typedef typename indices_t::iterator       index_iterator;
  typedef typename indices_t::const_iterator const_index_iterator;
  typedef IndType starts_t;
  typedef typename starts_t::iterator        starts_iterator;
  typedef typename starts_t::const_iterator  const_starts_iterator;

  typedef typename IndType::value_type size_type;

  typedef dimension<size_type> dim_type;

  /* Type Definitions */
  typedef internal_tag storage_loc;
  enum { M = 0, N = 0 };

  //: This vector reference is created on-the-fly as needed.
  class vec_ref {
  public:
    enum { N = 0 };
    typedef TT value_type;
    typedef TT* pointer;
    typedef typename IndType::value_type size_type; 
    typedef typename IndType::difference_type difference_type;
    typedef elt_ref<vec_ref> reference;
    typedef const_elt_ref<vec_ref> const_reference;
    typedef sparse_tag sparsity;
    typedef oned_tag dimension;
    typedef compressed_iter<0,values_t, indices_t, IND_OFFSET> iterator;
    typedef compressed_iter<1,values_t, indices_t, IND_OFFSET> const_iterator;

    typedef reverse_iter< iterator > reverse_iterator;
    typedef reverse_iter< const_iterator > const_reverse_iterator;

    typedef light1D<size_type,0,IND_OFFSET> IndexArrayRef;
    typedef light1D<size_type,0,IND_OFFSET> IndexArray;
    typedef vec_ref subrange_type;

    inline vec_ref(ValPtr val, IndPtr ind, IndPtr s, size_type major_)
      : values(val), indices(ind), starts(s), major(major_) {
    }
    inline vec_ref(const vec_ref& x)
      : values(x.values), indices(x.indices), starts(x.starts),
        major(x.major) { }

    inline iterator begin() {
      return iterator(values->begin(), indices->begin(),
                      (*starts)[major] + IND_OFFSET); /* F to C */
    }
    inline iterator end() { 
      return iterator(values->begin(), indices->begin(),
                      (*starts)[major+1] + IND_OFFSET); /* F to C */
    }
    inline const_iterator begin() const {
      return const_iterator(((const values_t*)values)->begin(), 
                            ((const indices_t*)indices)->begin(),
                            (*starts)[major] + IND_OFFSET); /* F to C */
    }
    inline const_iterator end() const {
      return const_iterator(((const values_t*)values)->begin(), 
                            ((const indices_t*)indices)->begin(),
                            (*starts)[major+1] + IND_OFFSET); /* F to C */
    }

    inline reverse_iterator rbegin() {
      return reverse_iterator(end()); 
    }
    inline reverse_iterator rend() {
      return reverse_iterator(begin()); 
    }
    inline const_reverse_iterator rbegin() const {
      return const_reverse_iterator(end());
    }
    inline const_reverse_iterator rend() const {
      return const_reverse_iterator(begin()); 
    }

    inline reference operator[](size_type i) {
      return reference(*this, i);
    }

    inline const_reference operator[](size_type i) const {
      return const_reference(*this, i);
    }

    //add this for BEAM
    value_type& get_ref(size_type i) {
      iterator iter = find(i);
      /*
      if ( iter != end() ) { 
	if ( iter.index() != i )
	  iter = insert(iter, i, value_type(0));
      } else 
	iter = insert(iter, i, value_type(0));
      */
	  assert(iter.index() == i);
      return *iter;
    }

    /* JGS this is wrong */
    inline size_type size() const { return nnz(); }

    inline size_type nnz() const { 
      return (*starts)[major+1] - (*starts)[major]; 
    }

    inline IndexArrayRef nz_struct() const {
      size_type first = (*starts)[major] + IND_OFFSET; /* F to C */
      size_type last = (*starts)[major + 1] + IND_OFFSET; /* F to C */
      return IndexArrayRef(&(*indices)[0] + first, last - first);
    }

    inline iterator find(size_type i) {
      difference_type first = (*starts)[major] + IND_OFFSET; /* F to C */
      difference_type last = (*starts)[major + 1] + IND_OFFSET; /* F to C */
      index_iterator indices_begin = indices->begin();
      index_iterator iter = lower_bound(indices_begin + first,
                                        indices_begin + last, 
                                        size_type(i - IND_OFFSET)); /* F to C */
      size_type n = iter - indices_begin;
      return iterator(values->begin(), indices_begin, n);
    }
    
    inline const_iterator find(size_type i) const {
      difference_type first = (*starts)[major] + IND_OFFSET; /* F to C */
      difference_type last = (*starts)[major + 1] + IND_OFFSET; /* F to C */
      const_index_iterator indices_begin =
        ((const indices_t*)indices)->begin();
      const_index_iterator iter = lower_bound(indices_begin + first,
                                              indices_begin + last, 
                                              i - IND_OFFSET); /* F to C */
      size_type n = iter - indices_begin;
      return const_iterator(((const values_t*)values)->begin(), 
                            indices_begin, n);
    }
    
    inline iterator insert(iterator iter, size_type i, TT v) {
      index_iterator ind = indices->insert(iter.index_iter(),
                                            size_type(i - IND_OFFSET)); /* F to C */
      values->insert(iter.value_iter(), v);
      size_type n = ind - indices->begin();
      increment_starts(major + 1);
      return iterator(values->begin(), indices->begin(), n);
    }
    
    inline void increment_starts(size_type i) {
      while (i < size_type(starts->size())) 
        ++(*starts)[i++];
    }

    inline void resize(size_type) { }

  private:

    ValPtr values;
    IndPtr indices;
    IndPtr starts;
    size_type major;
  };

public:

  //: The 1D container type
  typedef vec_ref value_type;
  //: A reference to the value type
  typedef vec_ref reference;
  //: A const reference to the value type
  typedef const vec_ref const_reference;

  /* old typedefs */
  typedef vec_ref MajorVector;
  typedef vec_ref MajorVectorRef;
  typedef const vec_ref ConstMajorVectorRef;
public:
  //: Specify that this matrix is sparse
  typedef sparse_tag sparsity;

  //: Specify that this matrix is not strideable (can not use rows(A), columns(A))  
  typedef not_strideable strideability;

  /* bogus types, need to change */
  typedef generic_comp2D<ValsType,ValPtr,IndType,IndPtr,IND_OFFSET> transpose_type;
  typedef generic_comp2D<ValsType,ValPtr,IndType,IndPtr,IND_OFFSET> submatrix_type;
  typedef generic_comp2D<ValsType,ValPtr,IndType,IndPtr,IND_OFFSET> banded_view_type;

#if !defined(__MWERKS__) // internal compiler errors

  //: The type for the iterators
  template <int isConst>
  class _iterator {
    typedef _iterator self;
  public:
    typedef typename IF<isConst, const_value_iterator, value_iterator>::RET
                 value_iter;
    
    typedef typename std::iterator_traits<value_iter>::iterator_category
                 iterator_category;
#if defined(_MSVCPP_)
    typedef typename std::iterator_traits<value_iter>::distance_type
                 difference_type;
    typedef difference_type distance_type;
#else
    typedef typename std::iterator_traits<value_iter>::difference_type
                 difference_type;
#endif

    typedef MajorVectorRef value_type;
    
    typedef typename IF<isConst, const MajorVectorRef, MajorVectorRef>::RET reference;
    typedef typename IF<isConst, const MajorVectorRef*, MajorVectorRef*>::RET pointer;

    typedef typename IF<isConst, const ValPtr, ValPtr>::RET myValPtr;
    typedef typename IF<isConst, const IndPtr, IndPtr>::RET myIndPtr;

    inline _iterator() : pos(0) { }
    inline _iterator(myValPtr val, myIndPtr ind, myIndPtr s, size_type p)
      : values(val), indices(ind), starts(s), pos(p) { }
    inline _iterator(const self& x)
      : values(x.values), indices(x.indices), starts(x.starts), pos(x.pos) { }

    inline size_type index() const { return pos; }

    inline reference operator*() const {
      return reference(values, indices, starts, pos);
    }
    inline self& operator++() { ++pos; return *this; }
    inline self& operator--() { --pos; return *this; }
    inline self& operator+=(size_type n) { pos += n; return *this; }
    inline self& operator-=(size_type n) { pos -= n; return *this; }
    inline difference_type operator-(const self& x) const { 
      return pos - x.pos; }
    inline bool operator<(const self& x) const { return pos < x.pos; }
    inline bool operator!=(const self& x) const { return pos != x.pos; }
    inline bool operator==(const self& x) const { return pos == x.pos; }
  private:
    ValPtr values;
    IndPtr indices;
    IndPtr starts;
    size_type pos;
  };

  typedef _iterator<0> iterator;
  typedef _iterator<1> const_iterator;
  
#else

  //: The type for the iterators
  class iterator {
    typedef iterator self;
  public:
    typedef value_iterator value_iter;
    
    typedef typename std::iterator_traits<value_iter>::iterator_category
                 iterator_category;
#if defined(_MSVCPP_)
    typedef typename std::iterator_traits<value_iter>::distance_type
                 difference_type;
        typedef difference_type distance_type;
#else
    typedef typename std::iterator_traits<value_iter>::difference_type
                 difference_type;
#endif

    typedef MajorVectorRef value_type;
    
    typedef MajorVectorRef reference;
    typedef MajorVectorRef pointer;

    typedef ValPtr myValPtr;
    typedef IndPtr myIndPtr;

    inline iterator() : pos(0) { }
    inline iterator(myValPtr val, myIndPtr ind, myIndPtr s, size_type p)
      : values(val), indices(ind), starts(s), pos(p) { }
    inline iterator(const self& x)
      : values(x.values), indices(x.indices), starts(x.starts), pos(x.pos) { }

    inline size_type index() const { return pos; }

    inline reference operator*() const {
      return reference(values, indices, starts, pos);
    }
    inline self& operator++() { ++pos; return *this; }
    inline self& operator--() { --pos; return *this; }
    inline self& operator+=(size_type n) { pos += n; return *this; }
    inline self& operator-=(size_type n) { pos -= n; return *this; }
    inline difference_type operator-(const self& x) const { 
      return pos - x.pos; }
    inline bool operator<(const self& x) const { return pos < x.pos; }
    inline bool operator!=(const self& x) const { return pos != x.pos; }
    inline bool operator==(const self& x) const { return pos == x.pos; }
  private:
    ValPtr values;
    IndPtr indices;
    IndPtr starts;
    size_type pos;
  };

  //: The type for the iterators
  class const_iterator {
    typedef const_iterator self;
  public:
    typedef const_value_iterator value_iter;
    typedef typename std::iterator_traits<value_iter>::iterator_category
                 iterator_category;
#if defined(_MSVCPP_)
    typedef typename std::iterator_traits<value_iter>::distance_type
                 difference_type;
        typedef difference_type distance_type;
#else
    typedef typename std::iterator_traits<value_iter>::difference_type
                 difference_type;
#endif

    typedef MajorVectorRef value_type;
    
    typedef const MajorVectorRef reference;
    typedef const MajorVectorRef* pointer;

    typedef const ValPtr myValPtr;
    typedef const IndPtr myIndPtr;

    inline const_iterator() : pos(-1) { }
    inline const_iterator(myValPtr val, myIndPtr ind, myIndPtr s, size_type p)
      : values(val), indices(ind), starts(s), pos(p) { }
    inline const_iterator(const self& x)
      : values(x.values), indices(x.indices), starts(x.starts), pos(x.pos) { }

    inline size_type index() const { return pos; }

    inline reference operator*() const {
      return reference(values, indices, starts, pos);
    }
    inline self& operator++() { ++pos; return *this; }
    inline self& operator--() { --pos; return *this; }
    inline self& operator+=(size_type n) { pos += n; return *this; }
    inline self& operator-=(size_type n) { pos -= n; return *this; }
    inline difference_type operator-(const self& x) const { 
      return pos - x.pos; }
    inline bool operator<(const self& x) const { return pos < x.pos; }
    inline bool operator!=(const self& x) const { return pos != x.pos; }
    inline bool operator==(const self& x) const { return pos == x.pos; }
  private:
    ValPtr values;
    IndPtr indices;
    IndPtr starts;
    size_type pos;
  };
#endif

  //: The type for the reverse iterators
  typedef reverse_iter< iterator > reverse_iterator;
  //: The type for the const reverse iterators
  typedef reverse_iter< const_iterator > const_reverse_iterator;

  /* Constructors */  

  //: Default Constructor
  inline generic_comp2D() : dim(0,0), values(0), indices(0), starts(0) { }

  //: External Storage Constructor
  inline generic_comp2D(dim_type d, ValPtr v, IndPtr ind, IndPtr s) 
    : dim(d), values(v), indices(ind), starts(s) { }

  //: Copy Constructor
  inline generic_comp2D(const self& x)
    : dim(x.dim), values(x.values), starts(x.starts), indices(x.indices) { }

  inline self& operator=(const self& x) {
    dim = x.dim; values = x.values; starts = x.starts; indices = x.indices;
    return *this;
  }

  /* Iterator Access Methods */