typedef size_t size_type;  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_typeoperator-(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 booloperator== (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 booloperator!= (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 booloperator< (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: TwoDStoragetemplate <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.  inline ~generic_dense2D() { }