RE: Re: RE: lots of problems VC7.1

// -*- 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_ONED_PART_
#define _MTL_ONED_PART_

/*
  A class template that allows for flexible indexing of a 1-D part of
  a matrix, depending on what shape and storage the matrix is in.
  The Indexer class is responsible for filling in how indexing is
  suppose to work.
*/

/* #include "mtl/iterator_adaptor.h" no longer used */
#include "mtl/reverse_iter.h"
#include "mtl/scaled1D.h"

namespace mtl {


template <class Vector, class VecRef, class Indexer>
class oned_part {
public:
  typedef typename Vector::iterator base_iterator;
  typedef typename Vector::const_iterator const_base_iterator;
  typedef oned_part self;
  enum { N = Vector::N };

  /***** associated types ******/
  typedef typename Vector::value_type value_type;
  typedef typename Vector::reference reference;
  typedef typename Vector::const_reference const_reference;
  typedef typename Vector::pointer pointer;
  typedef typename Vector::size_type size_type;
  typedef typename Vector::difference_type difference_type;
  typedef typename Vector::subrange_type inner_subrange;
  typedef oned_part<inner_subrange, inner_subrange, Indexer> subrange_type;

  typedef typename Vector::sparsity sparsity;
  typedef typename Vector::IndexArrayRef IndexArrayRef;
  typedef typename Vector::IndexArray IndexArray;
  typedef typename Vector::dimension dimension;
  typedef scaled1D< self > scaled_type;

  /***** iterator classes *****/  

  template <int isConst>
  class __iterator { 
    typedef __iterator self;
    typedef typename IF<isConst,const_base_iterator,base_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 std::iterator_traits<Iterator>::value_type value_type;
    typedef typename std::iterator_traits<Iterator>::iterator_category 
                               iterator_category;
#if !defined (_MSVCPP_)
    typedef typename std::iterator_traits<Iterator>::pointer pointer; 
    typedef typename std::iterator_traits<Iterator>::reference reference;
#else
    typedef typename Iterator::pointer pointer;
    typedef typename Iterator::reference reference;
#endif

    typedef difference_type Distance;
    typedef Iterator iterator_type;

    inline __iterator() { }

    inline __iterator(const 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 { return *iter; }

    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 reference operator[](Distance n) const {
      self tmp = (*this);
      return *(tmp += n);
    }
    
    inline Distance 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; 
    }
    
    inline bool operator<(const self& y) const {
      return iter < y.iter; 
    }  
    
    inline size_type row() const { return indexer.row(iter); }
    inline size_type column() const { return indexer.column(iter); }
    inline size_type index() const { return indexer.minor(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;

  /***** oned_part methods *****/

  inline oned_part(VecRef v, Indexer ind) : vec(v), indexer(ind) { }

  inline oned_part(const oned_part& x) : vec(x.vec), indexer(x.indexer) { }

#if !defined(_MSVCPP_)
  template <class OtherOneD>
  inline oned_part(const OtherOneD& x) : vec(x.vec), indexer(x.indexer) { }
#endif

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

  template <class OtherOneD>
  inline oned_part& operator=(const OtherOneD& x) {
    vec = x.vec; indexer = x.indexer; return *this;
  }
  
  inline iterator begin() { /* JGS possibly don't need indexer::begin() */
    return iterator(indexer.begin(vec.begin()), indexer);
  }
  inline iterator end() {
    return iterator(vec.end(), indexer);
  }
  inline const_iterator begin() const {
    return const_iterator(indexer.begin(vec.begin()), indexer);
  }
  inline const_iterator end() const {
    return const_iterator(vec.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()); 
  }
  
  inline reference operator[](size_type n) {
    return vec[indexer.at(n)]; 
  }
  inline const_reference operator[](size_type n) const { 
    return vec[indexer.at(n)]; 
  }

  //added by Rich
  inline value_type& get_ref(size_type j) {
    return vec.get_ref(j);
  }

  inline subrange_type operator()(size_type i, size_type j) const {
    return subrange_type(vec(indexer.at(i), indexer.at(j)), indexer); 
  }

  inline size_type size() const { return vec.size(); }

  //add by Rich
  inline void reserve(size_type i) {
    vec.reserve(i);
  }
  //add by Rich
  inline void resize(size_type i) {
    vec.resize(i);
  }

  inline void clear() { vec.clear(); }
  inline void push_back(size_type i, const value_type& t) {
    vec.push_back(i, t);
  }

  inline size_type nnz() const { return vec.nnz(); }

  inline IndexArrayRef nz_struct() const { return vec.nz_struct(); }

// trying to get mtl_classical_gram_schmidt.h to work
  inline value_type* data()	{ return &get_ref(0); }

  /*
    inline void adjust_index(size_type i) { vec.adjust_index(i); }
    */

  /* protected:*/
  VecRef vec;
  Indexer indexer;
};

} /* namespace mtl */
#endif

//
// 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_UTILS_H
#define MTL_UTILS_H

#include "mtl/mtl_config.h"
#include "mtl/not_at.h"
#include "mtl/entry.h"
#include "mtl/matrix_traits.h"
#include "mtl/dimension.h"

#include <utility>
#include "mtl/mtl_complex.h"
#include <algorithm>
#include <iostream>

namespace mtl {


using std::complex;

/* Utility Functions */


template <class Vector>
inline void
print_partitioned_vector(Vector x)
{
  for (typename Vector::iterator i = x.begin();
       i != x.end(); ++i) {
    mtl::print_vector(*i);
  }
}

template <class Matrix>
inline void
print_partitioned_matrix(const Matrix& A)
{
  std::cout << "Top mat: " << A.nrows() << "x" << A.ncols() << std::endl;
  typedef typename mtl::matrix_traits<Matrix>::size_type Int;
  Int i,j;
  for (i=0; i < A.nrows(); ++i) {
    for (j=0; j < A.ncols(); ++j)
      mtl::print_all_matrix( A(i,j) );
    std::cout << std::endl;
  }
}

template <class Matrix>
inline void
print_partitioned_by_row(const Matrix& A)
{
  std::cout << "Top mat: " << A.nrows() << "x" << A.ncols() << std::endl;
  typename Matrix::const_iterator A_kk;
  typename Matrix::Row::const_iterator A_i, A_iend;
  A_kk = A.begin();
  while (not_at(A_kk, A.end())) {
    A_i = (*A_kk).begin();
    A_iend = (*A_kk).end();
    while (not_at(A_i, A_iend)) {
      mtl::print_all_matrix( *A_i );
      ++A_i;
    }
    ++A_kk;
  } 
}

template <class Matrix>
inline void
print_partitioned_by_column(const Matrix& A)
{
  std::cout << "Top mat: " << A.nrows() << "x" << A.ncols() << std::endl;
  typename Matrix::const_iterator A_kk;
  typename Matrix::Column::const_iterator A_i, A_iend;
  A_kk = A.begin();
  while (not_at(A_kk, A.end())) {
    A_i = (*A_kk).begin();
    A_iend = (*A_kk).end();
    while (not_at(A_i, A_iend)) {
      mtl::print_all_matrix( *A_i );
      ++A_i;
    }
    ++A_kk;
  } 
}

//: utility function for banded
//  should use this in a couple more places
//!noindex:
template <class size_type>
dimension<size_type>
calc_start_fini(int i, int minor, dimension<size_type> bandwidth)
{
  int start = MTL_MAX(i - int(bandwidth.first()), 0);
  int fini = MTL_MIN(i + int(bandwidth.second()) + 1, minor);
  if (start > fini) start = fini;
  return dimension<size_type>(start,fini);
}


template <class Iterator>
inline void
print_vector(Iterator y, Iterator y_end)
{
  std::cout << "[";
  while (not_at(y, y_end)) {
    std::cout << *y << ",";
    ++y;
  }
  std::cout << "]" << std::endl;
}

///
template <class Vector>
inline void
print_vector(Vector x)
{
  typename Vector::iterator i = x.begin();
  std::cout << "[";
  while (not_at(i, x.end())) {
    std::cout << *i << ",";
    ++i;
  }
  std::cout << "]" << std::endl;
}

template <class Vector>
inline void
print_vector_index(Vector x)
{
  typename Vector::iterator i = x.begin();
  std::cout << "[";
  while (not_at(i, x.end())) {
    std::cout << i.index() << ",";
    ++i;
  }
  std::cout << "]" << std::endl;
}

///
template <class Matrix>
inline void
print_coord(const Matrix& A)
{
  typename Matrix::const_iterator A_kk;
  typename Matrix::OneD::const_iterator A_i;
  A_kk = A.begin();
  while (not_at(A_kk, A.end())) {
    A_i = (*A_kk).begin();
    while (not_at(A_i, (*A_kk).end())) {
      std::cout << "(" << A_kk.row() << "," << A_i.column() << ") = " << *A_i << std::endl;
      ++A_i;
    }
    ++A_kk;
  } 
}

///
template <class Matrix>
inline void
print_all_matrix(const Matrix& A)
{
  typedef typename matrix_traits<Matrix>::size_type Int;
  Int i,j;
  std::cout << A.nrows() << "x" << A.ncols() << std::endl;
  std::cout << "[" << std::endl;
  for (i=0; i < A.nrows(); ++i) {
    std::cout << "[";
    for (j=0; j < A.ncols(); ++j) {
      std::cout << A(i,j);
      if (j < A.ncols() - 1)
	std::cout << ",";
    }
    std::cout << "]";
    if (i < A.nrows() - 1)
      std::cout << "," << std::endl;
    else
      std::cout << std::endl;
  }
  std::cout << "]" << std::endl;
}

///
template <class Matrix>
inline void
print_all_banded(const Matrix& A, int lo, int up)
{ 
  typedef typename matrix_traits<Matrix>::size_type Int;
  Int i, j;
  std::cout << A.nrows() << "x" << A.ncols() << std::endl;
  std::cout << "[" << std::endl;
  for (i = 0; i < A.nrows(); ++i) {
    Int first = MTL_MAX(0, int(i) - lo);
    Int last = MTL_MIN(int(A.ncols()), int(i) + up + 1);
    std::cout << "[";
    for (j = 0; j < A.ncols(); ++j) {
      if (j < first || j >= last)
	std::cout << 0;
      else {
	std::cout << A(i,j);
      }
      if (j < A.ncols() - 1)
	std::cout << ",";
    }
    std::cout << "]";
    if (i < A.nrows() - 1)
      std::cout << "," << std::endl;
    else
      std::cout << std::endl;
  }
  std::cout << "]" << std::endl;
}

///
template <class Matrix>
void print_minor(const Matrix& A)
{
  typename Matrix::const_minor_iterator A_kk;
  typename Matrix::MinorVector::const_iterator A_i;
  A_kk = A.begin_minor();
  std::cout << "[" << std::endl;
  while (not_at(A_kk, A.end_minor())) {
    std::cout << "[";
    A_i = (*A_kk).begin();
    while (not_at(A_i, (*A_kk).end())) {
      std::cout << *A_i;
      ++A_i;
      std::cout << ",";
    }
    std::cout << "]";
    std::cout << "," << std::endl;
    ++A_kk;
  } 
  std::cout << "]" << std::endl;
}

///
template <class Matrix>
void print_row(const Matrix& A)
{
  typename Matrix::const_iterator A_kk;
  typename Matrix::Row::const_iterator A_i, A_iend;
  A_kk = A.begin();
  std::cout << "[" << std::endl;
  while (not_at(A_kk, A.end())) {
    std::cout << "[";
    A_i = (*A_kk).begin();
    A_iend = (*A_kk).end();
    while (not_at(A_i, A_iend)) {
      std::cout << *A_i;
      ++A_i;
      std::cout << ",";
    }
    std::cout << "]";
    std::cout << "," << std::endl;
    ++A_kk;
  } 
  std::cout << "]" << std::endl;
}

///
template <class Matrix>
void print_rev_row(const Matrix& A)
{
  typename Matrix::const_reverse_iterator A_kk;
  typename Matrix::Row::const_reverse_iterator A_i;
  typename Matrix::Row::const_reverse_iterator A_iend;
  A_kk = A.rbegin();
  std::cout << "[" << std::endl;
  while (not_at(A_kk, A.rend())) {
    std::cout << "[";
    A_i = (*A_kk).rbegin();
    A_iend = (*A_kk).rend();
    while (not_at(A_i, A_iend)) {
      std::cout << *A_i;
      ++A_i;
      std::cout << ",";
    }
    std::cout << "]";
    std::cout << "," << std::endl;
    ++A_kk;
  } 
  std::cout << "]" << std::endl;
}

///
template <class Matrix>
void print_column(const Matrix& A)
{
  typename Matrix::const_iterator A_kk;
  typename Matrix::Column::const_iterator A_i;
  A_kk = A.begin();
  std::cout << "[" << std::endl;
  while (not_at(A_kk, A.end())) {
    std::cout << "[";
    A_i = (*A_kk).begin();
    while (not_at(A_i, (*A_kk).end())) {
      std::cout << *A_i;
      ++A_i;
      std::cout << ",";
    }
    std::cout << "]";
    std::cout << "," << std::endl;
    ++A_kk;
  } 
  std::cout << "]" << std::endl;
}

///
template <class Matrix>
void print_rev_column(const Matrix& A)
{
  typename Matrix::const_reverse_iterator A_kk;
  typename Matrix::Column::const_reverse_iterator A_i;
  A_kk = A.rbegin();
  std::cout << "[" << std::endl;
  while (not_at(A_kk, A.rend())) {
    std::cout << "[";
    A_i = (*A_kk).rbegin();
    while (not_at(A_i, (*A_kk).rend())) {
      std::cout << *A_i;
      ++A_i;
      std::cout << ",";
    }
    std::cout << "]";
    std::cout << "," << std::endl;
    ++A_kk;
  } 
  std::cout << "]" << std::endl;
}

///
template <class Matrix>
void print_major(const Matrix& A)
{
  typename Matrix::const_iterator A_kk;
  typename Matrix::OneD::const_iterator A_i;
  A_kk = A.begin();
  std::cout << "[" << std::endl;
  while (not_at(A_kk, A.end())) {
    std::cout << "[";
    A_i = (*A_kk).begin();
    while (not_at(A_i, (*A_kk).end())) {
      std::cout << *A_i;
      ++A_i;
      std::cout << ",";
    }
    std::cout << "]";
    std::cout << "," << std::endl;
    ++A_kk;
  } 
  std::cout << "]" << std::endl;
}

#if 1
template <class MatrixA, class MatrixB>
inline bool
matrix_equal(const MatrixA& A, const MatrixB& B)
{
  typedef typename matrix_traits<MatrixA>::size_type Int;
  if (A.nrows() != Int(B.nrows()) || A.ncols() != Int(B.ncols()))
    return false;

  for (Int i = 0; i < A.nrows(); ++i)
    for (Int j = 0; j < A.ncols(); ++j)
      if (A(i,j) != B(i,j)) {
#if !defined(_MSVCPP_)
		  std::cout << "(" << i << "," << j << ") " 
	     << A(i,j) << " != " << B(i,j) << std::endl;
#endif
		  return false;
      }

  return true;
}

#else

template <class MatrixA, class MatrixB>
inline bool
__matrix_equal(const MatrixA& A, const MatrixB& B, row_tag)
{
  typename MatrixA::const_row_2Diterator A_ii = A.begin_rows();
  while (not_at(A_ii, A.end_rows())) {
    typename MatrixA::RowVector::const_iterator A_j = (*A_ii).begin();
     while (not_at(A_j, (*A_ii).end())) {
       if (abs(*A_j - B(A_ii.index(), A_j.index())) > 0.00001) {
#if !defined(_MSVCPP_)
	 std::cout << "(" << A_ii.index() << "," << A_j.index() << ")="
	      << *A_j << " != " << B(A_ii.index(), A_j.index()) << std::endl;
#endif
	return false;
       }
      ++A_j;
    }
    ++A_ii;
  }
  return true;
}


template <class MatrixA, class MatrixB>
inline bool
__matrix_equal(const MatrixA& A, const MatrixB& B, column_tag)
{
  typename MatrixA::const_column_2Diterator A_jj = A.begin_columns();
  while (A_jj < A.end_columns()) {
    typename MatrixA::ColumnVector::const_iterator A_i = (*A_jj).begin();
     while (A_i != (*A_jj).end()) {
       if (abs(*A_i - B(A_i.index(), A_jj.index())) > 0.00001) {
#if !defined(_MSVCPP_)
		   std::cout << "(" << A_i.index() << "," << A_jj.index() << ")=" 
	      << *A_i << " != " << B(A_i.index(), A_jj.index()) << std::endl;
#endif
		   return false;
       }
      ++A_i;
    }
    ++A_jj;
  }
  return true;
}

template <class MatrixA, class MatrixB>
inline bool
matrix_row_equal(const MatrixA& A, const MatrixB& B) {
  return __matrix_equal(A, B, row_tag());
}

template <class MatrixA, class MatrixB>
inline bool
matrix_column_equal(const MatrixA& A, const MatrixB& B) {
  return __matrix_equal(A, B, column_tag());
}

///
template <class MatrixA, class MatrixB>
inline bool
matrix_equal(const MatrixA& A, const MatrixB& B)
{
  typedef typename MatrixA::orientation orien;
  return __matrix_equal(A, B, orien());
}
#endif

/*********************************************************/
/* Functions to get random numbers to fill the matrices  */
/*********************************************************/
  

inline float
make_rand_element(float)
{
  float r = float(rand());
  return r/float(RAND_MAX*10.0);
}

inline double
make_rand_element(double)
{
  double r = double(rand());
  return r/double(RAND_MAX)*10.0L;
}

inline complex<float>
make_rand_element(complex<float>)
{
  return complex<float>(make_rand_element(float()), 
			make_rand_element(float()));
}

inline complex<double>
make_rand_element(complex<double>)
{
  return complex<double>(make_rand_element(double()),
			 make_rand_element(double()));
}

#if 0
template <class Orien>
inline void
print_matrix(const matrix<Orien>& mA)
{
#if 0
  const Orien& A = mA.two_d();
  Orien::const_row_2Diterator row_iter = A.begin_rows();
  std::cout << "[";
  while (row_iter < A.end_rows()) {
    Orien::RowVector::const_iterator ri = (*row_iter).begin();
    std::cout << "[";
     while (ri < (*row_iter).end()) {
      std::cout << *ri;
      std::cout << ",";
      ++ri;
    }
    std::cout << "]," << std::endl;;
    ++row_iter;
  }
  std::cout << "]" << std::endl;
#endif

#if 1
  // good one
  const Orien& A = mA.two_d();
  typename Orien::const_row_2Diterator row_iter = A.begin_rows();
  while (row_iter < A.end_rows()) {
    typename Orien::RowVector::const_iterator ri = (*row_iter).begin();
     while (ri < (*row_iter).end()) {
       if (*ri != 0)
	 std::cout << "(" << row_iter.index() << "," << ri.index()
	      << ") = " << *ri << std::endl;;
      ++ri;
    }
    ++row_iter;
  }
#endif

#if 0
  int i, j;
  for (i=0; i < mA.nrows(); ++i)
    for (j=0; j < mA.ncols(); ++j)
       if (mA(i,j) != 0)
	 std::cout << "(" << i << "," << j
	      << ") = " << mA(i,j) << std::endl;;
#endif

#if 0
  int i,j;
  std::cout << "[";
  for (i=0; i < A.nrows(); ++i) {
    std::cout << "[";
    for (j=0; j < A.ncols(); ++j) {
      std::cout << A(i,j);
      if (j < A.ncols() - 1)
	std::cout << ",";
    }
    std::cout << "]";
    if (i < A.nrows() - 1)
      std::cout << "," << std::endl;
  }
  std::cout << "]" << std::endl;
#endif
}

template <class Orien>
inline void
dense_init(matrix<Orien>& A)
{
  int i, j;
  for(i = 0; i < A.nrows(); ++i)
    for (j = 0; j < A.ncols(); ++j)
      A(i,j) = i * A.ncols() + j;
}

#endif

template <class Mat>
inline void
insert_zero_matrix(Mat& A)
{
  typedef typename matrix_traits<Mat>::size_type Int;
  // idea, for sparse matrices, just zero the non zeroes
  // which should remove them
  Int i, j;
  for (i = 0; i < A.nrows(); ++i) {
    Int first = MTL_MAX(0, int(i) - int(A.sub()));
    Int last = MTL_MIN(int(A.ncols()), int(i) + int(A.super()) + 1);
    for (j = 0; j < A.ncols(); ++j)
      if (j >= first && j < last)
	A(i,j) = 0.0;
  }

}

template <class Matrix>
inline void
zero_matrix(Matrix& A){
  typedef typename matrix_traits<Matrix>::value_type T;
  typename Matrix::iterator oneD_iter = A.begin();
  while (oneD_iter < A.end()) {
    typename Matrix::OneD::iterator i = (*oneD_iter).begin();
    while (i < (*oneD_iter).end()) {
      *i = T(0);
      ++i;
    }
    ++oneD_iter;
  }
}
  
#if 0
template <class Orien>
inline void
sparse_init(matrix<Orien>& A, int c)
{
  int i;
  for (i=0; i < A.nrows(); ++i)
    A(i, i) = 4;
  for (i=0; i < A.nrows() - 1; ++i) {
    A(i, i + 1) = -1.0;
    A(i + 1, i) = -1.0;
  }
  for (i=0; i < A.nrows() - 2; ++i) {
    A(i, i + 2) = -1.0;
    A(i + 2, i) = -1.0;
  }
  for (i=0; i < A.nrows() - 3; ++i) {
    A(i, i + 3) = -1.0;
    A(i + 3, i) = -1.0;
  }
  for (i=0; i < A.nrows() - c; ++i) {
    A(i, i + c) = -1.0;
    A(i + c, i) = -1.0;
  }
}


//
//  mat_mat_multiply
//

template <class OA, class OB, class OC>
inline void
calc_correct(const matrix<OA>& mA, const matrix<OB>& mB,
	     matrix<OC>& mC2)
{
  int i, j, k;
  for (k = 0; k < mA.ncols(); ++k)
    for (j = 0; j < mB.ncols(); ++j)
      for (i = 0; i < mA.nrows(); ++i) {
	mC2(i,j) += mA(i,k) * mB(k,j);	 
      }
}

template <class OA>
inline int calc_nonzeroes(const matrix<OA>& mA)
{
  int nonzeroes = 0;

  const OA& A = mA.two_d();

  typename OA::const_row_2Diterator row_iter = A.begin_rows();
  while (row_iter < A.end_rows()) {
    typename OA::RowVector::const_iterator ri = (*row_iter).begin();
    while (ri.not_at((*row_iter).end())) {
      if (*ri != 0.0)
	++nonzeroes;
      ++ri;
    } 
    ++row_iter;
  }

  return nonzeroes;
}
#endif


} /* namespace mtl */

#endif

#ifndef MTL_BLOCK2D_H
#define MTL_BLOCK2D_H

#include "mtl/dimension.h"
#include "mtl/meta_if.h"

namespace mtl {

template <class Block, class OffsetGen, int M, int N>
struct gen_block2D;

//: Block View TwoD Storage
//
// For use in blocked algorithms with rectangle dense matrices. The
// blocks all have the same size (vs. variable sizes as in a
// partitioned matrix). The matrix objects for each block are not
// stored, they are generated on the fly as they are requested, and
// they are lightweight object on the stack so no overhead is
// incurred.  <p> The blocking size must divide evenly into the
// original matrix size. One good way to ensure this is to partition
// the original matrix into a main region that divides evenly and into
// the blocks, and 3 others edge regions that do not get blocked.  <p>
// Use the block_view type constructor and the blocked function to
// create matrices of this type.
//
//!tparam: Block - The submatrix block, a dense external matrix.
//!tparam: OffsetGen - The Offset generator.
//!category: containers, adaptors
//!component: type
//!example: blocked_matrix.cc
//!definition: block2D.h
//
template <class Block, class OffsetGen>
class block2D {
public:
  typedef block2D<Block, OffsetGen> self;
  typedef typename Block::value_type T;

  enum { M = 0, N = 0, BM = Block::M, BN = Block::N };

  //: The 1D container type
  typedef typename Block::size_type size_type;
  //: The type for differences between iterators
  typedef typename Block::difference_type difference_type;
public:
  typedef dimension<size_type, BM, BN> block_dim_type;
  typedef typename OffsetGen:: MTL_TEMPLATE bind<size_type>::type Offset;
  typedef typename Offset::dim_type dim_type;
public:
  //: This is a dense 2D container
  typedef dense_tag sparsity;
  //: This has external storage
  typedef external_tag storage_loc;
  //: This is strideable
  typedef strideable strideability;

#if 0
  /* bogus ? */
  template <class SubMatrix>
  struct partitioned {
    typedef block2D<SubMatrix, OffsetGen> type;
    typedef gen_block2D<SubMatrix, OffsetGen> generator;
  };
#endif

  class block_vector {
  public:
    typedef Block reference;
    typedef const Block const_reference;
    typedef Block value_type;
    typedef Block* pointer;
    typedef typename Block::size_type size_type;
    typedef typename Block::difference_type difference_type;

    enum { M = 0, N = 0 };

    typedef block_vector subrange_type;
    typedef dense_tag sparsity;
    typedef block_vector IndexArray; /* bogus */
    typedef block_vector IndexArrayRef; /* bogus */
    
    typedef oned_tag dimension; /* bogus */

    template <int isConst>
    class __iterator {
      typedef __iterator self;
    public:
      typedef typename IF<	isConst, 
							typename block_vector::const_reference, 
							typename block_vector::reference>::RET reference;
      typedef typename block_vector::value_type value_type;
      typedef block_vector::pointer pointer;
      typedef typename block_vector::size_type size_type;
      typedef typename block_vector::difference_type difference_type;
      typedef std::random_access_iterator_tag iterator_category;

      inline __iterator(T* s, size_type p, size_type str, size_type ld_,
		      block_dim_type bd)
	: start(s), pos(p), stride(str), ld(ld_), bdim(bd) { }

      inline __iterator(const self& x)
	: start(x.start), pos(x.pos), stride(x.stride), 
	  ld(x.ld), bdim(x.bdim) { }

      inline self& operator=(const self& x) {
	start = x.start; pos = x.pos; stride = x.stride; 
	ld = x.ld; bdim = x.bdim; 
	return *this;
      }

      inline __iterator( ) : start(0), pos(0), stride(0), ld(0) { }

      inline reference operator*() const {
	typedef typename Block::orien orienter;
	return Block(start + pos * stride, 
		     orienter::map(bdim).first(),
		     orienter::map(bdim).second(), 
		     ld);
      }

      inline self& operator++() { ++pos; return *this; }
      inline self& operator+=(size_type n) { pos += n; return *this; }
      inline self operator++(int) { self t = *this; ++(*this); return t; }
      inline self& operator--() { --pos; return *this; }
      inline self& operator-=(size_type n) { pos -= n; return *this; }
      inline self operator--(int) { self t = *this; --(*this); return t; }
      inline bool operator!=(const self& x) const { return pos != x.pos; }
      inline bool operator==(const self& x) const { return pos == x.pos; }
      inline bool operator<(const self& x) const { return pos < x.pos; }
      inline size_type index() const { return pos; }

      T* start;
      size_type pos;
      size_type stride;
      size_type ld;
      block_dim_type bdim;
    };

    typedef __iterator<0> iterator;
    typedef __iterator<1> const_iterator;
  
    inline block_vector(T* s, size_type l, block_dim_type bd, 
		      size_type str, size_type ld_)
    : start(s), len(l), bdim(bd), stride(str * bd.second()), ld(ld_) { }
  
    inline block_vector(const block_vector& x)
      : start(x.start), len(x.len), bdim(x.bdim), 
	stride(x.stride), ld(x.ld) { }

    inline block_vector& operator=(const block_vector& x) {
      start = x.start; len = x.len;
      bdim = x.bdim; stride = x.stride;
      ld = x.ld;
      return *this;
    }

    inline block_vector()
      : start(0), len(0), stride(0) { }

    inline ~block_vector() { }

    inline reference operator[](size_type n) {
      typedef typename Block::orien orienter;
      return Block(start + n * stride, 
		   orienter::map(bdim).first(),
		   orienter::map(bdim).second(), 
		   ld);
    }

    inline const_reference operator[](size_type n) const {
      typedef typename Block::orien orienter;
      return Block(start + n * stride, 
		   orienter::map(bdim).first(),
		   orienter::map(bdim).second(), 
		   ld);
    }

    inline iterator begin() {
      return iterator(start, 0, stride, ld, bdim);
    }
    inline iterator end() {
      return iterator(start, len, stride, ld, bdim);
    }

    inline const_iterator begin() const {
      return const_iterator(start, 0, stride, ld, bdim);
    }
    inline const_iterator end() const {
      return const_iterator(start, len, stride, ld, bdim);
    }
    
    T* start;
    size_type len;
    block_dim_type bdim;
    size_type stride;
    size_type ld;
  };
  //: The 1D container type
  typedef block_vector value_type;
  //: A reference to the value type
  typedef block_vector reference;
  //: The type for pointers to the value type
  typedef block_vector* pointer;

  //: The iterator type
  template <int isConst>
  class __iterator {
    typedef __iterator self;
  public:
    typedef std::random_access_iterator_tag iterator_category;
    typedef block_vector value_type;
    typedef typename IF<isConst, const block_vector, block_vector>::RET reference;
    typedef value_type* pointer;
    typedef typename Block::size_type size_type;
    typedef typename Block::difference_type difference_type;
    
    inline __iterator(T* s, size_type ld_, size_type p, 
		    Offset os, block_dim_type bd)
      : start(s), ld(ld_), pos(p), offset(os), bdim(bd) { }
    
    inline __iterator() : start(0), ld(0), pos(0) { }

    inline __iterator(const self& x)
      : start(x.start), ld(x.ld), pos(x.pos), 
	offset(x.offset), bdim(x.bdim) { }

    inline self& operator=(const self& x) {
      start = x.start; ld = x.ld; pos = x.pos; 
      offset = x.offset; bdim = x.bdim;
      return *this;
    }

    inline reference operator*() const {
      return block_vector(start + offset.oned_offset(pos * bdim.first()),
			  //offset.oned_length(pos * bdim.first()), 
			  // JGS, problem
			  offset.minor(), 
			  bdim, 
			  offset.stride(),
			  ld);
    }
    inline reference operator[](size_type n) const {
      return block_vector((T*)start + 
			  offset.oned_offset((pos + n) * bdim.first()),
			  // offset.oned_length((pos + n) * bdim.first()), 
			  // JGS, problem
			  offset.minor(), 
			  bdim, 
			  offset.stride(),
			  ld);
    }
    
    inline self& operator++() { ++pos; return *this; }
    inline self& operator+=(size_type n) { pos += n; return *this; }
    inline self operator++(int) { self t = *this; ++(*this); return t; }
    inline self& operator--() { --pos; return *this; }
    inline self& operator-=(size_type n) { pos -= n; return *this; }
    inline self operator--(int) { self t = *this; --(*this); return t; }
    inline bool operator!=(const self& x) const { return pos != x.pos; }
    inline bool operator==(const self& x) const { return pos == x.pos; }
    inline bool operator<(const self& x) const { return pos < x.pos; }
    inline size_type index() const { return pos; }
    
  protected:
    T* start;
    size_type pos;
    size_type ld;
    Offset offset;
    block_dim_type bdim;
  };

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

  typedef dimension<size_type> dyn_dim;
  typedef block_dim_type bdt;

  //: Constructor from underlying 2D container
  template <class TwoD>
  inline block2D(TwoD& x, dyn_dim b)
    : data_((T*)x.data()), 
      ld_(x.ld()), block_dim(b), offset(x.major() / bdt(b).first(),  
					x.minor() / bdt(b).second(), 
					x.ld()) { }
  //: Copy Constructor
  inline block2D(const block2D& x)
    : data_(x.data_), ld_(x.ld_), block_dim(x.block_dim), offset(x.offset) { }

  inline const block2D& operator=(const block2D& x) {
    data_ = x.data_; ld_ = x.ld_;
    block_dim = x.block_dim; offset = x.offset;
    return *this;
  }
  //: Default Constructor
  inline block2D() : data_(0), ld_(0) { }
  //: Destructor
  inline ~block2D() { }

  //: Return an iterator pointing to the first 1D container
  inline iterator begin() {
    return iterator(data_, ld_, 0, offset, block_dim);
  }
  //: Return an iterator pointing past the end of the 2D container
  inline iterator end() {
    return iterator(data_, ld_, offset.major(), offset, block_dim);
  }

  //: Return a const iterator pointing to the first 1D container
  inline const_iterator begin() const {
    return const_iterator(data_, ld_, 0, offset, block_dim);
  }
  //: Return a const iterator pointing past the end of the 2D container
  inline const_iterator end() const {
    return const_iterator(data_, ld_, offset.major(), offset, block_dim);
  }

  //: Return a reference to the ith 1D container
  inline block_vector operator[](size_type i) {
    return block_vector(s, offset.oned_length(i * block_dim.first()),
			block_dim, offset.stride(), ld_);
  }

  //: Return a reference to the (i,j) element, where (i,j) is in the 2D coordinate system
  inline Block operator()(size_type i, size_type j) {
    typedef typename Block::orien orienter;
    return Block(data_ + offset.elt(i * block_dim.first(),
				    j * block_dim.second()), 
		 orienter::map(block_dim).first(), 
		 orienter::map(block_dim).second(), 
		 ld_);
  }

  //: Return a const reference to the (i,j) element, where (i,j) is in the 2D coordinate system  
  inline const Block operator()(size_type i, size_type j) const {
    typedef typename Block::orien orienter;
    return Block(data_ + offset.elt(i * block_dim.first(),
				    j * block_dim.second()), 
		 orienter::map(block_dim).first(), 
		 orienter::map(block_dim).second(), 
		 ld_);
  }
  //: The leading dimension
  inline size_type ld() const { return ld_; }



protected:
  T* data_;
  size_type ld_;
  block_dim_type block_dim;
  Offset offset;
};

template <class T, class OffsetGen, int M, int N>
struct gen_external2D;

//: blah
//!noindex:
template <class Block, class OffsetGen, int M, int N>
struct gen_block2D {

  typedef gen_block2D<Block,
           typename OffsetGen::transpose_type, M, N> transpose_type;
  typedef gen_external2D<Block, OffsetGen, M, N> submatrix_type;
  typedef gen_block2D<Block, typename OffsetGen::banded_view_type, M, N>
           banded_view_type;

  typedef block2D<Block, OffsetGen> type;

};


} /* namespace mtl */

#endif /* MTL_BLOCK2D_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_DENSE2D_H
#define MTL_DENSE2D_H

#include "mtl/mtl_iterator.h"
#include <utility>
#include <assert.h>
#include <vector>

#include "mtl/mtl_config.h"
#include "mtl/linalg_vec.h"
#include "mtl/strided1D.h"
#include "mtl/initialize.h"
#include "mtl/reverse_iter.h"
#include "mtl/matrix_traits.h"
#include "mtl/dimension.h"

#ifndef MTL_DISABLE_BLOCKING
#include "mtl/block2D.h"
#endif

namespace mtl {

template <class size_t, int MM, int NN>
class strided_offset;


struct strided_tag { enum { id = 1 }; };
struct not_strided_tag { enum { id = 0 }; };

template <class size_t, int MM, int NN>
class band_view_offset;

template <class size_t, int MM, int NN>
class strided_band_view_offset;

template <int M, int N> struct gen_rect_offset;
template <int M, int N> struct gen_strided_offset;
template <int M, int N> struct gen_banded_offset;
template <int M, int N> struct gen_banded_view_offset;
template <int M, int N> struct gen_strided_band_view_offset;
template <int M, int N> struct gen_packed_offset;


//: Rectangular Offset Class
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class rect_offset {
public:
#if !defined(_MSVCPP_)
  template <class Vec>
  struct bind_oned { 
    typedef Vec type; 
  };
#endif
  typedef not_strided_tag is_strided;
  typedef size_t size_type;
  enum { M = MM, N = NN, IS_STRIDED = 0 };
  typedef dimension<size_type, MM, NN> dim_type;
  typedef dimension<int> band_type;
  typedef strided_offset<size_type, MM, NN> transpose_type;
  typedef strideable strideability;
  // VC++ doesn't like this 
  //friend class transpose_type;

  //what is that for? -- llee
  //inline rect_offset() : dim(4444,4444), ld(4444) { }
  inline rect_offset() : dim(0,0), ld(0) { }
  inline rect_offset(const rect_offset& x) : dim(x.dim), ld(x.ld) { }
  inline rect_offset(size_type m, size_type n, size_type ld_)
    : dim(m, n), ld(ld_) { }
  inline rect_offset(size_type m, size_type n, size_type ld_, band_type)
    : dim(m, n), ld(ld_) { }
  rect_offset(const transpose_type& x); /* see below strided_offset for def */
  inline rect_offset& operator=(const rect_offset& x) {
    dim = x.dim; ld = x.ld; return *this;
  }
  inline size_type elt(size_type i, size_type j) const { return i * ld + j; }
  inline size_type oned_offset(size_type i) const { return i * ld; }
  inline size_type oned_length(size_type) const { return dim.second(); }
  inline size_type twod_length() const { return dim.first(); }
  inline size_type stride() const { return 1; }
  inline static size_type size(size_type m, size_type n,
                               size_type , size_type) { return m * n; }
  inline size_type major() const { return dim.first(); }
  inline size_type minor() const { return dim.second(); }
  /* private: */
  dim_type dim;
  size_type ld;
};


//: blah
//!noindex:
template <int M, int N>
struct gen_rect_offset {
#if defined( _MSVCPP_ )
   typedef rect_offset<unsigned int, M, N> type;
#else
  template <class size_type>
  struct bind {
    typedef rect_offset<size_type, M, N> type;
  };
#endif
  typedef gen_strided_offset<M,N> transpose_type;
  typedef gen_banded_view_offset<M,N> banded_view_type;

};

//: Strided Rectangular Offset Class
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class strided_offset {
public:
#if !defined(_MSVCPP_)
  template <class Vec>
  struct bind_oned { 
    typedef strided1D<Vec> type;
  };
#endif
  /*  typedef strided_band_view_offset<size_t,MM,NN> banded_view_type;
   */
  typedef strided_tag is_strided;
  enum { M = MM, N = NN, IS_STRIDED = 1 };
  typedef size_t size_type;
  typedef dimension<size_type, MM, NN> dim_type;
  typedef dimension<int> band_type;
  typedef rect_offset<size_type,MM,NN> transpose_type;
  typedef strideable strideability;
// VC++ doesn't like this
  //friend class transpose_type;
  inline strided_offset() : dim(0,0), ld(0) { }
  inline strided_offset(size_type m, size_type n, size_type ld_)
    : dim(m, n), ld(ld_) { }
  inline strided_offset(const transpose_type& x) : dim(x.dim), ld(x.ld) { }
  inline strided_offset& operator=(const strided_offset& x) {
    dim = x.dim; ld = x.ld; return *this;
  }
  inline size_type elt(size_type i, size_type j) const { return j * ld + i; }
  inline size_type oned_offset(size_type i) const { return i; }
  inline size_type oned_length(size_type) const { return dim.first() * ld; }
  inline size_type twod_length() const { return dim.second(); }
  inline size_type stride() const { return ld; }
  inline static size_type size(size_type m, size_type n,
                               size_type , size_type) { return m * n; }
  inline size_type major() const { return dim.first(); }
  inline size_type minor() const { return dim.second(); }
  /* private: */
  dim_type dim;
  size_type ld;
};

//: blah
//!noindex:
template <int M, int N>
struct gen_strided_offset {
#if defined( _MSVCPP_ )
  typedef strided_offset<unsigned int, M, N> type;
#else
  template <class size_type>
  struct bind {
    typedef strided_offset<size_type, M, N> type;
  };
#endif
  typedef gen_rect_offset<M,N> transpose_type;
  typedef gen_strided_band_view_offset<M,N> banded_view_type;
};

template <class size_t, int MM, int NN>
inline rect_offset<size_t,MM,NN>::rect_offset( typename const rect_offset<size_t,MM,NN>::transpose_type& x)
  : dim(x.dim), ld(x.ld) { }


//: Banded View Offset Class
// This creates a banded view into a full matrix.
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class banded_view_offset {
public:
#if !defined(_MSVCPP_)
  template <class Vec>
  struct bind_oned { 
    typedef Vec type; 
  };
#endif

  typedef not_strided_tag is_strided;
  enum { M = MM, N = NN, IS_STRIDED = 0 };
  typedef size_t size_type;
  typedef dimension<size_type, MM, NN> dim_type;
  typedef dimension<int, MM, NN> band_type;
  typedef strided_band_view_offset<size_type, MM, NN> transpose_type;

  typedef not_strideable strideability;
// VC++ doesn't like this
  //friend class transpose_type;
  inline banded_view_offset()
    : dim(0,0), ld(0), bw(std::make_pair(0,0)) { }
  inline banded_view_offset(size_type m, size_type n, size_type leading_dim,
                            band_type band)  
    : dim(m, n), ld(leading_dim), bw(band) { }
  inline banded_view_offset(size_type m, size_type n, size_type leading_dim)  
    : dim(m, n), ld(leading_dim), bw(band_type(0,0)) { }

  template <class Offset>
  inline banded_view_offset(Offset os, band_type band)  
    : dim(os.dim), ld(os.ld), bw(band) { }

  inline banded_view_offset& operator=(const banded_view_offset& x) {
    dim = x.dim; ld = x.ld; bw = x.bw; return *this;
  }

  inline size_type elt(size_type i, size_type j) const { 
    size_type start = MTL_MAX(int(i) - bw.first(), 0);
    return i * ld + j + start; 
  }
  inline size_type oned_offset(size_type i) const { 
    size_type start = MTL_MAX(int(i) - bw.first(), 0);
    return i * ld + start; 
  }
  inline size_type oned_length(size_type i) const { 
    return MTL_MAX(0, MTL_MIN(int(dim.second()), int(i) + bw.second() + 1)
               - MTL_MAX(0, int(i) - bw.first()));
  }
  inline size_type twod_length() const { return dim.first(); }
  inline int stride() const { return 1; }
  inline static size_type size(size_type m, size_type n, 
                               size_type , size_type) {
    return m * n; 
  }
  inline size_type major() const { return dim.first(); }
  inline size_type minor() const { return dim.second(); }
  
  /* private: */
  dim_type dim;
  size_type ld;
  band_type bw; /* bandwidth */
};


//: blah
//!noindex:
template <int M, int N>
struct gen_banded_view_offset {
#if defined( _MSVCPP_ )
  typedef banded_view_offset<unsigned int, M, N> type;
#else
  template <class size_type>
  struct bind {
    typedef banded_view_offset<size_type, M, N> type;
  };
#endif
  typedef gen_strided_band_view_offset<M,N> transpose_type;
  typedef gen_banded_view_offset<M,N> banded_view_type; /* bogus */
};


//: Strided Band View Offset Class
//
// This creates a strided band view into a full matrix.
// This class is to banded_view as strided_offset is to rect_offset.
//
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class strided_band_view_offset {
public:
#if !defined(_MSVCPP_)
  template <class Vec>
  struct bind_oned { 
    typedef strided1D<Vec> type; 
  };
#endif
  typedef strided_tag is_strided;
  enum { M = MM, N = NN, IS_STRIDED = 1 };
  typedef size_t size_type;
  typedef dimension<size_type, MM, NN> dim_type;
  typedef dimension<int, MM, NN> band_type;
  typedef banded_view_offset<size_type, MM, NN> transpose_type;

  typedef not_strideable strideability;
// VC++ doesn't like this
  //friend class transpose_type;
  inline strided_band_view_offset()
    : dim(0,0), ld(0), bw(std::make_pair(0,0)) { }
  inline strided_band_view_offset(size_type m, size_type n, 
				  size_type leading_dim,
				  band_type band)  
    : dim(m, n), ld(leading_dim), bw(band) { }

  template <class Offset>
  inline strided_band_view_offset(Offset os, band_type band)  
    : dim(os.dim), ld(os.ld), bw(band) { }

  inline strided_band_view_offset& 
  operator=(const strided_band_view_offset& x) {
    dim = x.dim; ld = x.ld; bw = x.bw; return *this;
  }

  inline size_type elt(size_type i, size_type j) const { 
    size_type start = MTL_MAX(int(i) - bw.first(), 0);
    return (j + start) * ld + i; 
  }
  inline size_type oned_offset(size_type i) const { 
    size_type start = MTL_MAX(int(i) - bw.first(), 0);
    return start * ld + i; 
  }
  inline size_type oned_length(size_type i) const { 
    /* use dim.first() here */
    size_type len =  MTL_MAX(0, MTL_MIN(int(dim.first()), int(i) + bw.second() + 1)
                          - MTL_MAX(0, int(i) - bw.first()));
    return len * ld;
  }
  inline size_type twod_length() const { return dim.second(); }
  inline int stride() const { return ld; }
  inline static size_type size(size_type m, size_type n, 
                               size_type , size_type) {
    return m * n; 
  }
  inline size_type major() const { return dim.first(); }
  inline size_type minor() const { return dim.second(); }
  
  /* private: */
  dim_type dim;
  size_type ld;
  band_type bw; /* bandwidth */
};


//: blah
//!noindex:
template <int M, int N>
struct gen_strided_band_view_offset {
#if defined( _MSVCPP_ )
  typedef strided_band_view_offset<unsigned int, M, N> type;
#else
  template <class size_type>
  struct bind {
    typedef strided_band_view_offset<size_type, M, N> type;
  };
#endif
  typedef gen_banded_view_offset<M,N> transpose_type;
  typedef gen_strided_band_view_offset<M,N> banded_view_type; /* bogus */
};


template <class size_t, int MM, int NN>
class packed_offset;


//: Banded Offset Class
// This cooresponds to lapack/blas banded storage format.
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class banded_offset {
public:
#if !defined(_MSVCPP_)
  template <class Vec>
  struct bind_oned { 
    typedef Vec type; 
  };
#endif
  typedef not_strided_tag is_strided;
  enum { M = MM, N = NN, IS_STRIDED = 0 };
  typedef size_t size_type;
  typedef dimension<size_type, MM, NN> dim_type;
  typedef dimension<int> band_type;
  typedef packed_offset<size_type,MM,NN> transpose_type; /* bogus */
  typedef not_strideable strideability;
  inline banded_offset()
    : dim(0,0), bw(band_type(0,0)), ndiag(0) { }

  inline banded_offset(size_type m, size_type n, size_type /* lead */, 
                       band_type band) 
    : dim(m,n), bw(band), ndiag(band.first() + band.second() + 1) { }

 inline banded_offset(size_type m, size_type n, size_type /* lead */) 
    : dim(m,n), bw(band_type(0,0)), ndiag(0) { }

  inline banded_offset& operator=(const banded_offset& x) {
    dim = x.dim; ndiag = x.ndiag; bw = x.bw; return *this;
  }
  inline size_type elt(size_type i, size_type j) const {
    return this->oned_offset(i) + j; 
  }
  inline size_type oned_offset(size_type i) const {
    return i * ndiag + MTL_MAX(0, bw.first() - int(i));
  }
  inline size_type oned_length(size_type i) const {
    return MTL_MAX(0, MTL_MIN(int(dim.second()), int(i) + bw.second() + 1)
               - MTL_MAX(0, int(i) - bw.first()));
  }
  inline size_type twod_length() const { return dim.first(); }

  inline int stride() const { return 1; }

  inline static size_type size(size_type m, size_type n, 
                               size_type low, size_type up) { 
    /* M' = number of diagonals = low + up + 1
       N' = min (m, n + low) */
    return (low + up + 1) * MTL_MIN(m, n + low);
  }
  inline size_type major() const { return dim.first(); }
  inline size_type minor() const { return dim.second(); }
private:
  dim_type dim;
  band_type bw; /* bandwidth */
  size_type ndiag;
};


//: blah
//!noindex:
template <int M, int N>
struct gen_banded_offset {
#if defined( _MSVCPP_ )
  typedef banded_offset<unsigned int, M, N> type;
#else
  template <class size_type>
  struct bind {
    typedef banded_offset<size_type, M, N> type;
  };
#endif
  typedef gen_packed_offset<M,N> transpose_type; // bogus
  typedef gen_banded_view_offset<M,N> banded_view_type; /* bogus */
};


//: Packed Offset Class
// This cooresponds to lapack/blas packed storage format 
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class packed_offset {
public:
#if !defined(_MSVCPP_)
  template <class Vec>
  struct bind_oned { 
    typedef Vec type;
  };
#endif
  typedef not_strided_tag is_strided;
  enum { M = MM, N = NN, IS_STRIDED = 0 };
  typedef size_t size_type;
  typedef dimension<size_type, MM, NN> dim_type;
  typedef dimension<int> band_type;
  typedef banded_offset<size_type, MM, NN> transpose_type; /* bogus */
  typedef not_strideable strideability;
  inline packed_offset() 
    : dim(0,0), bw(band_type(0,0)) { }

  inline packed_offset(size_type m, size_type n, size_type /* lead */,
                       band_type bandwidth)
    : dim(m,n), bw(bandwidth) { }

  inline packed_offset& operator=(const packed_offset& x) {
    dim = x.dim; bw = x.bw; return *this;
  }
  inline int elt(size_type i, size_type j) const { 
    return this->oned_offset(i) + j; 
  }

  inline int calc_low(int i, int low) const {
    int l = MTL_MIN(low, int(i));
    int lower_area = low * i;
    lower_area -= ( - l*l + 2*low*l + l) / 2;
    return lower_area;
  }
  inline int calc_up(int i, int up) const {
    int upper_area = up * i;
    int n = i + up - dim.second();
    if (n > 0) {
      int n1 = MTL_MAX(n - up, 0);
      int n2 = n - n1;
      upper_area -= n1 * up;
      upper_area -= ((n2 + 1) * n2) / 2;
    }
    return upper_area;
  }

  inline int oned_offset(size_type i) const { /* the ith major container */
    int low = bw.first();
    int up = bw.second();
    int upper_area, lower_area;

    if (up < -1)
      upper_area = - calc_low(i, - (up + 1));
    else if (up > 0)
      upper_area = calc_up(i, up);
    else
      upper_area = 0;

    if (low < -1)
      lower_area = - calc_up(i, - (low + 1));
    else if (low > 0)
      lower_area = calc_low(i, low);
    else
      lower_area = 0;

    size_type diagonal_len;
    if (up < 0 || low < 0)
      diagonal_len = 0;
    else
      diagonal_len = MTL_MIN(MTL_MIN(i, dim.first()), dim.second());

    size_type ret =  upper_area + lower_area + diagonal_len;
    return ret;
  }

  inline int stride() const { return 1; }

  inline size_type oned_length(size_type i) const {
    return MTL_MAX(0, MTL_MIN(int(dim.second()), int(i) + bw.second() + 1)
               - MTL_MAX(0, int(i) - bw.first()));
  }
  inline size_type twod_length() const { return dim.first(); }

  inline static size_type size(int m, int n, int low, int up) { 
    packed_offset offset(m, n, n, band_type(low, up));
    return offset.oned_offset(m);
  }
  inline size_type major() const { return dim.first(); }
  inline size_type minor() const { return dim.second(); }

private:
  dim_type dim;
  band_type bw; /* bandwidth */
};

//: blah
//!noindex:
template <int M, int N>
struct gen_packed_offset {
#if defined( _MSVCPP_ )
  typedef packed_offset<unsigned int, M, N> type;
#else
  template <class size_type>
  struct bind {
    typedef packed_offset<size_type, M, N> type;
  };
#endif
  typedef gen_banded_offset<M,N> transpose_type; /* bogus */
  typedef gen_banded_view_offset<M,N> banded_view_type; /* bogus */
};


/* egcs doesn't "see" the friend functions
 *   when the dense2D_iterator class is defined inside of dense2D
 */

//: blah
//!noindex:
template <int isConst, class T, class Offset, class InnerOneD, class OneD>
class dense2D_iterator {
public:
  typedef typename Offset::size_type size_type;
  typedef std::pair<size_type,size_type> pair_type;

  typedef typename IF<isConst, const T*,T*>::RET Iterator;
    
  typedef dense2D_iterator self;

  typedef int distance_type;
  typedef int difference_type;
    
  typedef std::random_access_iterator_tag iterator_category;

  typedef OneD*           pointer;
  typedef OneD            value_type;
  typedef OneD            reference;
  typedef difference_type Distance;
  typedef Iterator        iterator_type;
    
protected:
    
  Iterator start;
  size_type pos;
  size_type ld;    /* leading dimension */
  pair_type starts;
  Offset offset;
public:
    
  inline size_type index() const { return pos + starts.second; }
    
  inline dense2D_iterator () {}
    
  inline dense2D_iterator(const self& x)
    : start(x.start), pos(x.pos), 
      ld(x.ld), starts(x.starts), offset(x.offset) { }

  inline self& operator=(const self& x) {
    start = x.start;
    pos = x.pos;
    ld = x.ld;
    starts = x.starts;
    offset = x.offset;
    return *this;
  }

  inline explicit
  dense2D_iterator(Iterator x, size_type ld_, size_type p, pair_type s,
                   Offset os)
    : start(x), pos(p), ld(ld_), starts(s), offset(os) { }
    
  inline Iterator base () const { return start + pos; }
    
  inline reference deref(Distance pos, not_strided_tag) const { 
    return reference((T*)start + offset.oned_offset(pos), 
		offset.oned_length(pos), 
		starts.first);
  }
  inline reference deref(Distance pos, strided_tag) const { 
    InnerOneD vec((T*)start + offset.oned_offset(pos), 
                  offset.oned_length(pos), 
                  starts.first);
    return strided(vec, offset.stride());
  }
  inline reference operator*() const { 
    typedef typename Offset::is_strided Strided;
    return deref(pos, Strided());
  }
  inline reference operator[] (Distance n) const { 
    typedef typename Offset::is_strided Strided;
    return deref(pos + n, Strided());
  }

  /*  won't work, the OneD is temporary 
  pointer   operator-> () const { return & (operator* ()); }
  */
    
  inline self& operator++ () { ++pos; return *this; }
  inline self operator++ (int) { self tmp = *this; ++pos; return tmp; }
  inline self& operator-- () { --pos; return *this; }
  inline self operator-- (int) { self tmp = *this; --pos; return tmp; }
  inline self& operator+=(size_type n) { pos += n; return *this; }
  inline self operator+(size_type n) const {
    return self(start, ld, pos + n, starts); 
  }
  inline self& operator-=(size_type n) { pos -= n; return *this; }
    
    
    
};

template <int isConst, class T, class Offset, class InnerOneD, class OneD>
inline typename dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>::difference_type
operator-(const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& x,
          const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& y)
{
  return x.index() - y.index(); 
}
    
template <int isConst,class T, class Offset, class InnerOneD, class OneD>
inline bool
operator== (const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& x,
            const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& y)
{
  return x.index() == y.index();
}
    
template <int isConst, class T, class Offset, class InnerOneD, class OneD>
inline bool
operator!= (const dense2D_iterator<isConst, T,Offset,InnerOneD,OneD>& x,
            const dense2D_iterator<isConst, T,Offset,InnerOneD,OneD>& y)
{
  return x.index() != y.index();
}
    
template <int isConst,class T, class Offset, class InnerOneD, class OneD>
inline bool
operator< (const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& x,
           const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& y)
{
  return x.index() < y.index(); 
}



/*
  Workaround (g++ 2.91) helper class
 */

template <class Strided>
struct __bracket { };

template <>
struct __bracket<strided_tag> { 
  template <class OneD, class InnerOneD, class elt_type, class size_type>
  inline OneD
  operator()(elt_type* d, size_type len, size_type f, size_type ld,
             const OneD*, const InnerOneD*) {
    InnerOneD vec(d , len, f);
    return OneD(vec, ld);
  }
}; 
 
template <>
struct __bracket<not_strided_tag> { 
  template <class OneD, class InnerOneD, class elt_type, class size_type>
  inline OneD
  operator()(elt_type* d, size_type len, size_type f, size_type,
             const OneD*, const InnerOneD*) {
    return OneD(d, len, f);
  }
};



template<class T, class OffsetGen, int MM, int NN>
class dense2D;

template <class T, class OffsetGen, int MM, int NN>
class external2D;

//: Generic Dense 2-D Container
//!category: containers
//!component: type
//
// The generic_dense2D container implements sevaral of the MTL storage
// types.  They include dense, packed, banded, and banded_view.  The
// common theme here is that the matrix is stored in a contiguous
// piece of memory.  The differences in these storage types has to do
// with where to find the OneD segements in the linear
// memory. Caclulating these offsets is the job of the Offset concept,
// which has a model to handle each of the different storage types:
// rect_offset, strided_offset, banded_offset, packed_offset, and
// banded_view_offset.
//
// There are two derived classes of generic_dense2D that specify the
// memory management, dense2D and external2D. The dense2D version owns
// its memory, while the external2D imports its memory from somewhere
// else through a pointer (which allows for interoperability with
// other codes -- even with Fortran!).  <p>
//
//!definition: dense2D.h
//!tparam: RepType - The Container used to store the elements
//!tparam: RepPtr - The type used to reference to the container
//!tparam: OffsetGen - The generator that creates the Offset class
//!tparam: MM - For static sized matrix, the major dimension
//!tparam: NN - For static sized matrix, the minor dimension
//!models: TwoDStorage

template <class RepType, class RepPtr, class OffsetGen, int MM, int NN>
class generic_dense2D {
public:
  //: Static sizes (0 if dynamic)
  enum { M = MM, N = NN };

  //: The type for dimensions and indices
  typedef typename RepType::size_type size_type;
  //: The type for differences between iterators
  typedef typename RepType::difference_type difference_type;

protected:
  typedef std::pair<size_type,size_type> pair_type;
  typedef RepType reptype;
  typedef RepPtr rep_ptr;
  typedef typename RepType::value_type elt_type;

#if defined(_MSVCPP_)
  //JGS Nasty VC++ workaround
  typedef typename OffsetGen::type Offset;
#else
  typedef typename OffsetGen:: MTL_TEMPLATE bind<size_type>::type Offset;
#endif

  typedef dimension<elt_type> dyn_dim;
public:
  //: A pair type for dimensions
  typedef typename Offset::dim_type dim_type;

  //: A pair type for bandwidth
  typedef typename Offset::band_type band_type;

  /* Type Definitions */
  
  //: This is a dense matrix
  typedef dense_tag sparsity;

  typedef typename Offset::is_strided is_strided;

protected:
  typedef external_vec<elt_type, N> InnerOneD;

#if defined(_MSVCPP_)
  enum { offset_strided = Offset::IS_STRIDED };
  typedef typename IF<offset_strided, strided1D<InnerOneD>, InnerOneD>::RET OneD;
#else
  typedef typename Offset:: MTL_TEMPLATE bind_oned<InnerOneD>::type OneD;
#endif

  typedef OneD OneDRef;
  typedef OneD ConstOneDRef;
public:

  //: The 1D container type
  typedef OneD value_type;
  //: The type for a reference to value_type
  typedef value_type reference;
  //: The type for a const reference to value_type
  typedef value_type const_reference;

  //: The iterator type
  typedef dense2D_iterator<0,elt_type, Offset, InnerOneD, OneD> iterator;

  //: The const iterator type
  typedef dense2D_iterator<1,elt_type, Offset, InnerOneD, OneD> const_iterator;

  //: The reverse iterator type
  typedef reverse_iter<iterator> reverse_iterator;

  //: The const reverse iterator type
  typedef reverse_iter<const_iterator> const_reverse_iterator;

  //: The type for the transpose of this container
  typedef generic_dense2D<RepType, RepPtr,
             typename OffsetGen::transpose_type, MM, NN> transpose_type;

  //: The type for a banded view of this container
  typedef generic_dense2D<RepType, RepPtr,
             typename OffsetGen::banded_view_type, MM, NN> banded_view_type;

  //: The type for a sub-section of this 2D container
  typedef external2D<elt_type, OffsetGen, MM, NN> submatrix_type;

#ifndef MTL_DISABLE_BLOCKING
  template <class Block>
  struct blocked_view {
    typedef block2D<Block, OffsetGen> type;
  };
#endif

  //: This is a stridable container, can use rows(A), columns(A)
  typedef typename Offset::strideability strideability;

  /* Constructors */
  
  //: Default Constructor
  inline generic_dense2D()
    : ld_(0), data_(0), starts(std::make_pair(0,0)) { }

  //: Normal Constructor
  inline generic_dense2D(rep_ptr data, size_type m, size_type n, size_type ld)
    : ld_(ld), data_(data),
      starts(std::make_pair(0,0)), offset(m, n, ld) { }

  //: Constructor with non-zero upper-left corner indices
  inline generic_dense2D(rep_ptr data, size_type m, size_type n, 
			 size_type ld, dyn_dim s, char)
    : ld_(ld), data_(data),
      starts(std::make_pair(s.first(),s.second())), offset(m, n, ld) { }

  //: Static M, N constructor
  inline generic_dense2D(rep_ptr data, size_type ld)
    : ld_(ld), data_(data),
      starts(std::make_pair(0,0)), offset(M, N, ld) { }

  //: with bandwidth constructor
  inline generic_dense2D(rep_ptr data, size_type m, size_type n, size_type ld,
                         band_type bw)
    : ld_(ld), data_(data), starts(std::make_pair(0,0)),
      offset(m, n, ld, bw) { }

  //: Static M, N with bandwith?

  //: Copy Constructor
  inline generic_dense2D(const generic_dense2D& x)
    : ld_(x.ld_), data_(x.data_), starts(x.starts),
      offset(x.offset) { }

  //: Assignment Operator
  inline generic_dense2D& operator=(const generic_dense2D& x) {
    ld_ = x.ld_; data_ = x.data_; starts = x.starts; offset = x.offset;
    return *this;
  }

  //: Subclass Constructor
  inline generic_dense2D(rep_ptr d, const generic_dense2D& x)
    : ld_(x.ld_), data_(d), starts(x.starts), offset(x.offset) { }

  //: Transpose Constructor
  inline generic_dense2D(const transpose_type& x, do_transpose, do_transpose)
    : ld_(x.ld_), data_(x.data_), starts(x.starts), offset(x.offset) { }      

  /* JGS, remove stream constructor, not very necessary
     just have them call another constructor 
   */

  //: Matrix Stream Constructor
  template <class MatrixStream, class Orien>
  inline generic_dense2D(rep_ptr data, MatrixStream& s, Orien)
    : ld_(Orien::map(dim_type(s.nrows(),s.ncols())).second()),
      data_(data),
      starts(std::make_pair(0,0)),
      offset(Orien::map(dim_type(s.nrows(),s.ncols())).first(), 
             Orien::map(dim_type(s.nrows(),s.ncols())).second(),
             Orien::map(dim_type(s.nrows(),s.ncols())).second()) { }
  //: Banded Matrix Stream Constructor
  template <class MatrixStream, class Orien>
  inline generic_dense2D(rep_ptr data, MatrixStream& s, 
                         Orien, band_type bw)
    : ld_(Orien::map(dim_type(s.nrows(),s.ncols())).second()),
      data_(data), 
      starts(std::make_pair(0,0)),
      offset(Orien::map(dim_type(s.nrows(),s.ncols())).first(), 
             Orien::map(dim_type(s.nrows(),s.ncols())).second(),
             Orien::map(dim_type(s.nrows(),s.ncols())).second(),
             bw) { }

  //: Banded View Constructor
  template <class TwoD>
  inline generic_dense2D(rep_ptr data, const TwoD& x, band_type bw, banded_tag)
    : ld_(x.ld_),
      data_(data),
      starts(x.starts),
      offset(x.offset, bw) { }

// VC++ doesn't like this
  //friend class transpose_type;

  //: The destructor.
  inline ~generic_dense2D() { }

  /* Access Methods */
  

  /* Iterator Access Methods */
  
  //: Return an iterator pointing to the first 1D container
  inline iterator begin() {
    return iterator(data(), ld_, 0, starts, offset);
  }
  //: Return an iterator pointing past the end of the 2D container
  inline iterator end() {
    return iterator(data(), ld_, offset.twod_length(), starts, offset);
  }
  //: Return a const iterator pointing to the first 1D container
  inline const_iterator begin() const {
    return const_iterator(data(), ld_, 0, starts, offset);
  }
  //: Return a const iterator pointing past the end of the 2D container
  inline const_iterator end() const {
    return const_iterator(data(), ld_, offset.twod_length(), 
                          starts, offset);
  }

  /* reverse iterators */

  //: Return a reverse iterator pointing to the last 1D container
  inline reverse_iterator rbegin() {
    return reverse_iterator(end());
  }
  //: Return a reverse iterator pointing past the start of the 2D container
  inline reverse_iterator rend() {
    return reverse_iterator(begin());
  }
  //: Return a const reverse iterator pointing to the last 1D container
  inline const_reverse_iterator rbegin() const {
    return const_reverse_iterator(end());
  }
  //: Return a const reverse iterator pointing past the start of the 2D container
  inline const_reverse_iterator rend() const {
    return const_reverse_iterator(begin());
  }


  /* Element Access Methods */
  //: Return a reference to the (i,j) element, where (i,j) is in the 2D coordinate system  
  inline const elt_type& operator()(size_type i, size_type j) const {
    return *(data() + offset.elt(i, j));
  }
  //: Return a const reference to the (i,j) element, where (i,j) is in the 2D coordinate system  
  inline elt_type& operator()(size_type i, size_type j) {
    return *(data() + offset.elt(i, j));
  }
  
  /* Size Methods */

  //: Number of non-zeroes
  inline size_type nnz() const { return offset.major() * offset.minor(); }

  //: Capacity
  inline size_type capacity() const { return offset.major() * offset.minor(); }
  
  //: Major axis size
  inline size_type major() const { return offset.major(); }

  //: Minor axis size
  inline size_type minor() const { return offset.minor(); }

  //: Leading Dimension
  inline size_type ld() const { return ld_; }
 
  //: Memory Access
  inline const elt_type* data() const { return &(*data_)[0]; }
  inline elt_type* data() { return &(*data_)[0]; }

  /* obsolete
  inline const elt_type* get_contiguous() const { return data_->data(); }
  inline elt_type* get_contiguous() { return data_->data(); }

  inline void set_contiguous(elt_type*) { }
  */

  /* Vector Access Methods */
  
  //: OneD Access
  inline OneD operator[](size_type i) const {
    typedef OneD* oned_ptr;
    typedef InnerOneD* inner_oned_ptr;
    return __bracket<is_strided>()((elt_type*)data() + offset.oned_offset(i),
                                   offset.oned_length(i), 
                                   starts.first, ld_,
                                   oned_ptr(), 
                                   inner_oned_ptr());
  }



  /* All the submatrix stuff is in matrix_implementation for now

  inline submatrix_type sub_matrix(size_type m_start, size_type m_finish,
                              size_type n_start, size_type n_finish) {
    return submatrix_type(data_->data() + m_start * ld_ + n_start,
                     dim_type(m_finish - m_start, n_finish - n_start), ld_);
  }
  inline submatrix_type sub_matrix(size_type m_start, size_type n_start,
                              size_type m, size_type n) {
    return submatrix_type(data_->data() + m_start * ld_ + n_start, 
                     dim_type(m, n), ld_);
  }

  inline submatrix_type section(size_type m_start, size_type n_start,
                           size_type m, size_type n) {
    return submatrix_type(data_->data() + m_start * ld_ + n_start, 
                     dim_type(m, n), ld_, dim_type(m_start, n_start));
  }
  typedef range<size_type> Range;
  inline generic_dense2D operator()(Range m, Range n) {
    return generic_dense2D(data_->data() + m.start * ld_ + n.start, 
                m.finish - m.start, n.finish - n.start, ld_);    
  }

  inline OneD::subrange_type operator()(size_type i, Range n) {
    return operator[i](n);
  }
  typedef strided1D< InnerOneD > MinorVector;
  inline MinorVector minor_vector(size_type i) const {
    InnerOneD vec((elt_type*)data_->data() + i, 
                  offset.major() * ld_, starts.first);
    return MinorVector(vec, ld_);
  }
  inline MinorVector::subrange_type operator()(Range m, size_type j) {
    return minor_vector(j)(m);
  }

  */

  /*JGS friend not working for transpose constructor  
  protected:
  */
  size_type ld_;/* JGS redundant */
  rep_ptr data_;
  pair_type starts;
  Offset offset;
};


template <class T, class OffsetGen, int M, int N>
struct gen_dense2D;

/* why didn't I use std::vector here?
 or perhaps I should just use plain old memory here?
 will that work with the reference counting in terms
 of deallocating?
 */

//: Dense2D Storage Type
//
// Inherits from generic_dense2D. The class "owns" its data.
//
//!category: containers
//!component: type
//!tparam: T - the element type
//!tparam: OffsetGen - the Offset class generator
//!tparam: MM - For static sized matrix, the major dimension
//!tparam: NN - For static sized matrix, the minor dimension
//!models: TwoDStorage

template<class T, class OffsetGen, int MM = 0, int NN = 0>
class dense2D
 : public generic_dense2D< std::vector<T> ,
               refcnt_ptr< std::vector<T> >, OffsetGen, MM, NN >
/* : public generic_dense2D< bare_bones_array<T> ,
               refcnt_ptr< bare_bones_array<T> >, OffsetGen, MM, NN >
*/
{
public:
  typedef generic_dense2D< std::vector<T> ,
               refcnt_ptr< std::vector<T> >, OffsetGen, MM, NN> super;
  /*  typedef generic_dense2D< bare_bones_array<T> ,
                 refcnt_ptr< bare_bones_array<T> >, OffsetGen, MM, NN> super;
  */
  typedef typename super::Offset Offset;
  //: Pair type for dimension
  typedef typename Offset::dim_type dim_type;
  //: Pair type for bandwidth
  typedef typename Offset::band_type band_type;
  typedef typename super::reptype reptype;
  typedef typename super::rep_ptr rep_ptr;
  //: Unsigned integral type for dimensions and indices
  typedef typename super::size_type size_type;
  //: The transpose type
  typedef dense2D<T, typename OffsetGen::transpose_type,
                  MM, NN> transpose_type;
// VC++ doesn't like this
  //friend class transpose_type;
  //: This has internal storage
  typedef internal_tag storage_loc;

  //: Default Constructor
  inline dense2D() { }

  //: Constructor from Dimension Pair
  inline dense2D(dim_type dim)
    : super(new reptype(Offset::size(dim.first(), dim.second(), 0, 0)),
            dim.first(), 
            dim.second(),
            dim.second()) { }

  //: Constructor from Dimension and Bandwidth Pairs
  inline dense2D(dim_type dim, band_type bw)
    : super(new reptype(Offset::size(dim.first(),dim.second(),
                                     bw.first(), bw.second())), 
            dim.first(), 
            dim.second(), 
            dim.second(),
            bw) { }

  //: Copy Constructor
  inline dense2D(const dense2D& x)
    : super(x) { }

  //: Assignment Operator
  inline dense2D& operator=(const dense2D& x) {
    super::operator=(x);
    return *this;
  }

  //: Transpose Constructor
  inline dense2D(const transpose_type& x, do_transpose t, do_transpose)
    : super(x, t, t) { }

#if !defined(_MSVCPP_)
  // JGS, use actual stream types
  //: Matrix Stream Constructor
  template <class MatrixStream, class Orien>
  inline dense2D(MatrixStream& s, Orien)
    : super(new reptype(Offset::size(Orien::map(dim_type(s.nrows(), 
                                                         s.ncols())).first(),
                                     Orien::map(dim_type(s.nrows(), 
                                                         s.ncols())).second(),
                                     0, 0)),
            s,
            Orien()) { }

  //: Matrix Stream Constructor with bandwidth
  template <class MatrixStream, class Orien>
  inline dense2D(MatrixStream& s, Orien, band_type bw)
    : super(new reptype(Offset::size(Orien::map(dim_type(s.nrows(), 
                                                         s.ncols())).first(),
                                     Orien::map(dim_type(s.nrows(), 
                                                         s.ncols())).second(),
                                     bw.first(), bw.second())),
            s,
            Orien(), 
            bw) { }
#endif

#if 0
  // deprecated
  template <class SubMatrix>
  struct partitioned {
    typedef dense2D<SubMatrix, OffsetGen> type;
    typedef gen_dense2D<SubMatrix, OffsetGen> generator;
  };
#endif
  
#if 1 // This makes no sense. dense2D can not be a "view"
  //: banded view constructor
  template <class TwoD>
  inline dense2D(const TwoD& x, band_type bw, banded_tag)
    : super(x.data_, x, bw, banded_tag()) { }
#endif

  //: Destructor
  inline ~dense2D() { }

  inline void resize(size_type m, size_type n) {
    rep_ptr newdata = new reptype(Offset::size(m, n, 0, 0));
    size_type i, j;
    size_type M = MTL_MIN(m, offset.major());
    size_type N = MTL_MIN(n, offset.minor());
    for (i = 0; i < M; ++i)
      for (j = 0; j < N; ++j)
	(*newdata)[i * n + j] = (*this)(i,j);
    for (; i < m; ++i)
      for (; j < n; ++j)
      (*newdata)[i * n + j] = T();

    data_ = newdata;
    ld_ = n;
    offset.dim = dim_type(m, n);
    offset.ld = n;
  }

};

template <class T, class OffsetGen, int M, int N>
struct gen_external2D; 

#ifndef MTL_DISABLE_BLOCKING
template <class Block, class OffsetGen, int M, int N>
struct gen_block2D;
#endif

//: blah
//!noindex:
template <class T, class OffsetGen, int M, int N>
struct gen_dense2D {
  typedef gen_dense2D<T, typename OffsetGen::transpose_type,N,M> transpose_type;
  typedef gen_external2D<T, OffsetGen,M,N> submatrix_type;

#ifndef MTL_DISABLE_BLOCKING
  template <class Block>
  struct blocked_view {
    typedef gen_block2D<Block, OffsetGen, M, N> type;
  };
#endif

  typedef gen_dense2D<T, typename OffsetGen::banded_view_type,M,N>
           banded_view_type;

  typedef dense2D<T, OffsetGen, M, N> type;
};


//: External2D Storage Type
//
// Inherits from generic_dense2D. The class does not "own" its data.
//
//!category: containers
//!component: type
//!tparam: T - the element type
//!tparam: OffsetGen - the Offset class generator
//!tparam: MM - For static sized matrix, the major dimension
//!tparam: NN - For static sized matrix, the minor dimension
//!models: TwoDStorage
//
template <class T, class OffsetGen, int MM = 0, int NN = 0>
class external2D
 : public generic_dense2D< external_vec<T,NN>,
                           external_vec<T,NN>*, OffsetGen, MM, NN >
{
  typedef generic_dense2D< external_vec<T,NN>,
                           external_vec<T,NN>*, OffsetGen, MM, NN > super;
public:
  external_vec<T,NN> rep;
  typedef dimension<T> dyn_dim;
  typedef typename super::Offset Offset;
  //: Pair type for dimension
  typedef typename Offset::dim_type dim_type;
  //: Pair type for bandwidth
  typedef typename Offset::band_type band_type;

  typedef typename super::reptype reptype;

  //: Unsigned integral type for dimensions and indices
  typedef typename super::size_type size_type;
  //: Type for the transpose
  typedef external2D<T, typename OffsetGen::transpose_type, 
                     MM, NN> transpose_type;
// VC++ doesn't like this
  //friend class transpose_type;
  //: This has external storage
  typedef external_tag storage_loc;

  //: Default Constructor
  inline external2D() { }

  //: Construct from pointer and dimensions
  inline external2D(T* data, dim_type dim)
    : super(&rep, dim.first(), dim.second(), dim.second()), 
      rep(data, dim.first() * dim.second()) 
  { }

  //: Construct from pointer, dimensions, and leading dimension
  inline external2D(T* data, dim_type dim, size_type ld)
    : super(&rep, dim.first(), dim.second(), ld),
      rep(data, dim.first() * ld)
  { }

  //: non-zero indices in upper left corner
  inline external2D(T* data, dim_type dim, size_type ld, 
		    dyn_dim s, char)
    : super(&rep, dim.first(), dim.second(), ld, s, char()),
      rep(data, dim.first() * ld)
  { }

  //: Constructor with bandwith
  inline external2D(T* data, dim_type dim, band_type bw)
    : super(&rep, dim.first(), dim.second(), dim.second(), bw),
      rep(data, dim.first() * dim.second())
  { }
  //: Constructor with leading dimension and bandwith
  inline external2D(T* data, dim_type dim, size_type ld, band_type bw)
    : super(&rep, dim.first(), dim.second(), ld, bw),
      rep(data, dim.first() * ld)
  { }

  //: Copy Constructor
  inline external2D(const external2D& x)
    :  super(&rep, x), rep(x.rep)
  { }

  //: Assignment Operator
  inline external2D& operator=(const external2D& x) {
    rep = x.rep;
    super::operator=(x);
    data_ = &rep;
    return *this;
  }

  //: Transpose Constructor
  inline external2D(const transpose_type& x, do_transpose t, do_transpose)
    : super(x, t, t), rep(x.rep) { }

  /* JGS This conflicts with the external2D(T* data, dim_type dim,
     band_type bw) constructor, and I am not sure this is really
     needed anyway.

  //: Matrix Stream Constructor
  template <class MatrixStream, class Orien>
  inline external2D(T* data, MatrixStream& s, Orien o)
    : rep(data, s.nrows() * s.ncols()), super(&rep, s, o) { }
  template <class MatrixStream, class Orien>
  inline external2D(T* data, MatrixStream& s, Orien o, 
                    band_type bw)
    : super(&rep, s, o, bw),
      rep(data, s.nrows() * s.ncols()) 
  { }  
  */

  //: banded view constructor
  template <class TwoD>
  inline external2D(const TwoD& x, band_type bw, banded_tag) 
    : super(&rep, x, bw, banded_tag()), rep((T*)x.data(), x.major() * x.ld()) { }


  inline ~external2D() { }
#if 0
  // deprecated
  template <class SubMatrix>
  struct partitioned {
    typedef dense2D<SubMatrix, OffsetGen> type;
    typedef gen_dense2D<SubMatrix, OffsetGen> generator;
  };
#endif
};

//: blah
//!noindex:
template <class T, class OffsetGen, int M, int N>
struct gen_external2D {
  typedef gen_external2D<T, typename OffsetGen::transpose_type, N, M> transpose_type;
  typedef gen_external2D<T, OffsetGen,M,N> submatrix_type;
  typedef gen_external2D<T, typename OffsetGen::banded_view_type,M,N>
           banded_view_type;

#ifndef MTL_DISABLE_BLOCKING
  template <class Block>
  struct blocked_view {
    typedef gen_block2D<Block, OffsetGen, M, N> type;
  };
#endif

  typedef external2D<T, OffsetGen, M, N> type;
};



} /* namespace mtl */


#endif /* MTL_DENSE2D_H */

#ifndef MTL_LIGHT_MATRIX_H
#define MTL_LIGHT_MATRIX_H

#include "mtl/matrix_traits.h"
#include "mtl/dimension.h"
#include "mtl/meta_if.h"
#include "mtl/meta_equal.h"

namespace mtl {


template <int Orien>
struct TRANS {
  enum { RET = 0 };
};

template<>
struct TRANS<ROW_MAJOR> {
  enum { RET = COL_MAJOR };
};

template<>
struct TRANS<COL_MAJOR> {
  enum { RET = ROW_MAJOR };
};


template <class T, class SizeType, int Orien, int Strided>
class light_matrix {
public:
  typedef light_matrix self;
  typedef light_matrix light_matrix_t; // VC++ workaround
  typedef T* DataPtr;

  typedef rectangle_tag shape;
  typedef typename IF< EQUAL<Orien,ROW_MAJOR>::RET,
              row_tag, column_tag>::RET orientation; // mostly wrong

  typedef typename IF< EQUAL<Orien,ROW_MAJOR>::RET,
              row_orien, column_orien>::RET orien;

  typedef light_matrix<T, SizeType, TRANS<Orien>::RET, Strided> transpose_type;
  typedef light_matrix<T, SizeType, Orien, !Strided> strided_type;
  typedef light_matrix<T, SizeType, Orien, Strided> scaled_type;// wrong

  typedef light_matrix<T, SizeType, Orien, Strided> submatrix_type;

  typedef int DiffType;

  //: The size type
  typedef SizeType size_type;
  //: The type for differences between iterators
  typedef DiffType difference_type;

  typedef T value_type;
  typedef value_type& reference;
  typedef const value_type& const_reference;
  typedef value_type* pointer;

  enum { M = 0, N = 0 };

protected:

  static inline size_type& twod_pos(size_type& i, size_type& j) {
    if (Orien == ROW_MAJOR)
	return i;
    else
	return j;
  }

  static inline const size_type& twod_pos(const size_type& i, 
					  const size_type& j) {
    if (Orien == ROW_MAJOR)
	return i;
    else
	return j;
  }

  static inline size_type& oned_pos(size_type& i, size_type& j) {
    if (Orien == ROW_MAJOR)
	return j;
    else
	return i;
  }

  static inline const size_type& oned_pos(const size_type& i, 
					  const size_type& j) {
    if (Orien == ROW_MAJOR)
	return j;
    else
	return i;
  }

  // idea: completely separate stride/offset/positioning from indexing
  //  but encapsulate both somehow

public:

  //: This is a dense 2D container
  typedef dense_tag sparsity;
  //: This has external storage
  typedef external_tag storage_loc;
  //: This is strideable
  typedef strideable strideability;

  class oned {
  public:
    typedef T& reference;
    typedef const T& const_reference;
    typedef T value_type;
    typedef T* pointer;
    typedef SizeType size_type;
    typedef int difference_type;

    enum { M = 0, N = 0 };

    typedef oned subrange_type;
    typedef dense_tag sparsity;
    typedef oned IndexArray; /* bogus */
    typedef oned IndexArrayRef; /* bogus */

    typedef oned_tag dimension; /* bogus */

    template <int isConst>
    class __iterator {
      typedef __iterator self;
    public:
      typedef typename oned::value_type value_type;
      typedef oned::pointer pointer;
      typedef typename oned::size_type size_type;
      typedef oned::difference_type difference_type;

      typedef typename IF<isConst, oned::const_reference, oned::reference>::RET reference;

      typedef std::random_access_iterator_tag iterator_category;

      inline __iterator(DataPtr d, 
			size_type ii, size_type jj,
			size_type os, size_type s)
	: data(d), i(ii), j(jj), offset(os), stride(s) { }

      inline __iterator(const self& x)
	: data(x.data), i(x.i), j(x.j), offset(x.offset), stride(x.stride) { }

      inline self& operator=(const self& x) {
	data = x.data; i = x.i; j = x.j; offset = x.offset; stride = x.stride;
	return *this;
      }

      inline __iterator() : data(0), i(0), j(0), offset(0), stride(0) { }

      inline reference operator*() const { return data[offset]; }
      inline self& operator++() { ++pos(); offset += stride; return *this; }
      inline self& operator+=(size_type n) { 
	pos() += n; offset += stride*n; return *this; 
      }
      inline self operator++(int) { self t = *this; ++(*this); return t; }
      inline self& operator--() { --pos(); offset -= stride; return *this; }
      inline self& operator-=(size_type n) { 
	pos() -= n; offset -= stride*n; return *this; }
      inline self operator--(int) { self t = *this; --(*this); return t; }
      inline bool operator!=(const self& x) const { return pos() != x.pos(); }
      inline bool operator==(const self& x) const { return pos() == x.pos(); }
      inline bool operator<(const self& x) const { return pos() < x.pos(); }
      inline size_type index() const { return pos(); }

      inline size_type& pos() { return oned_pos(i,j); }
      inline const size_type& pos() const { return oned_pos(i,j); }

      inline size_type row() const { return i; }
      inline size_type column() const { return j; }
    protected:
      DataPtr data;
      size_type i, j;
      size_type offset;
      size_type stride;
    };

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

    inline oned(DataPtr d, size_type ii, size_type jj, 
		size_type ie, size_type je,
		size_type os, size_type ld)
      : data(d), i(ii), j(jj), iend(ie), jend(je), 
	offset(os), ldim(ld) { }

    inline oned(const oned& x)
      : data(x.data), i(x.i), j(x.j), 
	iend(x.iend), jend(x.jend), 
	offset(x.offset), ldim(x.ldim) { }

    inline oned& operator=(const oned& x) {
      data = x.data; i = x.i; j = x.j; 
      iend = x.iend; jend = x.jend; 
      offset = x.offset; ldim = x.ldim;
      return *this;
    }
    inline oned()
      : data(0), i(0), j(0), iend(0), jend(0), offset(0), ldim(0) { }

    inline ~oned() { }

    inline reference operator[](size_type n) {
      return data[ Strided ? offset + n * ldim : offset + n];
    }

    inline const_reference operator[](size_type n) const {
      return data[ Strided ? offset + n * ldim : offset + n];
    }

    inline iterator begin() { 
      return iterator(data, i, j, offset, Strided ? ldim : 1); 
    }
    inline iterator end() { 
      size_type iiend, jjend;
      if (Orien == ROW_MAJOR) { iiend = i; jjend = jend; }
      else { iiend = iend; jjend = j; }

      return iterator(data, iiend, jjend, offset, Strided ? ldim: 1); 
    }

    inline const_iterator begin() const { 
      return const_iterator(data, i, j, offset, Strided ? ldim : 1); 
    }
    inline const_iterator end() const { 
      size_type iiend, jjend;
      if (Orien == ROW_MAJOR) { iiend = i; jjend = jend; }
      else { iiend = iend; jjend = j; }

      return const_iterator(data, iiend, jjend, offset, Strided ? ldim : 1); 
    }

  protected:
    DataPtr data;
    size_type i, j;
    size_type iend, jend;
    size_type offset;
    size_type ldim;
  };

  typedef oned OneD;
  typedef OneD OneDRef;
  typedef OneD Row;
  typedef OneD RowRef;
  typedef OneD Column;
  typedef OneD ColumnRef;

  //: The iterator type
  template <int Const>
  class __iterator {
    typedef __iterator self;
  public:
    typedef std::random_access_iterator_tag iterator_category;
    typedef oned value_type;
    typedef value_type* pointer;

#if defined(_MSVCPP_)
    typedef typename light_matrix_t::size_type size_type;
    typedef typename light_matrix_t::difference_type difference_type;
#else
    typedef SizeType size_type;
    typedef DiffType difference_type;
#endif

    typedef typename IF<Const, const oned, oned>::RET reference;
    
    inline __iterator(DataPtr d, size_type ii, size_type jj,
		      size_type ie, size_type je, size_type ld)
      : data(d), i(ii), j(jj), iend(ie), jend(je), offset(0), ldim(ld) { 
	if (Strided) stride = 1; else stride = ldim;
    }

    inline __iterator() : data(0), i(0), j(0), iend(0), jend(0),
	offset(0), stride(0), ldim(0) { }

    inline __iterator(const self& x)
      : data(x.data), i(x.i), j(x.j), 
	iend(x.iend), jend(x.jend), offset(x.offset), 
	stride(x.stride), ldim(x.ldim) { }

    inline self& operator=(const self& x) {
      data = x.data; i = x.i; j = x.j; 
      iend = x.iend; jend = x.jend; 
      offset = x.offset; stride = x.stride; ldim = x.ldim;
      return *this;
    }
    inline reference operator*() const {
      return oned(data, i, j, iend, jend, offset, ldim);
    }
    
    inline self& operator++() { ++pos(); offset += stride; return *this; }
    inline self& operator+=(size_type n) { 
      pos() += n; offset += stride*n; return *this; }
    inline self operator++(int) { self t = *this; ++(*this); return t; }
    inline self& operator--() { --pos(); offset -= stride; return *this; }
    inline self& operator-=(size_type n) { 
      pos() -= n; offset -= stride*n; return *this; }
    inline self operator--(int) { self t = *this; --(*this); return t; }
    inline bool operator!=(const self& x) const { return pos() != x.pos(); }
    inline bool operator==(const self& x) const { return pos() == x.pos(); }
    inline bool operator<(const self& x) const { return pos() < x.pos(); }
    inline size_type index() const { return pos(); }

    inline size_type& pos() { return twod_pos(i,j); }
    inline const size_type& pos() const { return twod_pos(i,j); }
    
    inline size_type row() const { return i; }
    inline size_type column() const { return j; }

  protected:
    DataPtr data;
    size_type i, j;
    size_type iend, jend;
    size_type offset;
    size_type stride;
    size_type ldim;
  };

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

  //: Standard Constructor
  inline light_matrix(DataPtr d, size_type m, size_type n, size_type ld)
    : data_(d), nrows_(m), ncols_(n), ldim(ld) { }

  inline light_matrix(DataPtr d, size_type m, size_type n)
    : data_(d), nrows_(m), ncols_(n), ldim(Orien == ROW_MAJOR ? n : m) { }

  //: Copy Constructor
  inline light_matrix(const light_matrix& x)
    : data_(x.data_), nrows_(x.nrows_), ncols_(x.ncols_), ldim(x.ldim) { }

  //: Assignment Operator
  inline const light_matrix& operator=(const light_matrix& x) {
    data_ = x.data_; nrows_ = x.nrows_; ncols_ = x.ncols_; ldim = x.ldim;
    return *this;
  }
  //: Default Constructor
  inline light_matrix() : data_(0), nrows_(0), ncols_(0), ldim(0) { }

  inline light_matrix(const strided_type& x, do_strided s)
    : data_(x.data_), nrows_(x.nrows_), ncols_(x.ncols_), ldim(x.ldim) { }

  template <class StridedType>
  inline light_matrix(const StridedType& x, do_strided s)
    : data_(x.data_), nrows_(x.nrows_), ncols_(x.ncols_), ldim(x.ldim) { }

  //: Destructor
  inline ~light_matrix() { }

  //: Return an iterator pointing to the first 1D container
  inline iterator begin() {
    return iterator(data_, 0, 0, nrows_, ncols_, ldim);
  }
  //: Return an iterator pointing past the end of the 2D container
  inline iterator end() {
    return iterator(data_, nrows_, ncols_, nrows_, ncols_, ldim);
  }

  //: Return a const iterator pointing to the first 1D container
  inline const_iterator begin() const {
    return const_iterator(data_, 0, 0, nrows_, ncols_, ldim);
  }
  //: Return a const iterator pointing past the end of the 2D container
  inline const_iterator end() const {
    return const_iterator(data_, nrows_, ncols_, nrows_, ncols_, ldim);
  }

  //: Return a reference to the ith 1D container
  inline oned operator[](size_type n) {
    if (Orien == ROW_MAJOR)
      return oned(data_, n, 0, nrows_, ncols_, Strided ? n : ldim * n, ldim);
    else
      return oned(data_, 0, n, nrows_, ncols_, Strided ? n : ldim * n, ldim);
  }

  inline const oned operator[](size_type n) const {
    if (Orien == ROW_MAJOR)
      return oned(data_, n, 0, nrows_, ncols_, Strided ? n : ldim * n, ldim);
    else
      return oned(data_, 0, n, nrows_, ncols_, Strided ? n : ldim * n, ldim);
  }

  //: Return a reference to the (i,j) element, where (i,j) is in the 2D coordinate system
  inline reference operator()(size_type i, size_type j) {
    return Orien == ROW_MAJOR ? operator[](i)[j] : operator[](j)[i];
  }

  //: Return a const reference to the (i,j) element, where (i,j) is in the 2D coordinate system  
  inline const_reference operator()(size_type i, size_type j) const {
    return Orien == ROW_MAJOR ? operator[](i)[j] : operator[](j)[i];
  }

  inline size_type nrows() const { return nrows_; }
  inline size_type ncols() const { return ncols_; }


  inline submatrix_type sub_matrix(size_type i, size_type iend,
				   size_type j, size_type jend) const
  {
    if (Strided)
      return submatrix_type(data_ + oned_pos(i,j) * ldim + twod_pos(iend,jend),
			    iend - i, jend - j, ldim);
    else
      return submatrix_type(data_ + twod_pos(i,j) * ldim + oned_pos(iend,jend),
			    iend - i, jend - j, ldim);
  }


  DataPtr data_;
  size_type nrows_, ncols_;
  size_type ldim;
};


} /* namespace mtl */

#endif /* MTL_LIGHT_MATRIX_H */

/* mtl/mtl_config.h.  Generated automatically by configure.  */
/* mtl/mtl_config.h.in.  Generated automatically from configure.in by autoheader.  */
#ifndef MTL_CONFIG_H
#define MTL_CONFIG_H


#define HAVE_COPYSIGN 1
#define HAVE_STD_ABS 1
#define MTL_EXCEPTIONS 0
/* #undef USE_BLAIS */
#define MTL_CMPLR_HAS_LIMITS 1

/* Define if you have the daxpy function.  */
/* #undef HAVE_DAXPY */

/* Define if you have the dgetrf function.  */
/* #undef HAVE_DGETRF */

/* Name of package */
#define PACKAGE "mtl"

/* Version number of package */
#define VERSION "2.1.2-21"


#if defined(_MSC_VER) && !defined(__MWERKS__) 
#	if	( _MSC_VER < 1300 )
#		define _MSVCPP_ _MSC_VER
#	elif(  _MSC_VER < 1310 )
#		define _MSVCPP7_ _MSC_VER
#	else
#		define _MSVCPP71_ _MSC_VER
#	endif
#endif

/* Assuming that configure is not run for Visual C++ and
   Metrowerks Codewarrior, so the above macros are not
   yet defined.
  */


#if !defined ( _MSVCPP_ )
#define STD_REVERSE_ITER 1
#else
#define STD_REVERSE_ITER 0
#endif


#if defined(__sgi) && !defined(__GNUC__) && (_COMPILER_VERSION < 730)
#  define MTL_NO_TEMPLATE
#elif defined(__GNUC__) && (__GNUC__ < 9)
#  define MTL_NO_TEMPLATE
#elif defined(__KCC_VERSION) && (__KCC_VERSION < 3400)
#  define MTL_NO_TEMPLATE
#elif defined(__MWERKS__)
#  define MTL_NO_TEMPLATE
#endif


#ifdef MTL_NO_TEMPLATE
#define MTL_TEMPLATE
#else
#define MTL_TEMPLATE template
#endif

#if defined ( _MSVCPP_ )
#define MTL_MAX(X,Y) (X > Y ? X : Y)
#define MTL_MIN(X,Y) (X < Y ? X : Y)
#else
#define MTL_MIN(X,Y) std::min(X,Y)
#define MTL_MAX(X,Y) std::max(X,Y)
#endif

#define MTL_CONJ std::conj
#define MTL_ABS std::abs


#if defined ( _MSVCPP_ ) || defined ( _MSVCPP7_ )
#define MTL_PARTIAL_SPEC 0
#else
#define MTL_PARTIAL_SPEC 1
#endif


#if defined(__GNUC__) || defined(_MSVCPP_) || defined(_MSVCPP7_) || defined( __MWERKS__)
#define MTL_DISABLE_BLOCKING
#endif

#if defined(MTL_FORTRAN_SYMBOLS_WITH_DOUBLE_TRAILING_UNDERSCORES)
#define MTL_FSYMB(X) X##__
#elif defined(MTL_FORTRAN_SYMBOLS_WITH_TRAILING_UNDERSCORES)
#define MTL_FSYMB(X) X##_
#else
#define MTL_FSYMB(X) X
#endif


#define PACKAGE "mtl"
#define VERSION "2.1.2-21"

#endif