// -*- 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_DENSE2D_H#define MTL_DENSE2D_H#include "mtl/mtl_iterator.h"#include <utility>#include <assert.h>#include <vector>#include "mtl/mtl_config.h"#include "mtl/linalg_vec.h"#include "mtl/strided1D.h"#include "mtl/initialize.h"#include "mtl/reverse_iter.h"#include "mtl/matrix_traits.h"#include "mtl/dimension.h"#ifndef MTL_DISABLE_BLOCKING#include "mtl/block2D.h"#endifnamespace mtl {template <class size_t, int MM, int NN>class strided_offset;struct strided_tag { enum { id = 1 }; };struct not_strided_tag { enum { id = 0 }; };template <class size_t, int MM, int NN>class band_view_offset;template <class size_t, int MM, int NN>class strided_band_view_offset;template <int M, int N> struct gen_rect_offset;template <int M, int N> struct gen_strided_offset;template <int M, int N> struct gen_banded_offset;template <int M, int N> struct gen_banded_view_offset;template <int M, int N> struct gen_strided_band_view_offset;template <int M, int N> struct gen_packed_offset;//: Rectangular Offset Class//!models: Offset//!category: utilities//!component: typetemplate <class size_t, int MM, int NN>class rect_offset {public:#if !defined(_MSVCPP_)  template <class Vec>  struct bind_oned {     typedef Vec type;   };#endif  typedef not_strided_tag is_strided;  typedef size_t size_type;  enum { M = MM, N = NN, IS_STRIDED = 0 };  typedef dimension<size_type, MM, NN> dim_type;  typedef dimension<int> band_type;  typedef strided_offset<size_type, MM, NN> transpose_type;  typedef strideable strideability;  // VC++ doesn't like this   //friend class transpose_type;  //what is that for? -- llee  //inline rect_offset() : dim(4444,4444), ld(4444) { }  inline rect_offset() : dim(0,0), ld(0) { }  inline rect_offset(const rect_offset& x) : dim(x.dim), ld(x.ld) { }  inline rect_offset(size_type m, size_type n, size_type ld_)    : dim(m, n), ld(ld_) { }  inline rect_offset(size_type m, size_type n, size_type ld_, band_type)    : dim(m, n), ld(ld_) { }  rect_offset(const transpose_type& x); /* see below strided_offset for def */  inline rect_offset& operator=(const rect_offset& x) {    dim = x.dim; ld = x.ld; return *this;  }  inline size_type elt(size_type i, size_type j) const { return i * ld + j; }  inline size_type oned_offset(size_type i) const { return i * ld; }  inline size_type oned_length(size_type) const { return dim.second(); }  inline size_type twod_length() const { return dim.first(); }  inline size_type stride() const { return 1; }  inline static size_type size(size_type m, size_type n,                               size_type , size_type) { return m * n; }  inline size_type major() const { return dim.first(); }  inline size_type minor() const { return dim.second(); }  /* private: */  dim_type dim;  size_type ld;};//: blah//!noindex:template <int M, int N>struct gen_rect_offset {#if defined( _MSVCPP_ )   typedef rect_offset<unsigned int, M, N> type;#else  template <class size_type>  struct bind {    typedef rect_offset<size_type, M, N> type;  };#endif  typedef gen_strided_offset<M,N> transpose_type;  typedef gen_banded_view_offset<M,N> banded_view_type;};//: Strided Rectangular Offset Class//!models: Offset//!category: utilities//!component: typetemplate <class size_t, int MM, int NN>class strided_offset {public:#if !defined(_MSVCPP_)  template <class Vec>  struct bind_oned {     typedef strided1D<Vec> type;  };#endif  /*  typedef strided_band_view_offset<size_t,MM,NN> banded_view_type;   */  typedef strided_tag is_strided;  enum { M = MM, N = NN, IS_STRIDED = 1 };  typedef size_t size_type;  typedef dimension<size_type, MM, NN> dim_type;  typedef dimension<int> band_type;  typedef rect_offset<size_type,MM,NN> transpose_type;  typedef strideable strideability;// VC++ doesn't like this  //friend class transpose_type;  inline strided_offset() : dim(0,0), ld(0) { }  inline strided_offset(size_type m, size_type n, size_type ld_)    : dim(m, n), ld(ld_) { }  inline strided_offset(const transpose_type& x) : dim(x.dim), ld(x.ld) { }  inline strided_offset& operator=(const strided_offset& x) {    dim = x.dim; ld = x.ld; return *this;  }  inline size_type elt(size_type i, size_type j) const { return j * ld + i; }  inline size_type oned_offset(size_type i) const { return i; }  inline size_type oned_length(size_type) const { return dim.first() * ld; }  inline size_type twod_length() const { return dim.second(); }  inline size_type stride() const { return ld; }  inline static size_type size(size_type m, size_type n,                               size_type , size_type) { return m * n; }  inline size_type major() const { return dim.first(); }  inline size_type minor() const { return dim.second(); }  /* private: */  dim_type dim;  size_type ld;};//: blah//!noindex:template <int M, int N>struct gen_strided_offset {#if defined( _MSVCPP_ )  typedef strided_offset<unsigned int, M, N> type;#else  template <class size_type>  struct bind {    typedef strided_offset<size_type, M, N> type;  };#endif  typedef gen_rect_offset<M,N> transpose_type;  typedef gen_strided_band_view_offset<M,N> banded_view_type;};template <class size_t, int MM, int NN>inline rect_offset<size_t,MM,NN>::rect_offset(const rect_offset<size_t,MM,NN>::transpose_type& x)  : dim(x.dim), ld(x.ld) { }//: Banded View Offset Class// This creates a banded view into a full matrix.//!models: Offset//!category: utilities//!component: typetemplate <class size_t, int MM, int NN>class banded_view_offset {public:#if !defined(_MSVCPP_)  template <class Vec>  struct bind_oned {     typedef Vec type;   };#endif  typedef not_strided_tag is_strided;  enum { M = MM, N = NN, IS_STRIDED = 0 };  typedef size_t size_type;  typedef dimension<size_type, MM, NN> dim_type;  typedef dimension<int, MM, NN> band_type;  typedef strided_band_view_offset<size_type, MM, NN> transpose_type;  typedef not_strideable strideability;// VC++ doesn't like this  //friend class transpose_type;  inline banded_view_offset()    : dim(0,0), ld(0), bw(std::make_pair(0,0)) { }  inline banded_view_offset(size_type m, size_type n, size_type leading_dim,                            band_type band)      : dim(m, n), ld(leading_dim), bw(band) { }  inline banded_view_offset(size_type m, size_type n, size_type leading_dim)      : dim(m, n), ld(leading_dim), bw(band_type(0,0)) { }  template <class Offset>  inline banded_view_offset(Offset os, band_type band)      : dim(os.dim), ld(os.ld), bw(band) { }  inline banded_view_offset& operator=(const banded_view_offset& x) {    dim = x.dim; ld = x.ld; bw = x.bw; return *this;  }  inline size_type elt(size_type i, size_type j) const {     size_type start = MTL_MAX(int(i) - bw.first(), 0);    return i * ld + j + start;   }  inline size_type oned_offset(size_type i) const {     size_type start = MTL_MAX(int(i) - bw.first(), 0);    return i * ld + start;   }  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 int stride() const { return 1; }  inline static size_type size(size_type m, size_type n,                                size_type , size_type) {    return m * n;   }  inline size_type major() const { return dim.first(); }  inline size_type minor() const { return dim.second(); }    /* private: */  dim_type dim;  size_type ld;  band_type bw; /* bandwidth */};//: blah//!noindex:template <int M, int N>struct gen_banded_view_offset {#if defined( _MSVCPP_ )  typedef banded_view_offset<unsigned int, M, N> type;#else  template <class size_type>  struct bind {    typedef banded_view_offset<size_type, M, N> type;  };#endif  typedef gen_strided_band_view_offset<M,N> transpose_type;  typedef gen_banded_view_offset<M,N> banded_view_type; /* bogus */};//: Strided Band View Offset Class//// This creates a strided band view into a full matrix.// This class is to banded_view as strided_offset is to rect_offset.////!models: Offset//!category: utilities//!component: typetemplate <class size_t, int MM, int NN>class strided_band_view_offset {public:#if !defined(_MSVCPP_)  template <class Vec>  struct bind_oned {     typedef strided1D<Vec> type;   };#endif  typedef strided_tag is_strided;  enum { M = MM, N = NN, IS_STRIDED = 1 };  typedef size_t size_type;  typedef dimension<size_type, MM, NN> dim_type;  typedef dimension<int, MM, NN> band_type;  typedef banded_view_offset<size_type, MM, NN> transpose_type;  typedef not_strideable strideability;// VC++ doesn't like this  //friend class transpose_type;  inline strided_band_view_offset()    : dim(0,0), ld(0), bw(std::make_pair(0,0)) { }  inline strided_band_view_offset(size_type m, size_type n, 				  size_type leading_dim,				  band_type band)      : dim(m, n), ld(leading_dim), bw(band) { }  template <class Offset>  inline strided_band_view_offset(Offset os, band_type band)      : dim(os.dim), ld(os.ld), bw(band) { }  inline strided_band_view_offset&   operator=(const strided_band_view_offset& x) {    dim = x.dim; ld = x.ld; bw = x.bw; return *this;  }  inline size_type elt(size_type i, size_type j) const {     size_type start = MTL_MAX(int(i) - bw.first(), 0);    return (j + start) * ld + i;   }  inline size_type oned_offset(size_type i) const {     size_type start = MTL_MAX(int(i) - bw.first(), 0);    return start * ld + i;   }  inline size_type oned_length(size_type i) const {     /* use dim.first() here */    size_type len =  MTL_MAX(0, MTL_MIN(int(dim.first()), int(i) + bw.second() + 1)                          - MTL_MAX(0, int(i) - bw.first()));    return len * ld;  }  inline size_type twod_length() const { return dim.second(); }  inline int stride() const { return ld; }  inline static size_type size(size_type m, size_type n,                                size_type , size_type) {    return m * n;   }  inline size_type major() const { return dim.first(); }  inline size_type minor() const { return dim.second(); }    /* private: */  dim_type dim;  size_type ld;  band_type bw; /* bandwidth */};//: blah//!noindex:template <int M, int N>struct gen_strided_band_view_offset {#if defined( _MSVCPP_ )  typedef strided_band_view_offset<unsigned int, M, N> type;#else  template <class size_type>  struct bind {    typedef strided_band_view_offset<size_type, M, N> type;  };#endif  typedef gen_banded_view_offset<M,N> transpose_type;  typedef gen_strided_band_view_offset<M,N> banded_view_type; /* bogus */};template <class size_t, int MM, int NN>class packed_offset;//: Banded Offset Class// This cooresponds to lapack/blas banded storage format.//!models: Offset//!category: utilities//!component: typetemplate <class size_t, int MM, int NN>class banded_offset {public:#if !defined(_MSVCPP_)  template <class Vec>  struct bind_oned {     typedef Vec type;   };#endif  typedef not_strided_tag is_strided;  enum { M = MM, N = NN, IS_STRIDED = 0 };  typedef size_t size_type;  typedef dimension<size_type, MM, NN> dim_type;  typedef dimension<int> band_type;  typedef packed_offset<size_type,MM,NN> transpose_type; /* bogus */  typedef not_strideable strideability;  inline banded_offset()    : dim(0,0), bw(band_type(0,0)), ndiag(0) { }  inline banded_offset(size_type m, size_type n, size_type /* lead */,                        band_type band)     : dim(m,n), bw(band), ndiag(band.first() + band.second() + 1) { } inline banded_offset(size_type m, size_type n, size_type /* lead */)     : dim(m,n), bw(band_type(0,0)), ndiag(0) { }  inline banded_offset& operator=(const banded_offset& x) {    dim = x.dim; ndiag = x.ndiag; bw = x.bw; return *this;  }  inline size_type elt(size_type i, size_type j) const {    return this->oned_offset(i) + j;   }  inline size_type oned_offset(size_type i) const {    return i * ndiag + MTL_MAX(0, bw.first() - int(i));  }  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 int stride() const { return 1; }  inline static size_type size(size_type m, size_type n,                                size_type low, size_type up) {     /* M' = number of diagonals = low + up + 1       N' = min (m, n + low) */    return (low + up + 1) * MTL_MIN(m, n + low);  }  inline size_type major() const { return dim.first(); }  inline size_type minor() const { return dim.second(); }private:  dim_type dim;  band_type bw; /* bandwidth */  size_type ndiag;};//: blah//!noindex:template <int M, int N>struct gen_banded_offset {#if defined( _MSVCPP_ )  typedef banded_offset<unsigned int, M, N> type;#else  template <class size_type>  struct bind {    typedef banded_offset<size_type, M, N> type;  };#endif  typedef gen_packed_offset<M,N> transpose_type; // bogus  typedef gen_banded_view_offset<M,N> banded_view_type; /* bogus */};//: Packed Offset Class// This cooresponds to lapack/blas packed storage format //!models: Offset//!category: utilities//!component: typetemplate <class size_t, int MM, int NN>class packed_offset {public:#if !defined(_MSVCPP_)  template <class Vec>  struct bind_oned {     typedef Vec type;  };#endif  typedef not_strided_tag is_strided;  enum { M = MM, N = NN, IS_STRIDED = 0 };