dense2d.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_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"
#endif

namespace 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: type
template <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;
  inline rect_offset() : dim(4444,4444), ld(4444) { }
  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: type
template <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: type
template <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: type
template <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: type
template <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: type
template <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 };
  typedef size_t size_type;