matrix_implementation.h

强大的矩阵模版类

// -*- 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_MATRIX_IMPLEMENTATION_H_
#define _MTL_MATRIX_IMPLEMENTATION_H_

#include "mtl/mtl_iterator.h" /* for advance() */
#include <utility>
#include "mtl/reverse_iter.h"

#include "mtl/oned_part.h"
#include "mtl/scaled2D.h"
#include "mtl/meta_if.h"
#include "mtl/meta_equal.h"
#include "mtl/dimension.h"
#include "mtl/matrix_stream.h"
#include "mtl/linalg_vec.h"
#include "mtl/banded_indexer.h"
#include "mtl/partition.h"
#include "mtl/light_matrix.h"

namespace mtl {

template < class T, class Shape, class Storage, class Orientation>
struct matrix;

template <int MM, int NN> 
class rectangle;

template <int External>
struct dense;


//: The main MTL matrix implementation type.
// This class synthesizes all of the various components
// that go into a matrix and presents the functionality
// to the user with the proper interface.
// The other matrix implementation types derive from 
// this class. This class is not used directly.
//
//!category: container
//!component: type
template <class TwoDGen, class IndexerGen>
class matrix_implementation {
  typedef matrix_implementation<TwoDGen,IndexerGen> self;
public:
  typedef typename IndexerGen::twod_dim_type twoddim;

  enum { M = twoddim::M, N = twoddim::N };
  typedef typename TwoDGen::type TwoD;

  typedef typename TwoD::size_type size_type;
  typedef mtl::dimension<size_type> dyn_dim;

  typedef twod_tag dimension; /* name clash a problem */

  typedef typename IndexerGen::type Indexer;
  
  typedef typename TwoD::value_type OldOneD;
  typedef typename TwoD::reference OldOneDRef;
  typedef typename TwoD::iterator oned_iterator;
  typedef typename TwoD::const_iterator const_oned_iterator;

  typedef typename Indexer::dim_type dim_type;
  typedef typename Indexer::band_type band_type;

  typedef oned_part<OldOneD, OldOneD, typename Indexer::OneDIndexer> OneD;
  typedef oned_part<OldOneD, OldOneDRef, typename Indexer::OneDIndexer> OneDRef;
  typedef typename TwoD::value_type TwoD_value_type;
  typedef typename TwoD_value_type::value_type value_type;
  typedef typename TwoD_value_type::reference reference;
  typedef typename TwoD_value_type::const_reference const_reference;
  typedef typename TwoD_value_type::pointer pointer;
  typedef typename TwoD_value_type::difference_type difference_type;

  typedef typename Indexer::shape shape;
  typedef typename Indexer::orientation orientation;
  typedef typename TwoD::sparsity sparsity;
  typedef typename TwoD::storage_loc storage_loc;

  typedef matrix_implementation< TwoDGen, 
                         typename IndexerGen::transpose_type > transpose_type;

  typedef matrix_implementation< typename TwoDGen::transpose_type, 
                         typename IndexerGen::strided_type > strided_type;

  typedef matrix_implementation< typename TwoDGen::banded_view_type, 
                         typename IndexerGen::strided_type > banded_view_type;

  typedef typename TwoD::strideability strideability;

  typedef matrix_implementation< gen_scaled2D<TwoD, value_type>,
                                 IndexerGen > scaled_type;

  typedef OneD NewOneD;
  
  template <int isConst>
  class _iterator {
    typedef _iterator self;
    typedef typename IF<isConst, const_oned_iterator,oned_iterator>::RET Iterator;
  public:
#if !defined( _MSVCPP_ )
	  typedef typename std::iterator_traits<Iterator>::difference_type
                               difference_type;
#else
	  typedef typename std::iterator_traits<Iterator>::distance_type difference_type;
	  typedef difference_type distance_type;
#endif

    typedef typename IF<isConst, const NewOneD, NewOneD>::RET value_type;
    typedef typename IF<isConst, const NewOneD, NewOneD>::RET reference;
    typedef typename IF<isConst, const NewOneD*, NewOneD*>::RET pointer;

    typedef typename std::iterator_traits<Iterator>::iterator_category 
                               iterator_category;

    typedef difference_type Distance;
    typedef Iterator iterator_type;

    inline difference_type index() const { return iter.index(); }
  
    inline _iterator() { }

    inline _iterator(Iterator x, Indexer ind)
      : iter(x), indexer(ind) { }
    
    inline _iterator(const self& x) 
      : iter(x.iter), indexer(x.indexer) { }

    inline self& operator=(const self& x) { 
      iter = x.iter; indexer = x.indexer; return *this; 
    }

    inline operator Iterator() { return iter; }

    inline Iterator base() const { return iter; }

    inline reference operator*() const {
      typename Indexer::OneDIndexer oned_indexer = indexer.deref(iter);
#if 0
	  typename Iterator::value_type v = *iter; /* VC++ */
#endif
      return reference(*iter, oned_indexer);
	}

    inline self& operator++() { ++iter; return *this; }

    inline self operator++(int) { 
      self tmp = (*this);
      ++(*this); 
      return tmp; 
    }

    inline self& operator--() { --iter; return *this; }

    inline self operator--(int) { 
      self tmp = (*this);
      --(*this); 
      return tmp; 
    }

    inline self operator+(Distance n) const {
      self tmp = (*this);
      tmp += n;
      return tmp;
    }

    inline self& operator+=(Distance n) {
      iter += n; return (*this); 
    }

    inline self operator-(Distance n) const {
      self tmp = (*this); 
      tmp -= n;
      return tmp;
    }

    inline self& operator-=(Distance n) { 
      iter -= n; return (*this); 
    }

    inline value_type operator[](Distance n) const {
      self tmp = (*this);
      return *(tmp += n);
    }

    inline Distance operator-(const self& y) {
      return iter - y.iter;
    }

    inline bool operator==(const self& y) const {
      return iter == y.iter; 
    }

    inline bool operator!=(const self& y) const {
      return iter != y.iter; 
    }

    inline bool operator<(const self& y) const {
      return iter < y.iter; 
    }  

  protected:
    Iterator iter;
    Indexer indexer;
  };

  typedef _iterator<0> iterator;
  typedef _iterator<1> const_iterator;

  typedef reverse_iter<iterator> reverse_iterator;
  typedef reverse_iter<const_iterator> const_reverse_iterator;

  typedef typename Indexer::orienter orien;

  /* constructors */
  inline matrix_implementation() { }

  inline matrix_implementation(dim_type dim, size_type nnz_max)
    : twod(Indexer::twod_dim(orien::map(dim)), nnz_max), 
      indexer(orien::map(dim)) { }

  inline matrix_implementation(dim_type dim)
    : twod(Indexer::twod_dim(orien::map(dim))),
      indexer(orien::map(dim)) { }

  inline matrix_implementation(dim_type dim, band_type bw) 
    : twod(Indexer::twod_dim(orien::map(dim), orien::map(bw)),
           Indexer::twod_band(orien::map(dim), orien::map(bw))), 
      indexer(orien::map(dim), orien::map(bw)) { }

  inline matrix_implementation(const TwoD& x, Indexer ind) 
    : twod(x), indexer(ind) { }

  // copy constructor
  inline matrix_implementation(const matrix_implementation& x)
    : twod(x.twod), indexer(x.indexer) { }

  inline matrix_implementation(const matrix_implementation& x,
                               do_strided s) 
    : twod(x.twod), indexer(x.indexer) { }

  inline self& operator=(const self& x) {
    twod = x.twod; indexer = x.indexer;
    return *this;
  }


  /* char added to fix compiler error with KCC with multiply
     matching constructors */
  inline matrix_implementation(const transpose_type& x,
                               do_transpose, do_transpose) 
    : twod(x.get_twod()), indexer(x.get_indexer(), not_strideable()) { }

  inline matrix_implementation(const strided_type& x, do_strided, do_strided)
    : twod(x.get_twod(), do_transpose(), do_transpose()), 
      indexer(x.get_indexer(), strideable()) { }

  template <class MatrixT, class ScalarT>
  inline matrix_implementation(const MatrixT& x, const ScalarT& y, do_scaled) 
    : twod(x.get_twod(), y), indexer(x.get_indexer()) { }

  /* construct from external (pre-existing) memory,
   *    only for use with external2D
   */

  inline matrix_implementation(pointer data, dim_type dim, char)
    : twod(data, orien::map(dim)), indexer(orien::map(dim)) { }

  inline matrix_implementation(pointer data, dim_type dim, size_type ld)
    : twod(data, orien::map(dim), ld), indexer(orien::map(dim)) { }

  //: With non-zero upper-left corner starts
  inline matrix_implementation(pointer data, dim_type dim, size_type ld,
                               dyn_dim starts, char)
    : twod(data, orien::map(dim), ld, orien::map(starts), char()),
      indexer(orien::map(dim)) { }

  inline matrix_implementation(pointer data, dim_type dim,
                               band_type bw)
    : twod(data, orien::map(dim), orien::map(bw)),
      indexer(orien::map(dim), orien::map(bw)) { }

  inline matrix_implementation(pointer data, dim_type dim, size_type ld, 
                               band_type bw)
    : twod(data, orien::map(dim), ld, orien::map(bw)),
      indexer(orien::map(dim)) { }

  //: compressed2D external data constructor
  inline matrix_implementation(dim_type dim, size_type nnz, 
                        pointer val, size_type* ptrs, size_type* inds)
    : twod(orien::map(dim), nnz, val, ptrs, inds),
      indexer(orien::map(dim)) { }

  //: banded view constructor 
  template <class MatrixT>
  inline matrix_implementation(const MatrixT& x, band_type bw)
    : twod(x.twod, orien::map(bw), banded_tag()),
      indexer(orien::map(dim_type(x.nrows(), x.ncols())), orien::map(bw)) { }

  //: block view constructor
  template <class Matrix, class ST, int BM, int BN>
  inline matrix_implementation(const Matrix& x, 
                               mtl::dimension<ST,BM,BN> bd, char)
    : twod(x.twod, orien::map(bd)), 
      indexer(orien::map(dim_type(x.nrows()/bd.first(), 
                                  x.ncols()/bd.second()))) { }

  //: Static M, N Constructor
  inline matrix_implementation(pointer data)
    : twod(data, orien::map(dim_type(0, 0))), 
      indexer(orien::map(dim_type(0, 0))) { }

  inline matrix_implementation(pointer data, size_type ld)
    : twod(data, orien::map(dim_type(0, 0)), ld),
      indexer(orien::map(dim_type(0, 0))) { }

  //: matrix stream constructor
  typedef matrix_market_stream<value_type> mmstream;
  typedef harwell_boeing_stream<value_type> hbstream;

  template <class Me>
  inline matrix_implementation(mmstream& m_in, Me& me)
    : twod(m_in, orien()), indexer(orien::map(dim_type(m_in.nrows(),
                                                       m_in.ncols())))
  { 
    mtl::initialize(me, m_in);
  }
  template <class Me>
  inline matrix_implementation(hbstream& m_in, Me& me)
    : twod(m_in, orien()), indexer(orien::map(dim_type(m_in.nrows(),
                                                       m_in.ncols())))
  { 
    mtl::initialize(me, m_in);
  }
  template <class Me>
  inline matrix_implementation(mmstream& m_in, band_type bw, Me& me)
    : twod(m_in, orien(), 
           Indexer::twod_band(orien::map(dim_type(m_in.nrows(),
                                                  m_in.ncols())), 
                              orien::map(bw))),
      indexer(orien::map(dim_type(m_in.nrows(), m_in.ncols())),
              orien::map(bw))
  { 
    mtl::initialize(me, m_in);
  }
  template <class Me>
  inline matrix_implementation(hbstream& m_in, band_type bw, Me& me)
    : twod(m_in, orien(), 
           Indexer::twod_band(orien::map(dim_type(m_in.nrows(),
                                                  m_in.ncols())), 
                              orien::map(bw))),
      indexer(orien::map(dim_type(m_in.nrows(), m_in.ncols())),
              orien::map(bw))
  { 
    mtl::initialize(me, m_in);
  }

  inline ~matrix_implementation() { }

  /* iterators */
  inline iterator begin() { return iterator(twod.begin(), indexer); }
  inline iterator end() { return iterator(twod.end(), indexer); }

  inline const_iterator begin() const {
    return const_iterator(twod.begin(),indexer); 
  }
  inline const_iterator end() const { 
    return const_iterator(twod.end(),indexer); 
  }

  inline reverse_iterator rbegin() { return reverse_iterator(end()); }
  inline reverse_iterator rend() { return reverse_iterator(begin()); }

  inline const_reverse_iterator rbegin() const {
    return const_reverse_iterator(end()); 
  }
  inline const_reverse_iterator rend() const {
    return const_reverse_iterator(begin()); 
  }

  /* element access */
  inline typename OneD::reference operator()(size_type i, size_type j) {
    dyn_dim p = indexer.at(dyn_dim(i, j));
    return twod(p.first(), p.second());
  }
  inline typename OneD::const_reference operator()(size_type i,
                                                   size_type j) const {
    dyn_dim p = indexer.at(dyn_dim(i, j));
    return twod(p.first(), p.second());
  }

  /* OneD access */
  inline OneDRef operator[](size_type n) {
    return OneDRef(twod[n], indexer.deref(n));
  }
  inline const OneDRef operator[](size_type n) const {
    return OneDRef((OldOneDRef)twod[n], indexer.deref(n));
  }

  /* size */
  inline size_type nrows() const { return indexer.nrows(); }
  inline size_type ncols() const { return indexer.ncols(); }
  inline size_type noneds() const { return twod.major(); }
  inline size_type major() const { return twod.major(); }
  inline size_type minor() const { return twod.minor(); }
  inline size_type nnz() const { return twod.nnz(); }
  inline size_type capacity() const { return twod.capacity(); }

  /* bandwidth */
  inline int sub() const { return indexer.sub(); }
  inline int super() const { return indexer.super(); }

  /* some shape properties */
  inline bool is_upper() const { return false; }
  inline bool is_lower() const { return false; }
  inline bool is_unit() const { return false; }

  inline const TwoD& get_twod() const { return twod; }