// -*- c++ -*-
//
// $COPYRIGHT$
//
//===========================================================================

#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>
#include <mtl/mtl_set.h>
#include <mtl/matrix_helpers.h>
#include <mxtl/fwd.h>
#include <mxtl/traits.h>
#include <mtl/conj.h>

#include <mtl/dense_matrix.h>
#include <mtl/banded_matrix.h>
#include <mtl/triangle_matrix.h>
#include <mtl/triangle_banded_matrix.h>
#include <mtl/triangle_packed_matrix.h>
#include <mtl/symmetric_matrix.h>
#include <mtl/symmetric_banded_matrix.h>
#include <mtl/symmetric_packed_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 IndexerGen::type Indexer;

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

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

  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 TwoD::sparsity sparsity;
  typedef typename TwoD::storage_loc storage_loc;

  typedef typename Indexer::dim_type dim_type;
  typedef typename Indexer::band_type band_type;
  typedef typename Indexer::shape shape;
  typedef typename Indexer::orientation orientation;
  typedef typename orien_trans<orientation>::tag TransOrien;

  typedef oned_part<OldOneD, OldOneD,
                    typename Indexer::OneDIndexer, orientation> OneD;
  typedef const OneD const_OneD;
  typedef oned_part<OldOneD, OldOneDRef,
                    typename Indexer::OneDIndexer, orientation> OneDRef;
  typedef typename TwoD::value_type TwoD_value_type;

#if 0
  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;
#endif

  typedef typename TwoD_value_type::value_type element_type;
  typedef typename TwoD_value_type::reference element_ref;
  typedef typename TwoD_value_type::const_reference const_element_ref;
  typedef typename TwoD_value_type::pointer element_pointer;
  typedef typename TwoD_value_type::const_pointer const_element_pointer;

  typedef typename TwoD_value_type::difference_type difference_type;

  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 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;

    typedef typename TwoD::size_type size_type;

    inline size_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 Mat, class Scalar>
  inline matrix_implementation(const Mat& x, const Scalar& 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(element_pointer data, dim_type dim, char)
    : twod(data, orien::map(dim)), indexer(orien::map(dim)) { }

  inline matrix_implementation(element_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(element_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(element_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(element_pointer data, dim_type dim, size_type ld, 
                               band_type bw)
    : twod(data, orien::map(dim), ld, orien::map(bw)),
      indexer(orien::map(dim), orien::map(bw)) { }

  //: compressed2D external data constructor
  inline matrix_implementation(dim_type dim, size_type nnz, 
			       element_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 Matrix>
  inline matrix_implementation(const Matrix& 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(element_pointer data)
    : twod(data, orien::map(dim_type(0, 0))), 
      indexer(orien::map(dim_type(0, 0))) { }

  inline matrix_implementation(element_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<element_type> mmstream;
  typedef harwell_boeing_stream<element_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(); }

#if 0
  typedef Cons<size_type, Cons<size_type,Nil> > dims_type;
  inline dims_type dims() const { return cons(nrows(), cons(ncols(), Nil())); }
#else
  inline std::pair<size_type,size_type> dims() const { 
    return std::make_pair(nrows(), ncols()); 
  }
#endif
  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; }
  inline const Indexer& get_indexer() const { return indexer; }

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

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

  /* compressed2D external storage interface */
  inline element_type* get_val() { return twod.get_val(); }
  inline const element_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>,
  public expr< row_matrix<TwoDGen,IndexerGen> >
{
  typedef row_matrix self;
  typedef matrix_implementation<TwoDGen, IndexerGen> Base;
public:
  typedef typename Base::size_type size_type;
#if 0
  typedef typename Base::value_type value_type;
  typedef typename Base::pointer pointer;
  typedef typename Base::reference reference;
  typedef typename Base::const_reference const_reference;
#else
  typedef typename Base::element_type element_type;
  typedef typename Base::element_pointer element_pointer;
  typedef typename Base::const_element_pointer const_element_pointer;
  typedef typename Base::element_ref element_ref;
  typedef typename Base::const_element_ref const_element_ref;
#endif
  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 typename Base::orientation orientation;
  typedef typename Base::TransOrien TransOrien;

  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 */
  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

  typedef matrix_tag category;
  enum { CATEGORY = MATRIX,
	 NROWS = Indexer::M, NCOLS = Indexer::N,
	 ORIENTATION = orientation::id, 
	 SHAPE = shape::id,
	 SPARSITY = sparsity::id,
	 DIMENSION = 2 };

  inline row_matrix() { }

  inline row_matrix& operator=(const row_matrix& x) {
    Base::operator=(x);
    return *this;
  }

  template <class Expr>
  inline row_matrix& operator=(const expr<Expr>& e) {
    assign(*this, e);
    return *this;
  }

  inline row_matrix& operator=(const element_type& t) {
    mtl::set(*this, t);
    return *this;
  }

  inline strided_type get_strided(TransOrien) const {
    return strided_type(*this, do_strided());
  }
  inline self get_strided(orientation) const { return *this; }

  inline transpose_type get_transpose() const {
    return transpose_type(*this, do_transpose());
  }

  //: 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 constructor
  inline row_matrix(element_pointer data, size_type m, size_type n)
    : Base(data, dim_type(m, n), char()) { }

  //: external data constructor with leading dimension
  inline row_matrix(element_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(element_pointer data, size_type m, size_type n, 
                    size_type ld, dyn_dim starts)
    : Base(data, dim_type(m, n), ld, starts, char()) { }

  //: external data banded matrix constructor
  inline row_matrix(element_pointer data, size_type m, size_type n,
                    int sub, int super)
    : Base(data, dim_type(m, n), 
           band_type(sub, super)) { }

  //: external data banded matrix constructor with leading dimension
  inline row_matrix(element_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(element_pointer data)
    : Base(data) { }

  //: compressed2D external data constructor
  inline row_matrix(size_type m, size_type n, size_type nnz,
             element_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<element_type> mmstream;
  typedef harwell_boeing_stream<element_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 Mat, class Scalar>
  inline row_matrix(const Mat& x, const Scalar& y, do_scaled s) 
    : Base(x, y, s) { }

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

  inline ~row_matrix() { }

  /* 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((element_type*)twod.data() 
                            + starts.second() * twod.ld() + starts.first(),
                            m, n, twod.ld());
    } else {
      return submatrix_type((element_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>,
  public expr< column_matrix<TwoDGen, IndexerGen> >
{
  typedef column_matrix self;
  typedef matrix_implementation<TwoDGen, IndexerGen> Base;
public:

#if 0
  typedef typename Base::value_type value_type;
  typedef typename Base::pointer pointer;
  typedef typename Base::reference reference;
  typedef typename Base::const_reference const_reference;
#endif
  typedef typename Base::element_type element_type;
  typedef typename Base::element_pointer element_pointer;
  typedef typename Base::const_element_pointer const_element_pointer;
  typedef typename Base::element_ref element_ref;
  typedef typename Base::const_element_ref const_element_ref;

  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::orientation orientation;
  typedef typename Base::TransOrien TransOrien;

  typedef OneD Column;
  typedef OneDRef ColumnRef;
  enum { M = Indexer::M, N = Indexer::N };
  
  typedef matrix_tag category;
  enum { CATEGORY = MATRIX,
	 NROWS = Indexer::M, NCOLS = Indexer::N,
	 ORIENTATION = orientation::id, 
	 SHAPE = shape::id,
	 SPARSITY = sparsity::id,
	 DIMENSION = 2 };
  
  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 */
  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 strided_type get_strided(TransOrien) const {
    return strided_type(*this, do_strided());
  }
  inline self get_strided(orientation) const { return *this; }

  inline transpose_type get_transpose() const {
    return transpose_type(*this, do_transpose());
  }

  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;
  }

  template <class Expr>
  inline column_matrix& operator=(const Expression<Expr>& e) {
    evaluate(*this, e);
    return *this;
  }

  inline column_matrix& operator=(const element_type& t) {
    mtl::set(*this, t);
    return *this;
  }

  //: external data
  inline column_matrix(element_pointer data, size_type m, size_type n)
    : Base(data, dim_type(m, n), char()) { }

  inline column_matrix(element_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(element_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(element_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(element_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(element_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,
                element_pointer val, size_type* ptrs, size_type* inds)
    : Base(dim_type(m, n), nnz, val, ptrs, inds) { }

  //: streams
  typedef matrix_market_stream<element_type> mmstream;
  typedef harwell_boeing_stream<element_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 Mat, class Scalar>
  inline column_matrix(const Mat& x, const Scalar& y, do_scaled s) 
    : Base(x, y, s) { }

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

  inline ~column_matrix() { }

  /* 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((element_type*)twod.data() 
                            + starts.second() * twod.ld() + starts.first(),
                            m, n, twod.ld());
    } else {
      return submatrix_type((element_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

};

#if 0
// moved to matrix_helpers.h
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());
}
#endif

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


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

  typedef diagonal_matrix_old<TwoDGen,
                          typename IndexerGen::transpose_type> transpose_type;
  /* these are bogus */
#if 0
  typedef diagonal_matrix_old<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 strided_type get_strided(TransOrien) const {
    return strided_type(*this, do_strided());
  }
  inline self get_strided(orientation) const { return *this; }

  inline transpose_type get_transpose() const {
    return transpose_type(*this, do_transpose());
  }

  template <class RHS_>
  inline diagonal_matrix_old& operator=(const Expression<RHS_>& rhs) { 
    typedef Expression<RHS_> RHS;
    typedef typename CreateLeaf<RHS>::Leaf_t Leaf_t;
    evaluate(*this, OpAssign(),
             MakeReturn<Leaf_t>::make(CreateLeaf<RHS>::make(rhs)));
    return *this; 
  }

  inline diagonal_matrix_old() { }

  //: user callable constructors
  inline diagonal_matrix_old(size_type m, size_type n) 
    : Base(dim_type(m, n)) { }
  inline diagonal_matrix_old(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_old(pointer d, size_type m, size_type n)
    : Base(d, dim_type(m, n)) { }
  inline diagonal_matrix_old(pointer d, size_type m, size_type n, 
                  int sub, int super)
    : Base(d, dim_type(m, n), band_type(sub, super)) { }

  inline diagonal_matrix_old(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)) { }


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

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

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

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

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

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

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

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

  inline ~diagonal_matrix_old() { }
};

//: triangle
// example and documentation
template <class Base_, class Uplo>
class triangle_matrix_old : public Base_, 
  public expr< triangle_matrix_old<Base_,Uplo> >
{
  typedef triangle_matrix_old self;
  typedef Base_ Base;
public:
  Uplo uplo;
  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 typename Base::orientation orientation;
  typedef typename Base::TransOrien TransOrien;

  typedef triangle_tag shape;

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


  template <class RHS_>
  inline triangle_matrix_old& operator=(const Expression<RHS_>& rhs) { 
    typedef Expression<RHS_> RHS;
    typedef typename CreateLeaf<RHS>::Leaf_t Leaf_t;
    evaluate(*this, OpAssign(),
             MakeReturn<Leaf_t>::make(CreateLeaf<RHS>::make(rhs)));
    return *this; 
  }

  inline strided_type get_strided(TransOrien) const {
    return strided_type(*this, do_strided());
  }
  inline self get_strided(orien) const { return *this; }

  inline transpose_type get_transpose() const {
    return transpose_type(*this, do_transpose());
  }

  inline triangle_matrix_old() { }

  inline triangle_matrix_old(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_old(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_old(size_type m, size_type n, uplo_e uplo_)
    : Base(m, n, Uplo::bandwidth(uplo_, m, n)) , uplo(uplo_) { }

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

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

  inline triangle_matrix_old(pointer d, size_type m, size_type n, size_type lda,
                         uplo_e uplo_, int sub)
    : Base(d, m, n, lda,
           Uplo::bandwidth(uplo_,sub,sub).first,
           Uplo::bandwidth(uplo_,sub,sub).second), uplo(uplo_) { }

  //: triangular view constructor
  template <class Matrix>
  inline triangle_matrix_old(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_old(const transpose_type& x, do_transpose t)
    : Base(x, t), uplo(x.uplo.trans()) { }

  //: rows(A), columns(A)
  inline triangle_matrix_old(const strided_type& x, do_strided s)
    : Base(x, s), uplo(x.uplo.trans()) { }

  typedef matrix_market_stream<value_type> mmstream;
  typedef harwell_boeing_stream<value_type> hbstream;
  inline triangle_matrix_old(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_old(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 Mat, class Scalar>
  inline triangle_matrix_old(const Mat& x, const Scalar& y, do_scaled s) 
    : Base(x, y, s), uplo(x.uplo) { }

  inline ~triangle_matrix_old() { }

  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 and hermitian
template <class Base_, class Uplo, bool IsHerm >
class symmetric_matrix_old
  : public Base_, public expr< symmetric_matrix_old<Base_, Uplo, IsHerm> >
{
  typedef symmetric_matrix_old self;
  typedef Base_ Base;
public:
  Uplo uplo;
  typedef typename Base::value_type value_type;
  typedef typename IF<IsHerm, hermitian_tag, symmetric_tag>::RET shape;
  typedef typename IF<IsHerm, 
                      conj_proxy<value_type>,
                      typename Base::reference>::RET reference;
  typedef typename IF<IsHerm,
                      value_type,
                      typename Base::const_reference>::RET const_reference;
  typedef typename Base::pointer pointer;
  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::TransOrien TransOrien;
protected:
  bool is_herm(symmetric_tag) const { return false; }
  bool is_herm(hermitian_tag) const { return true; }
public:
  bool is_hermitian() const { return is_herm(shape()); }  

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

  inline strided_type get_strided(TransOrien) const {
    return strided_type(*this, do_strided());
  }
  inline self get_strided(orien) const { return *this; }

  inline transpose_type get_transpose() const {
    return transpose_type(*this, do_transpose());
  }

  template <class RHS_>
  inline symmetric_matrix_old& operator=(const Expression<RHS_>& rhs) { 
    typedef Expression<RHS_> RHS;
    typedef typename CreateLeaf<RHS>::Leaf_t Leaf_t;
    evaluate(*this, OpAssign(),
             MakeReturn<Leaf_t>::make(CreateLeaf<RHS>::make(rhs)));
    return *this; 
  }

  inline symmetric_matrix_old() { }

  inline symmetric_matrix_old(size_type n)
    : Base(n, n, 
           uplo.bandwidth(n-1,n-1).first,
           uplo.bandwidth(n-1,n-1).second) { }
  inline symmetric_matrix_old(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_old(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_old(size_type n, uplo_e uplo_, int sub)
    : Base(n, n, Uplo::bandwidth(uplo_, sub, sub)), uplo(uplo_) { }

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

  //: external data, lda, and dynamic uplo
  inline symmetric_matrix_old(pointer d, size_type n, int lda, uplo_e uplo_, int sub)
    : Base(d, n, n, lda,
           Uplo::bandwidth(uplo_, sub, sub).first,
           Uplo::bandwidth(uplo_, sub, sub).second), uplo(uplo_) { }

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

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

  inline symmetric_matrix_old(const transpose_type& x, do_transpose t)
    : Base(x, t), uplo(x.uplo.trans()) { }

  inline symmetric_matrix_old(const strided_type& x, do_strided s)
    : Base(x, s), uplo(x.uplo.trans()) { }

  typedef matrix_market_stream<value_type> mmstream;
  typedef harwell_boeing_stream<value_type> hbstream;
  inline symmetric_matrix_old(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_old(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 Mat, class Scalar>
  inline symmetric_matrix_old(const Mat& x, const Scalar& y, do_scaled s) 
    : Base(x, y, s) { }

  inline ~symmetric_matrix_old() { }

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

  //: element access
  inline reference
  elt_access(size_type row, size_type col, symmetric_tag) {
    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
  elt_access(size_type row, size_type col, symmetric_tag) 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);
    }
  }
  inline reference
  elt_access(size_type row, size_type col, hermitian_tag) {
    if (uplo.is_upper()) {
      if (col > row)
        return reference(Base::operator()(row, col), false);
      else
        return reference(Base::operator()(col, row), true);
    } else {
      if (row > col)
        return reference(Base::operator()(row, col), false);
      else
        return reference(Base::operator()(col, row), true);
    }
  }
  //:
  inline const_reference
  elt_access(size_type row, size_type col, hermitian_tag) const {
    if (uplo.is_upper()) {
      if (col > row)
        return Base::operator()(row, col);
      else
        return MTL_CONJ(Base::operator()(col, row));
    } else {
      if (row > col)
        return Base::operator()(row, col);
      else
        return MTL_CONJ(Base::operator()(col, row));
    }
  }
  inline reference
  operator()(size_type row, size_type col) {
    return elt_access(row,col, shape());
  }
  inline const_reference
  operator()(size_type row, size_type col) const {
    return elt_access(row,col, shape());
  }
  /* 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());
  }

};

template <class TwoDGen, class IndexerGen>
struct linalg_category< mtl::row_matrix<TwoDGen,IndexerGen> >
{
  typedef matrix_tag type;
  enum { RET = MATRIX };
};

template <class TwoDGen, class IndexerGen>
struct linalg_category< mtl::column_matrix<TwoDGen,IndexerGen> >
{
  typedef matrix_tag type;
  enum { RET = MATRIX };
};

} /* namespace mtl */


#if 0
template <class TwoDGen, class IndexerGen>
struct expr_traits< mtl::row_matrix<TwoDGen,IndexerGen> >
{
  typedef mtl::row_matrix<TwoDGen,IndexerGen> Mat;
  enum { orien = ROW_MAJOR };
  enum { is_scalar = false, is_expr = false, 
	 dim = 2, 
	 size = dynamic_size, 
	 oned_size = dynamic_size,
	 dynamic_sized = true };
  typedef typename Mat::value_type scalar_type;
  //  typedef typename Mat::sparsity sparsity;
  // 2D part always is dense
  typedef mtl::dense_tag sparsity;
};

template <class TwoDGen, class IndexerGen>
struct expr_traits< mtl::column_matrix<TwoDGen,IndexerGen> >
{
  typedef mtl::column_matrix<TwoDGen,IndexerGen> Mat;
  enum { orien = COL_MAJOR };
  enum { is_scalar = false, is_expr = false, 
	 dim = 2, 
	 size = dynamic_size, 
	 oned_size = dynamic_size,
	 dynamic_sized = true };
  typedef typename Mat::value_type scalar_type;
  //  typedef typename Mat::sparsity sparsity;
  typedef mtl::dense_tag sparsity;
};
#endif

#endif