matrix_implementation.h

强大的矩阵模版类

			    m, n, twod.ld());
    } else {
      return submatrix_type((value_type*)twod.data() 
			    + starts.first() * twod.ld() + starts.second(),
			    m, n, twod.ld());
    }
  }

#if 0
  // JGS see row_matrix::partitioned
  /* partition */
  typedef column_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

};

template <class Matrix>
struct rows_type {
  typedef typename Matrix::orientation orien;
  enum { orienid = orien::id }; // VC++ workaround
  typedef typename IF<EQUAL<orienid,ROW_MAJOR>::RET,
                      Matrix, typename Matrix::strided_type>::RET type;
};

template <class Matrix>
struct columns_type {
  typedef typename Matrix::orientation orien;
  enum { orienid = orien::id }; // VC++ workaround
  typedef typename IF<EQUAL<orienid,COL_MAJOR>::RET,
                       Matrix, typename Matrix::strided_type>::RET type;
};

//: Access the row-wise view of the matrix
//
//  For matrix A, A[i] now gives you the ith row
//  and A.begin() gives you an iterator over rows
//
//!example: swap_rows.cc
//!component: function
//!category: containers
//!tparam: Matrix - The Matrix to access row-wise. Matrix must be dense.
template<class Matrix>
inline typename rows_type<Matrix>::type
rows(const Matrix& A) { 
  return rows_type<Matrix>::type(A, do_strided());
}

//: Access the column-wise view of the matrix
//
//  For matrix A, A[i] now gives you the ith column and A.begin()
//  gives you an iterator over columns. See rows for an example.
//  
//
//!component: function
//!category: containers
//!tparam: Matrix - The Matrix to access column-wise. Matrix must be dense.
template<class Matrix>
inline typename columns_type<Matrix>::type
columns(const Matrix& A) { 
  return columns_type<Matrix>::type(A, do_strided());
}

//: Swap the orientation of a matrix.
//
// Swap the orientation of a matrix (i.e., from row-major to
// column-major). In essence this transposes the matrix. This
// operation occurs at compile time.
//
//!component: function
//!category: containers
//!tparam: Matrix - The Matrix to transpose
//!example: trans_mult.cc
template <class Matrix>
inline typename Matrix::transpose_type
trans(const Matrix& A) {
  typedef typename Matrix::transpose_type Trans;
  return Trans(A, do_transpose());
}

//: Diagonal Matrix
//
//  This class implements a DiagonalMatrix.
//
//!models: DiagonalMatrix
//!componont: type
//!category: container
template <class TwoDGen, class IndexerGen>
class diagonal_matrix : public matrix_implementation<TwoDGen, IndexerGen> {
  typedef matrix_implementation<TwoDGen, IndexerGen> Base;
  typedef diagonal_matrix self;
public:
  typedef typename Base::pointer pointer;
  typedef typename Base::reference 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 typename Base::const_reference const_reference;


  typedef OneD Diagonal;
  typedef OneDRef DiagonalRef;
  typedef diagonal_tag shape;

  typedef diagonal_matrix<TwoDGen,
                          typename IndexerGen::transpose_type> transpose_type;
  /* these are bogus */
#if 0
  typedef diagonal_matrix<typename TwoDGen::transpose_type,
                          typename IndexerGen::strided_type> strided_type;
#else
  // VC++ workaround
  typedef row_matrix<typename TwoDGen::transpose_type,
                          typename IndexerGen::strided_type> strided_type;
#endif
  inline diagonal_matrix() { }

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

  //: constructor for external data
  inline diagonal_matrix(pointer d, size_type m, size_type n)
    : Base(d, dim_type(m, n)) { }
  inline diagonal_matrix(pointer d, size_type m, size_type n, 
                  int sub, int super)
    : Base(d, dim_type(m, n), band_type(sub, super)) { }

  typedef matrix_market_stream<value_type> mmstream;
  typedef harwell_boeing_stream<value_type> hbstream;

  //: stream constructors
  inline diagonal_matrix(mmstream& m_in) : Base(m_in) { }

  //: stream constructors
  inline diagonal_matrix(hbstream& m_in) : Base(m_in) { }

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

  //: called by strided(A) helper function
  inline diagonal_matrix(const self& x, do_strided s)
    : Base(x) { }

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

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

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

  inline ~diagonal_matrix() { }
};

//: triangle
// example and documentation
template <class Base_, class Uplo>
class triangle_matrix : public Base_ {
  typedef Base_ Base;
  Uplo uplo;
public:
  typedef typename Base::pointer pointer;
  typedef typename Base::reference 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 typename Base::orien orien;
  typedef typename Base::const_reference const_reference;

  typedef triangle_tag shape;

  typedef typename Uplo::transpose_type Uplotrans;
  typedef triangle_matrix<typename Base::transpose_type,
                          Uplotrans> transpose_type;
  /* JGS strided type is bogus */
  typedef triangle_matrix<typename Base::strided_type, Uplo> strided_type;


  inline triangle_matrix() { }

  inline triangle_matrix(size_type m, size_type n)
    : Base(m, n, 
           uplo.bandwidth(m-1,n-1).first,
           uplo.bandwidth(m-1,n-1).second) { }

  //: constructor for external data
  inline triangle_matrix(pointer d, size_type m, size_type n)
    : Base(d, m, n, 
           uplo.bandwidth(m-1, n-1).first,
           uplo.bandwidth(m-1, n-1).second) { }

  //: dynamic uplo constructor
  inline triangle_matrix(size_type m, size_type n, int uplo_)
    : Base(m, n, Uplo::bandwidth(uplo_, m, n)) , uplo(uplo_) { }

  //: constructor for external data with dynamic uplo
  inline triangle_matrix(pointer d, size_type m, size_type n, int uplo_)
    : Base(d, m, n, 
           Uplo::bandwidth(uplo_, m-1, n-1).first,
           Uplo::bandwidth(uplo_, m-1, n-1).second) { }

  //: deprecated
  inline triangle_matrix(pointer d, size_type m, size_type n, 
                         int sub, int super)
    : Base(d, m, n, 
           uplo.bandwidth(sub,super).first,
           uplo.bandwidth(sub,super).second) { }

  //: triangular view constructor
  template <class Matrix>
  inline triangle_matrix(const Matrix& x) 
    : Base(uplo.bandwidth(x.nrows()-1, x.ncols()-1).first,
           uplo.bandwidth(x.nrows()-1, x.ncols()-1).second, x) { }

  //: trans(A)
  inline triangle_matrix(const transpose_type& x, do_transpose t)
    : Base(x, t) { }

  //: rows(A), columns(A) This currently doesn't work
  inline triangle_matrix(const strided_type& x, do_strided s)
    : Base(x, s) { }

  inline triangle_matrix(const triangle_matrix& x, do_strided)
    : Base(x) { }


  typedef matrix_market_stream<value_type> mmstream;
  typedef harwell_boeing_stream<value_type> hbstream;
  inline triangle_matrix(mmstream& m_in) 
    : Base(m_in,
           uplo.bandwidth(m_in.nrows()-1,m_in.ncols()-1).first,
           uplo.bandwidth(m_in.nrows()-1,m_in.ncols()-1).second,
           *this) { }
  inline triangle_matrix(hbstream& m_in) 
    : Base(m_in,
           uplo.bandwidth(m_in.nrows()-1,m_in.ncols()-1).first,
           uplo.bandwidth(m_in.nrows()-1,m_in.ncols()-1).second,
           *this) { }

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

  inline ~triangle_matrix() { }

  inline bool is_upper() const { return uplo.is_upper(); }
  inline bool is_lower() const { return ! uplo.is_upper(); }
  inline bool is_unit() const { return uplo.is_unit(); }

};


//: symmetric

template <class Base_, class Uplo>
class symmetric_matrix : public Base_ {
  typedef Base_ Base;
  Uplo uplo;
public:
  typedef typename Base::pointer pointer;
  typedef typename Base::reference 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 typename Base::orien orien;
  typedef typename Base::const_reference const_reference;

  typedef symmetric_tag shape;

  typedef typename Uplo::transpose_type Uplotrans;
  typedef symmetric_matrix<typename Base::transpose_type,
                           Uplotrans> transpose_type;
  /* JGS strided type is bogus */
  typedef symmetric_matrix<typename Base::strided_type, Uplo> strided_type;

  inline symmetric_matrix() { }

  inline symmetric_matrix(size_type n)
    : Base(n, n, 
           uplo.bandwidth(n-1,n-1).first,
           uplo.bandwidth(n-1,n-1).second) { }
  inline symmetric_matrix(size_type n, int sub)
    : Base(n, n, 
           uplo.bandwidth(sub,sub).first,
           uplo.bandwidth(sub,sub).second) { }

  //: constructor for external data
  inline symmetric_matrix(pointer d, size_type n)
    : Base(d, n, n, 
           uplo.bandwidth(n-1,n-1).first,
           uplo.bandwidth(n-1,n-1).second) { }

  //: dynamic uplo constructor
  inline symmetric_matrix(size_type n, int uplo_, int sub)
    : Base(m, n, Uplo::bandwidth(uplo_, sub, sub)) , uplo(uplo_) { }

  //: constructor for external data with dynamic uplo
  inline symmetric_matrix(pointer d, size_type n, int uplo_, int sub)
    : Base(d, m, n, 
           Uplo::bandwidth(uplo_, sub, sub).first,
           Uplo::bandwidth(uplo_, sub, sub).second) { }

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

  //: deprecated 
  inline symmetric_matrix(pointer d, size_type n, int sub)
    : Base(d, n, n, 
           uplo.bandwidth(sub,sub).first,
           uplo.bandwidth(sub,sub).second) { }

  inline symmetric_matrix(const transpose_type& x, do_transpose t)
    : Base(x, t) { }

  inline symmetric_matrix(const strided_type& x, do_strided s)
    : Base(x, s) { }

  typedef matrix_market_stream<value_type> mmstream;
  typedef harwell_boeing_stream<value_type> hbstream;
  inline symmetric_matrix(mmstream& m_in) 
    : Base(m_in,
           uplo.bandwidth(m_in.nrows()-1,m_in.ncols()-1).first,
           uplo.bandwidth(m_in.nrows()-1,m_in.ncols()-1).second,
           *this) { }
  inline symmetric_matrix(hbstream& m_in) 
    : Base(m_in,
           uplo.bandwidth(m_in.nrows()-1,m_in.ncols()-1).first,
           uplo.bandwidth(m_in.nrows()-1,m_in.ncols()-1).second,
           *this) { }

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

  inline ~symmetric_matrix() { }

  inline bool is_upper() const { return uplo.is_upper(); }
  inline bool is_lower() const { return ! uplo.is_upper(); }

  //: element access
  inline reference
  operator()(size_type row, size_type col) {
    if (uplo.is_upper()) {
      if (col > row)
        return Base::operator()(row, col);
      else
        return Base::operator()(col, row);
    } else {
      if (row > col)
        return Base::operator()(row, col);
      else
        return Base::operator()(col, row);
    }
  }
  //:
  inline const_reference
  operator()(size_type row, size_type col) const {
    if (uplo.is_upper()) {
      if (col > row)
        return Base::operator()(row, col);
      else
        return Base::operator()(col, row);
    } else {
      if (row > col)
        return Base::operator()(row, col);
      else
        return Base::operator()(col, row);
    }
  }
  
  /* bandwidth (is symmetric too) */
  inline int sub() const {
    return MTL_MAX(indexer.super(), indexer.sub());
  }
  inline int super() const { 
    return MTL_MAX(indexer.super(), indexer.sub());
  }

};


/* hermitian (on the TODO list) */


} /* namespace mtl */
#endif