matrix_implementation.h

强大的矩阵模版类

  inline const Indexer& get_indexer() const { return indexer; }

  inline void print() const { twod.print(); }

  /* external storage interface */
  inline value_type* data() { return twod.data(); }
  inline const value_type* data() const { return twod.data(); }

  /* compressed2D external storage interface */
  inline value_type* get_val() { return twod.get_val(); }
  inline const value_type* get_val() const { return twod.get_val(); }
  inline size_type* get_ind() { return twod.get_ind(); }
  inline const size_type* get_ind() const { return twod.get_ind(); }
  inline size_type* get_ptr() { return twod.get_ptr(); }
  inline const size_type* get_ptr() const { return twod.get_ptr(); }

  template <class Matrix>
  inline void fast_copy(const Matrix& x) { twod.fast_copy(x); }


  /* protected: */
  TwoD twod;
  Indexer indexer;
};


template <class TwoDGen, class IndexerGen>
class column_matrix;


//: row matrix
//
// This class derives from the matrix_implementation class.
// The main purpose of this class is merely to add the "Row"
// type definition.

template <class TwoDGen, class IndexerGen>
class row_matrix : public matrix_implementation<TwoDGen, IndexerGen> {
  typedef matrix_implementation<TwoDGen, IndexerGen> Base;
public:
  typedef typename Base::value_type value_type;
  typedef typename Base::size_type size_type;
  typedef typename Base::pointer pointer;
  typedef typename Base::reference reference;
  typedef typename Base::const_reference const_reference;
  typedef typename Base::dim_type dim_type;
  typedef typename Base::dyn_dim dyn_dim;
  typedef typename Base::TwoD TwoD;
  typedef typename Base::OneD OneD;
  typedef typename Base::OneDRef OneDRef;
  typedef typename Base::Indexer Indexer;
  typedef OneD Row;
  typedef OneDRef RowRef;
  enum { M = Indexer::M, N = Indexer::N };
  typedef column_matrix<TwoDGen, typename IndexerGen::transpose_type>
             transpose_type;
  typedef column_matrix<typename TwoDGen::transpose_type,
                        typename IndexerGen::strided_type> strided_type;

  typedef row_matrix<typename TwoDGen::banded_view_type,
        gen_banded_indexer<typename IndexerGen::orienter,M,N,size_type> > banded_view_type;


  //#if !defined( __GNUC__ ) && !defined( _MSVCPP_ )  /* internal compiler error */
#if !defined( MTL_DISABLE_BLOCKING )
  template <int BM, int BN>
  struct blocked_view {
#if 0
    typedef matrix<value_type, rectangle<BM,BN>,
                   dense<external>, row_major>::type Block;
#else
    typedef typename TwoD::is_strided IsStrided;

    //    typedef light_matrix<value_type, size_type, ROW_MAJOR, IsStrided::id> Block;
    typedef light_matrix<value_type, int, ROW_MAJOR, IsStrided::id> Block;
    //    typedef light_matrix<value_type, 1, TwoD::is_strided::id,BM,BN> Block;
#endif

    typedef typename TwoDGen:: MTL_TEMPLATE blocked_view<Block>::type BlockTwoDGen;
    typedef row_matrix<BlockTwoDGen, IndexerGen> type;
  };
#endif

  inline row_matrix() { }

  inline row_matrix& operator=(const row_matrix& x) {
    Base::operator=(x);
    return *this;
  }
  //: user callable constructors
  inline row_matrix(size_type m, size_type n) 
    : Base(dim_type(m, n)) { }
  inline row_matrix(size_type m, size_type n, size_type nnz_max) 
    : Base(dim_type(m, n), nnz_max) { }
  inline row_matrix(size_type m, size_type n, int sub, int super) 
    : Base(dim_type(m, n), band_type(sub, super)) { }

  //: external data constructors
  inline row_matrix(pointer data, size_type m, size_type n)
    : Base(data, dim_type(m, n), char()) { }
  inline row_matrix(pointer data, size_type m, size_type n, size_type ld)
    : Base(data, dim_type(m, n), ld) { }

  //: non zero index upper left corner
  inline row_matrix(pointer data, size_type m, size_type n, 
                    size_type ld, dyn_dim starts)
    : Base(data, dim_type(m, n), ld, starts, char()) { }

  inline row_matrix(pointer data, size_type m, size_type n,
                    int sub, int super)
    : Base(data, dim_type(m, n), 
           band_type(sub, super)) { }
  inline row_matrix(pointer data, size_type m, size_type n, size_type ld, 
                    int sub, int super)
    : Base(data, dim_type(m, n), ld,
           band_type(sub, super)) { }
  //: Static M, N Constructor
  inline row_matrix(pointer data)
    : Base(data) { }

  //: compressed2D external data constructor
  inline row_matrix(size_type m, size_type n, size_type nnz,
             pointer val, size_type* ptrs, size_type* inds)
    : Base(dim_type(m, n), nnz, val, ptrs, inds) { }

  //: banded view constructor
  template <class Matrix>
  inline row_matrix(int sub, int super, const Matrix& x)
    : Base(x, band_type(sub, super)) { }

  //: block view constructor
  template <class Matrix, class ST, int BM, int BN>
  inline row_matrix(const Matrix& x, mtl::dimension<ST,BM,BN> bdim)
    : Base(x, bdim, char()) { }

  //: stream constructors
  typedef matrix_market_stream<value_type> mmstream;
  typedef harwell_boeing_stream<value_type> hbstream;
  inline row_matrix(mmstream & m_in) : Base(m_in, *this) { }
  inline row_matrix(hbstream & m_in) : Base(m_in, *this) { }
  inline row_matrix(mmstream& m_in, int sub, int super)
    : Base(m_in, band_type(sub, super), *this) { }
  inline row_matrix(hbstream& m_in, int sub, int super)
    : Base(m_in, band_type(sub, super), *this) { }

  template <class Subclass>
  inline row_matrix(mmstream & m_in, Subclass& s) : Base(m_in, s) { }
  template <class Subclass>
  inline row_matrix(hbstream & m_in, Subclass& s) : Base(m_in, s) { }
  template <class Subclass>
  inline row_matrix(mmstream& m_in, int sub, int super, Subclass& s)
    : Base(m_in, band_type(sub, super), s) { }
  template <class Subclass>
  inline row_matrix(hbstream& m_in, int sub, int super, Subclass& s)
    : Base(m_in, band_type(sub, super), s) { }

  //: copy constructor
  inline row_matrix(const row_matrix& x)
    : Base(x) { }

  //:: called by rows(A) and columns(A)
  inline row_matrix(const row_matrix& x, do_strided)
    : Base(x) { }

  //: called by trans helper function
  inline row_matrix(const transpose_type& x, do_transpose t)
    : Base(x, t, t) { }

  //: called by rows and columns helper functions
  inline row_matrix(const strided_type& x, do_strided s)
    : Base(x, s, s) { }

  //: called by scaled helper function
  template <class MatrixT, class ScalarT>
  inline row_matrix(const MatrixT& x, const ScalarT& y, do_scaled s) 
    : Base(x, y, s) { }

  inline row_matrix(const TwoD& x, Indexer ind)
    : Base(x, ind) { }

  inline ~row_matrix() { }

  inline void resize(size_type m, size_type n) { 
    twod.resize(m, n); 
    indexer.dim = dim_type(m, n);
  }

  /* submatrix */
#if 1
  typedef typename TwoDGen::submatrix_type SubTwoDGen;  
  typedef typename SubTwoDGen::type::size_type new_sizeT;
#if defined(_MSVCPP_)
  typedef IndexerGen SubMatIndexerGen;
#else
  typedef typename IndexerGen:: MTL_TEMPLATE bind<new_sizeT>::other SubMatIndexerGen;
#endif
  typedef row_matrix<SubTwoDGen, SubMatIndexerGen> submatrix_type;

  // problem, The new TwoDGen::submatrix_type may have a difference
  // size_type then previously. Since IndexerGen still has the
  // old size_type, there is a conflict
#else
  //  typedef light_matrix<value_type, 1, TwoD::is_strided::id, 0, 0> submatrix_type;
  typedef typename TwoD::is_strided IsStrided;
  //correct  typedef light_matrix<value_type, size_type, ROW_MAJOR, IsStrided::id> submatrix_type;
  typedef light_matrix<value_type, int, ROW_MAJOR, IsStrided::id> submatrix_type;
#endif

  //: create a submatrix view into this matrix
  // Currently only works with dense rectangular matrices
  inline submatrix_type sub_matrix(size_type row_start, size_type row_finish,
                                   size_type col_start, size_type col_finish) const
  {
    dim_type starts = indexer.at(dim_type(row_start, col_start));
    size_type m = row_finish - row_start;
    size_type n = col_finish - col_start;
    typedef typename TwoD::is_strided IsStrided;
    if (IsStrided::id) {
      return submatrix_type((value_type*)twod.data() 
			    + starts.second() * twod.ld() + starts.first(),
			    m, n, twod.ld());
    } else {
      return submatrix_type((value_type*)twod.data() 
			    + starts.first() * twod.ld() + starts.second(),
			    m, n, twod.ld());
    }
  }

#if 0
  // deprecated, derive partitioned matrix type from submatrix_type instead
  // and do it in partition.h
  /* partition */
  typedef row_matrix< typename TwoD::template partitioned<submatrix_type>::generator,
                      IndexerGen> partitioned;
#endif
  
#if 0 // deprecated, use functions in partition.h instead
  template <class Sequence1, class Sequence2>
  inline partitioned
  partition(const Sequence1& prows, const Sequence2& pcols) const
  {
    return partition_matrix(prows, pcols, *this);
  }
  inline partitioned
  subdivide(size_type split_row, size_type split_col) const
  {
    dense1D<size_type> row_splits(1);
    dense1D<size_type> col_splits(1);

    row_splits[0] = split_row;
    col_splits[0] = split_col;

    return partition_matrix(row_splits, col_splits, *this);
  }
#endif
};

//: column matrix
//
// This class derives from the matrix_implementation class.
// The main purpose of this class is merely to add the "Column"
// type definition.

template <class TwoDGen, class IndexerGen>
class column_matrix : public matrix_implementation<TwoDGen, IndexerGen> {
  typedef matrix_implementation<TwoDGen, IndexerGen> Base;
public:
  typedef typename Base::pointer pointer;
  typedef typename Base::reference reference;
  typedef typename Base::const_reference const_reference;
  typedef typename Base::value_type value_type;
  typedef typename Base::size_type size_type;
  typedef typename Base::dim_type dim_type;
  typedef typename Base::dyn_dim dyn_dim;
  typedef typename Base::band_type band_type;
  typedef typename Base::TwoD TwoD;
  typedef typename Base::OneD OneD;
  typedef typename Base::OneDRef OneDRef;
  typedef typename Base::Indexer Indexer;
  typedef OneD Column;
  typedef OneDRef ColumnRef;
  enum { M = Indexer::M, N = Indexer::N };
  typedef row_matrix<TwoDGen,
                     typename IndexerGen::transpose_type> transpose_type;
  typedef row_matrix<typename TwoDGen::transpose_type,
                     typename IndexerGen::strided_type> strided_type;

  typedef column_matrix<typename TwoDGen::banded_view_type,
        gen_banded_indexer<typename IndexerGen::orienter,M,N,size_type> > banded_view_type;

  //#if !defined( __GNUC__) && !defined( _MSVCPP_ ) /* internal compiler error */
#if !defined( MTL_DISABLE_BLOCKING )
  template <int BM, int BN>
  struct blocked_view {
#if 0
    typedef matrix<value_type, rectangle<BM,BN>, 
                   dense<external>, column_major>::type Block;
#else
    typedef typename TwoD::is_strided IsStrided;
    typedef light_matrix<value_type, size_type, COL_MAJOR, IsStrided::id> Block;
    //    typedef light_matrix<value_type, 0, TwoD::is_strided::id,BM,BN> Block;
#endif
    typedef column_matrix<typename TwoDGen:: MTL_TEMPLATE blocked_view<Block>::type, 
                          IndexerGen> type;
  };
#endif

  inline column_matrix() { }

  //: user callable constructors
  inline column_matrix(size_type m, size_type n) 
    : Base(dim_type(m, n)) { }
  inline column_matrix(size_type m, size_type n, size_type nnz_max) 
    : Base(dim_type(m, n), nnz_max) { }
  inline column_matrix(size_type m, size_type n, int sub, int super) 
    : Base(dim_type(m, n), 
           band_type(sub, super)) { }

  inline column_matrix& operator=(const column_matrix& x) {
    Base::operator=(x);
    return *this;
  }
  //: external data
  inline column_matrix(pointer data, size_type m, size_type n)
    : Base(data, dim_type(m, n), char()) { }
  inline column_matrix(pointer data, size_type m, size_type n, size_type ld)
    : Base(data, dim_type(m, n), ld) { }

  //: non zero index upper left corner
  inline column_matrix(pointer data, size_type m, size_type n, 
                       size_type ld, dyn_dim starts)
    : Base(data, dim_type(m, n), ld, starts, char()) { }

  inline column_matrix(pointer data, size_type m, size_type n,
                       int sub, int super)
    : Base(data, dim_type(m, n), 
           band_type(sub, super)) { }
  inline column_matrix(pointer data, size_type m, size_type n, size_type ld, 
                       int sub, int super)
    : Base(data, dim_type(m, n), ld,
           band_type(sub, super)) { }

  //: Static M, N Constructor
  inline column_matrix(pointer data)
    : Base(data) { }

  //: banded view constructor
  template <class Matrix>
  inline column_matrix(int sub, int super, const Matrix& x) 
    : Base(x, band_type(sub, super)) { }

  //: block view constructor
  template <class Matrix, class ST, int BM, int BN>
  inline column_matrix(const Matrix& x, mtl::dimension<ST,BM,BN> bdim)
    : Base(x, bdim, char()) { }


  //: compressed2D external data constructor
  inline column_matrix(size_type m, size_type n, size_type nnz,
                pointer val, size_type* ptrs, size_type* inds)
    : Base(dim_type(m, n), nnz, val, ptrs, inds) { }

  //: streams
  typedef matrix_market_stream<value_type> mmstream;
  typedef harwell_boeing_stream<value_type> hbstream;
  inline column_matrix(mmstream& m_in) : Base(m_in, *this) { }
  inline column_matrix(mmstream& m_in, int sub, int super)
    : Base(m_in, band_type(sub,super), *this) { }
  inline column_matrix(hbstream& m_in) : Base(m_in, *this) { }
  inline column_matrix(hbstream& m_in, int sub, int super)
    : Base(m_in, band_type(sub,super), *this) { }

  template <class Subclass>
  inline column_matrix(mmstream & m_in, Subclass& s) : Base(m_in, s) { }
  template <class Subclass>
  inline column_matrix(hbstream & m_in, Subclass& s) : Base(m_in, s) { }
  template <class Subclass>
  inline column_matrix(mmstream& m_in, int sub, int super, Subclass& s)
    : Base(m_in, band_type(sub, super), s) { }
  template <class Subclass>
  inline column_matrix(hbstream& m_in, int sub, int super, Subclass& s)
    : Base(m_in, band_type(sub, super), s) { }

  //: copy constructor
  inline column_matrix(const column_matrix& x)
    : Base(x) { }

  //: called by rows(A), columns(A)
  inline column_matrix(const column_matrix& x, do_strided)
    : Base(x) { }

  //: called by trans(A) helper function
  inline column_matrix(const transpose_type& x, do_transpose t)
    : Base(x, t, t) { }

  //: called by rows(A) and columns(A) helper functions
  inline column_matrix(const strided_type& x, do_strided s)
    : Base(x, s, s) { }

  //: scaled
  template <class MatrixT, class ScalarT>
  inline column_matrix(const MatrixT& x, const ScalarT& y, do_scaled s) 
    : Base(x, y, s) { }

  inline column_matrix(const TwoD& x, Indexer ind)
    : Base(x, ind) { }

  inline ~column_matrix() { }

  inline void resize(size_type m, size_type n) { 
    twod.resize(n, m); 
    indexer.dim = dim_type(n, m);
  }

  /* submatrix */
#if 1
  typedef typename TwoDGen::submatrix_type SubTwoDGen;  
  typedef typename SubTwoDGen::type::size_type new_sizeT;
#if defined(_MSVCPP_)
  typedef IndexerGen SubMatIndexerGen;
#else
  typedef typename IndexerGen:: MTL_TEMPLATE bind<new_sizeT>::other SubMatIndexerGen;
#endif
  typedef column_matrix<SubTwoDGen, SubMatIndexerGen> submatrix_type;
                       
#else
  //  typedef light_matrix<value_type, 0, TwoD::is_strided::id, 0, 0> submatrix_type;
  typedef typename TwoD::is_strided IsStrided;
  typedef light_matrix<value_type, size_type, COL_MAJOR, IsStrided::id> submatrix_type;
#endif

  //: create a submatrix view into this matrix
  // Currently only works with dense rectangular matrices
  inline submatrix_type sub_matrix(size_type row_start, size_type row_finish,
                             size_type col_start, size_type col_finish) const
  {
    dim_type starts = indexer.at(dim_type(row_start, col_start));
    size_type m = row_finish - row_start;
    size_type n = col_finish - col_start;
    typedef typename TwoD::is_strided IsStrided;
    if (IsStrided::id) {
      return submatrix_type((value_type*)twod.data() 
			    + starts.second() * twod.ld() + starts.first(),