dense2d.h

MTL C++ Numeric Library

  typedef dimension<size_type, MM, NN> dim_type;
  typedef dimension<int> band_type;
  typedef banded_offset<size_type, MM, NN> transpose_type; /* bogus */
  typedef not_strideable strideability;
  inline packed_offset() 
    : dim(0,0), bw(band_type(0,0)) { }

  inline packed_offset(size_type m, size_type n, size_type /* lead */,
                       band_type bandwidth)
    : dim(m,n), bw(bandwidth) { }

  inline packed_offset& operator=(const packed_offset& x) {
    dim = x.dim; bw = x.bw; return *this;
  }
  inline int elt(size_type i, size_type j) const { 
    return this->oned_offset(i) + j; 
  }

  inline int calc_low(int i, int low) const {
    int l = MTL_MIN(low, int(i));
    int lower_area = low * i;
    lower_area -= ( - l*l + 2*low*l + l) / 2;
    return lower_area;
  }
  inline int calc_up(int i, int up) const {
    int upper_area = up * i;
    int n = i + up - dim.second();
    if (n > 0) {
      int n1 = MTL_MAX(n - up, 0);
      int n2 = n - n1;
      upper_area -= n1 * up;
      upper_area -= ((n2 + 1) * n2) / 2;
    }
    return upper_area;
  }

  inline int oned_offset(size_type i) const { /* the ith major container */
    int low = bw.first();
    int up = bw.second();
    int upper_area, lower_area;

    if (up < -1)
      upper_area = - calc_low(i, - (up + 1));
    else if (up > 0)
      upper_area = calc_up(i, up);
    else
      upper_area = 0;

    if (low < -1)
      lower_area = - calc_up(i, - (low + 1));
    else if (low > 0)
      lower_area = calc_low(i, low);
    else
      lower_area = 0;

    size_type diagonal_len;
    if (up < 0 || low < 0)
      diagonal_len = 0;
    else
      diagonal_len = MTL_MIN(MTL_MIN(i, dim.first()), dim.second());

    size_type ret =  upper_area + lower_area + diagonal_len;
    return ret;
  }

  inline int stride() const { return 1; }

  inline size_type oned_length(size_type i) const {
    return MTL_MAX(0, MTL_MIN(int(dim.second()), int(i) + bw.second() + 1)
               - MTL_MAX(0, int(i) - bw.first()));
  }
  inline size_type twod_length() const { return dim.first(); }

  inline static size_type size(int m, int n, int low, int up) { 
    packed_offset offset(m, n, n, band_type(low, up));
    return offset.oned_offset(m);
  }
  inline size_type major() const { return dim.first(); }
  inline size_type minor() const { return dim.second(); }

private:
  dim_type dim;
  band_type bw; /* bandwidth */
};

//: blah
//!noindex:
template <int M, int N>
struct gen_packed_offset {
#if defined( _MSVCPP_ )
  typedef packed_offset<unsigned int, M, N> type;
#else
  template <class size_type>
  struct bind {
    typedef packed_offset<size_type, M, N> type;
  };
#endif
  typedef gen_banded_offset<M,N> transpose_type; /* bogus */
  typedef gen_banded_view_offset<M,N> banded_view_type; /* bogus */
};


/* egcs doesn't "see" the friend functions
 *   when the dense2D_iterator class is defined inside of dense2D
 */

//: blah
//!noindex:
template <int isConst, class T, class Offset, class InnerOneD, class OneD>
class dense2D_iterator {
public:
  typedef typename Offset::size_type size_type;
  typedef std::pair<size_type,size_type> pair_type;

  typedef typename IF<isConst, const T*,T*>::RET Iterator;
    
  typedef dense2D_iterator self;

  typedef int distance_type;
  typedef int difference_type;
    
  typedef std::random_access_iterator_tag iterator_category;

  typedef OneD*           pointer;
  typedef OneD            value_type;
  typedef OneD            reference;
  typedef difference_type Distance;
  typedef Iterator        iterator_type;
    
protected:
    
  Iterator start;
  size_type pos;
  size_type ld;    /* leading dimension */
  pair_type starts;
  Offset offset;
public:
    
  inline size_type index() const { return pos + starts.second; }
    
  inline dense2D_iterator () {}
    
  inline dense2D_iterator(const self& x)
    : start(x.start), pos(x.pos), 
      ld(x.ld), starts(x.starts), offset(x.offset) { }

  inline self& operator=(const self& x) {
    start = x.start;
    pos = x.pos;
    ld = x.ld;
    starts = x.starts;
    offset = x.offset;
    return *this;
  }

  inline explicit
  dense2D_iterator(Iterator x, size_type ld_, size_type p, pair_type s,
                   Offset os)
    : start(x), pos(p), ld(ld_), starts(s), offset(os) { }
    
  inline Iterator base () const { return start + pos; }
    
  inline reference deref(Distance pos, not_strided_tag) const { 
    return reference((T*)start + offset.oned_offset(pos), 
		offset.oned_length(pos), 
		starts.first);
  }
  inline reference deref(Distance pos, strided_tag) const { 
    InnerOneD vec((T*)start + offset.oned_offset(pos), 
                  offset.oned_length(pos), 
                  starts.first);
    return strided(vec, offset.stride());
  }
  inline reference operator*() const { 
    typedef typename Offset::is_strided Strided;
    return deref(pos, Strided());
  }
  inline reference operator[] (Distance n) const { 
    typedef typename Offset::is_strided Strided;
    return deref(pos + n, Strided());
  }

  /*  won't work, the OneD is temporary 
  pointer   operator-> () const { return & (operator* ()); }
  */
    
  inline self& operator++ () { ++pos; return *this; }
  inline self operator++ (int) { self tmp = *this; ++pos; return tmp; }
  inline self& operator-- () { --pos; return *this; }
  inline self operator-- (int) { self tmp = *this; --pos; return tmp; }
  inline self& operator+=(size_type n) { pos += n; return *this; }
  inline self operator+(size_type n) const {
    return self(start, ld, pos + n, starts); 
  }
  inline self& operator-=(size_type n) { pos -= n; return *this; }
    
    
    
};

template <int isConst, class T, class Offset, class InnerOneD, class OneD>
inline typename dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>::difference_type
operator-(const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& x,
          const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& y)
{
  return x.index() - y.index(); 
}
    
template <int isConst,class T, class Offset, class InnerOneD, class OneD>
inline bool
operator== (const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& x,
            const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& y)
{
  return x.index() == y.index();
}
    
template <int isConst, class T, class Offset, class InnerOneD, class OneD>
inline bool
operator!= (const dense2D_iterator<isConst, T,Offset,InnerOneD,OneD>& x,
            const dense2D_iterator<isConst, T,Offset,InnerOneD,OneD>& y)
{
  return x.index() != y.index();
}
    
template <int isConst,class T, class Offset, class InnerOneD, class OneD>
inline bool
operator< (const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& x,
           const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& y)
{
  return x.index() < y.index(); 
}



/*
  Workaround (g++ 2.91) helper class
 */

template <class Strided>
struct __bracket { };

template <>
struct __bracket<strided_tag> { 
  template <class OneD, class InnerOneD, class elt_type, class size_type>
  inline OneD
  operator()(elt_type* d, size_type len, size_type f, size_type ld,
             const OneD*, const InnerOneD*) {
    InnerOneD vec(d , len, f);
    return OneD(vec, ld);
  }
}; 
 
template <>
struct __bracket<not_strided_tag> { 
  template <class OneD, class InnerOneD, class elt_type, class size_type>
  inline OneD
  operator()(elt_type* d, size_type len, size_type f, size_type,
             const OneD*, const InnerOneD*) {
    return OneD(d, len, f);
  }
};



template<class T, class OffsetGen, int MM, int NN>
class dense2D;

template <class T, class OffsetGen, int MM, int NN>
class external2D;

//: Generic Dense 2-D Container
//!category: containers
//!component: type
//
// The generic_dense2D container implements sevaral of the MTL storage
// types.  They include dense, packed, banded, and banded_view.  The
// common theme here is that the matrix is stored in a contiguous
// piece of memory.  The differences in these storage types has to do
// with where to find the OneD segements in the linear
// memory. Caclulating these offsets is the job of the Offset concept,
// which has a model to handle each of the different storage types:
// rect_offset, strided_offset, banded_offset, packed_offset, and
// banded_view_offset.
//
// There are two derived classes of generic_dense2D that specify the
// memory management, dense2D and external2D. The dense2D version owns
// its memory, while the external2D imports its memory from somewhere
// else through a pointer (which allows for interoperability with
// other codes -- even with Fortran!).  <p>
//
//!definition: dense2D.h
//!tparam: RepType - The Container used to store the elements
//!tparam: RepPtr - The type used to reference to the container
//!tparam: OffsetGen - The generator that creates the Offset class
//!tparam: MM - For static sized matrix, the major dimension
//!tparam: NN - For static sized matrix, the minor dimension
//!models: TwoDStorage

template <class RepType, class RepPtr, class OffsetGen, int MM, int NN>
class generic_dense2D {
public:
  //: Static sizes (0 if dynamic)
  enum { M = MM, N = NN };

  //: The type for dimensions and indices
  typedef typename RepType::size_type size_type;
  //: The type for differences between iterators
  typedef typename RepType::difference_type difference_type;

protected:
  typedef std::pair<size_type,size_type> pair_type;
  typedef RepType reptype;
  typedef RepPtr rep_ptr;
  typedef typename RepType::value_type elt_type;

#if defined(_MSVCPP_)
  //JGS Nasty VC++ workaround
  typedef typename OffsetGen::type Offset;
#else
  typedef typename OffsetGen:: MTL_TEMPLATE bind<size_type>::type Offset;
#endif

  typedef dimension<elt_type> dyn_dim;
public:
  //: A pair type for dimensions
  typedef typename Offset::dim_type dim_type;

  //: A pair type for bandwidth
  typedef typename Offset::band_type band_type;

  /* Type Definitions */
  
  //: This is a dense matrix
  typedef dense_tag sparsity;

  typedef typename Offset::is_strided is_strided;

protected:
  typedef external_vec<elt_type, N> InnerOneD;

#if defined(_MSVCPP_)
  enum { offset_strided = Offset::IS_STRIDED };
  typedef typename IF<offset_strided, strided1D<InnerOneD>, InnerOneD>::RET OneD;
#else
  typedef typename Offset:: MTL_TEMPLATE bind_oned<InnerOneD>::type OneD;
#endif

  typedef OneD OneDRef;
  typedef OneD ConstOneDRef;
public:

  //: The 1D container type
  typedef OneD value_type;
  //: The type for a reference to value_type
  typedef value_type reference;
  //: The type for a const reference to value_type
  typedef value_type const_reference;

  //: The iterator type
  typedef dense2D_iterator<0,elt_type, Offset, InnerOneD, OneD> iterator;

  //: The const iterator type
  typedef dense2D_iterator<1,elt_type, Offset, InnerOneD, OneD> const_iterator;

  //: The reverse iterator type
  typedef reverse_iter<iterator> reverse_iterator;

  //: The const reverse iterator type
  typedef reverse_iter<const_iterator> const_reverse_iterator;

  //: The type for the transpose of this container
  typedef generic_dense2D<RepType, RepPtr,
             typename OffsetGen::transpose_type, MM, NN> transpose_type;

  //: The type for a banded view of this container
  typedef generic_dense2D<RepType, RepPtr,
             typename OffsetGen::banded_view_type, MM, NN> banded_view_type;

  //: The type for a sub-section of this 2D container
  typedef external2D<elt_type, OffsetGen, MM, NN> submatrix_type;

#ifndef MTL_DISABLE_BLOCKING
  template <class Block>
  struct blocked_view {
    typedef block2D<Block, OffsetGen> type;
  };
#endif

  //: This is a stridable container, can use rows(A), columns(A)
  typedef typename Offset::strideability strideability;

  /* Constructors */
  
  //: Default Constructor
  inline generic_dense2D()
    : ld_(0), data_(0), starts(std::make_pair(0,0)) { }

  //: Normal Constructor
  inline generic_dense2D(rep_ptr data, size_type m, size_type n, size_type ld)
    : ld_(ld), data_(data),
      starts(std::make_pair(0,0)), offset(m, n, ld) { }

  //: Constructor with non-zero upper-left corner indices
  inline generic_dense2D(rep_ptr data, size_type m, size_type n, 
			 size_type ld, dyn_dim s, char)
    : ld_(ld), data_(data),
      starts(std::make_pair(s.first(),s.second())), offset(m, n, ld) { }

  //: Static M, N constructor
  inline generic_dense2D(rep_ptr data, size_type ld)
    : ld_(ld), data_(data),
      starts(std::make_pair(0,0)), offset(M, N, ld) { }

  //: with bandwidth constructor
  inline generic_dense2D(rep_ptr data, size_type m, size_type n, size_type ld,
                         band_type bw)
    : ld_(ld), data_(data), starts(std::make_pair(0,0)),
      offset(m, n, ld, bw) { }

  //: Static M, N with bandwith?

  //: Copy Constructor
  inline generic_dense2D(const generic_dense2D& x)
    : ld_(x.ld_), data_(x.data_), starts(x.starts),
      offset(x.offset) { }

  //: Assignment Operator
  inline generic_dense2D& operator=(const generic_dense2D& x) {
    ld_ = x.ld_; data_ = x.data_; starts = x.starts; offset = x.offset;
    return *this;
  }

  //: Subclass Constructor
  inline generic_dense2D(rep_ptr d, const generic_dense2D& x)
    : ld_(x.ld_), data_(d), starts(x.starts), offset(x.offset) { }

  //: Transpose Constructor
  inline generic_dense2D(const transpose_type& x, do_transpose, do_transpose)
    : ld_(x.ld_), data_(x.data_), starts(x.starts), offset(x.offset) { }      

  /* JGS, remove stream constructor, not very necessary
     just have them call another constructor 
   */

  //: Matrix Stream Constructor
  template <class MatrixStream, class Orien>
  inline generic_dense2D(rep_ptr data, MatrixStream& s, Orien)
    : ld_(Orien::map(dim_type(s.nrows(),s.ncols())).second()),
      data_(data),
      starts(std::make_pair(0,0)),
      offset(Orien::map(dim_type(s.nrows(),s.ncols())).first(), 
             Orien::map(dim_type(s.nrows(),s.ncols())).second(),
             Orien::map(dim_type(s.nrows(),s.ncols())).second()) { }
  //: Banded Matrix Stream Constructor
  template <class MatrixStream, class Orien>
  inline generic_dense2D(rep_ptr data, MatrixStream& s, 
                         Orien, band_type bw)
    : ld_(Orien::map(dim_type(s.nrows(),s.ncols())).second()),
      data_(data), 
      starts(std::make_pair(0,0)),
      offset(Orien::map(dim_type(s.nrows(),s.ncols())).first(), 
             Orien::map(dim_type(s.nrows(),s.ncols())).second(),
             Orien::map(dim_type(s.nrows(),s.ncols())).second(),
             bw) { }

  //: Banded View Constructor
  template <class TwoD>
  inline generic_dense2D(rep_ptr data, const TwoD& x, band_type bw, banded_tag)
    : ld_(x.ld_),
      data_(data),
      starts(x.starts),
      offset(x.offset, bw) { }

// VC++ doesn't like this
  //friend class transpose_type;

  //: The destructor.