/* Copyright (c) 1997-2006
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
   http://www.math.tu-berlin.de/polymake,  mailto:polymake@math.tu-berlin.de

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version: http://www.gnu.org/licenses/gpl.txt.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
*/

#ifndef _POLYMAKE_GENERIC_MATRIX_H
#define _POLYMAKE_GENERIC_MATRIX_H "$Project: polymake $$Id: GenericMatrix.h,v 1.38 2006/03/09 16:21:35 gawrilow Exp $"

#ifndef _POLYMAKE_VECTOR_H
#  include <Vector.h>
#endif
#ifndef _POLYMAKE_CASCADED_CONTAINER_H
#  include <CascadedContainer.h>
#endif
#ifndef _POLYMAKE_MATRIX_METHODS_H
#  include <matrix_methods.h>
#endif

namespace pm {

/* ------------------------------------------
 *  Matrix masquerade: ConcatRows, ConcatCols
 * ------------------------------------------ */

template <typename _Matrix>
class ConcatRows_default
   : public cascade_impl< ConcatRows_default<_Matrix>,
               cons< Container< Rows<_Matrix> >,
               cons< CascadeDepth< int2type<2> >,
		     Hidden< _Matrix > > >
            > {
   typedef cascade_impl<ConcatRows_default> _super;
public:
   int dim() const
   {
      return this->hidden().rows() * this->hidden().cols();
   }
   int size() const
   {
      return check_container_feature<_Matrix, sparse>::answer
	 ? _super::size() : dim();
   }
};

template <typename _Matrix>
class ConcatRows
   : public ConcatRows_default<_Matrix>,
     public GenericVector< ConcatRows<_Matrix>,
			   typename object_traits<typename _Matrix::value_type>::persistent_type > {
protected:
   ~ConcatRows();
public:
   ConcatRows& operator= (const ConcatRows& s) { return GenericVector<ConcatRows>::operator=(s); }
   using ConcatRows::generic_type::operator=;
   using ConcatRows_default<_Matrix>::dim;
};

template <typename _Matrix>
class ConcatCols
   : public ConcatRows< Transposed<_Matrix> > { };

template <typename _Matrix>
struct check_container_feature<ConcatRows<_Matrix>, sparse>
   : check_container_feature<_Matrix, sparse> { };

template <typename _Matrix>
struct check_container_feature<ConcatRows<_Matrix>, pure_sparse>
   : check_container_feature<_Matrix, pure_sparse> { };

template <typename _Matrix>
struct check_container_feature<ConcatCols<_Matrix>, sparse>
   : check_container_feature<_Matrix, sparse> { };

template <typename _Matrix>
struct check_container_feature<ConcatCols<_Matrix>, pure_sparse>
   : check_container_feature<_Matrix, pure_sparse> { };

template <typename _Matrix> inline
typename masquerade_matrix_helper<ConcatRows, _Matrix>::result_type&
concat_rows(_Matrix& m)
{
   return masquerade_matrix<ConcatRows>(m);
}

template <typename _Matrix> inline
typename masquerade_matrix_helper<ConcatRows, const _Matrix>::result_type&
concat_rows(const _Matrix& m)
{
   return masquerade_matrix<ConcatRows>(m);
}

template <typename _Matrix> inline
typename masquerade_matrix_helper<ConcatCols, _Matrix>::result_type&
concat_cols(_Matrix& m)
{
   return masquerade_matrix<ConcatCols>(m);
}

template <typename _Matrix> inline
typename masquerade_matrix_helper<ConcatCols, const _Matrix>::result_type&
concat_cols(const _Matrix& m)
{
   return masquerade_matrix<ConcatCols>(m);
}

/* ----------------
 *  Generic matrix
 * ---------------- */

template <typename E> class Matrix;
template <typename E> class SparseMatrix;
template <typename MatrixRef, typename Operation> class LazyMatrix1;
template <typename MatrixRef1, typename MatrixRef2, typename Operation> class LazyMatrix2;
template <typename E> class SparseMatrix2x2;

template <typename _Matrix, typename E=typename _Matrix::element_type>
class GenericMatrix
   : public Generic<_Matrix>,  public matrix_methods<_Matrix, E>,  public operators::base {
protected:
   GenericMatrix() { }
   GenericMatrix(const GenericMatrix&) { }
public:
   static const bool
      is_sparse= check_container_feature<_Matrix,sparse>::answer,
      prefer_concat=false;

   typedef typename if_else<is_sparse, SparseMatrix<E>, Matrix<E> >::result persistent_type;
   typedef GenericMatrix generic_type;
protected:
   bool is_shared() const { return false; }

   template <typename _Matrix2>
   void _assign(const GenericMatrix<_Matrix2>& M, True*)
   {
      concat_rows(*this)=concat_rows(M);
   }

   template <typename _Matrix2>
   void _assign(const GenericMatrix<_Matrix2>& M, False*)
   {
      copy(pm::rows(M).begin(), entire(pm::rows(*this)));
   }

   template <typename _Matrix2>
   void assign(const GenericMatrix<_Matrix2>& M)
   {
      DIMENSION_CHECK((M.rows() != this->rows() || M.cols() != this->cols()), "operator= - matrix dimension mismatch");
      _assign(M, (bool2type<_Matrix::prefer_concat && _Matrix2::prefer_concat>*)0);
   }

   template <typename Operation>
   void _assign_op(const Operation& op, True*)
   {
      concat_rows(*this).assign_op(op);
   }

   template <typename Operation>
   void _assign_op(const Operation& op, False*)
   {
      for (typename Entire< Rows<_Matrix> >::iterator r_i=entire(pm::rows(*this)); !r_i.at_end(); ++r_i)
	 r_i->assign_op(op);
   }

   template <typename _Matrix2, typename Operation>
   void _assign_op(const GenericMatrix<_Matrix2>& M, const Operation& op, True*)
   {
      concat_rows(*this).assign_op(concat_rows(M), op);
   }

   template <typename _Matrix2, typename Operation>
   void _assign_op(const GenericMatrix<_Matrix2>& M, const Operation& op, False*)
   {
      typename Rows<_Matrix2>::const_iterator m_i=pm::rows(M).begin();
      for (typename Entire< Rows<_Matrix> >::iterator r_i=entire(pm::rows(*this)); !r_i.at_end(); ++r_i, ++m_i)
	 r_i->assign_op(*m_i,op);
   }

   template <typename E2>
   void _fill(const E2& x, True*)
   {
      concat_rows(this->top()).fill(x);
   }

   template <typename E2>
   void _fill(const E2& x, False*)
   {
      for (typename Entire< Rows<_Matrix> >::iterator r_i=entire(pm::rows(*this)); !r_i.at_end(); ++r_i)
	 r_i->fill(x);
   }

   bool _zero_check(True*) const
   {
      return !concat_rows(*this);
   }

   bool _zero_check(False*) const
   {
      for (typename Entire< Rows<_Matrix> >::const_iterator r_i=entire(pm::rows(*this)); !r_i.at_end(); ++r_i)
	 if (*r_i) return false;
      return true;
   }

   template <typename _Matrix2>
   bool trivial_assignment(const GenericMatrix<_Matrix2,E>&) const { return false; }

   bool trivial_assignment(const GenericMatrix& M) const { return this==&M; }

public:
   template <typename Operation>
   void assign_op(const Operation& op)
   {
      if (this->top().is_shared())
	 this->top() = LazyMatrix1<const _Matrix&, Operation>(this->top(),op);
      else
	 _assign_op(op, (bool2type<_Matrix::prefer_concat>*)0);
   }

   template <typename _Matrix2, typename Operation>
   void assign_op(const GenericMatrix<_Matrix2>& M, const Operation& op)
   {
      if (this->top().is_shared())
	 this->top() = LazyMatrix2<const _Matrix&, const _Matrix2&, Operation>(this->top(),M.top(),op);
      else
	 _assign_op(M, op, (bool2type<_Matrix::prefer_concat && _Matrix2::prefer_concat>*)0);
   }

   _Matrix& operator= (const GenericMatrix& M)
   {
      if (!trivial_assignment(M))
	  this->top().assign(M);
      return this->top();
   }

   template <typename _Matrix2>
   _Matrix& operator= (const GenericMatrix<_Matrix2>& M)
   {
      this->top().assign(M);
      return this->top();
   }

   template <typename _Matrix2>
   void swap(GenericMatrix<_Matrix2,E>& M)
   {
      if (trivial_assignment(M)) return;
      DIMENSION_CHECK((this->rows() != M.rows() || this->cols() != M.cols()), "swap - matrix dimension mismatch");
      std::swap_ranges(entire(pm::rows(*this)), pm::rows(M).begin());
   }

   template <typename E2>
   void fill(const E2& x)
   {
      this->top()._fill(x, (bool2type<_Matrix::prefer_concat>*)0);
   }

   _Matrix& negate()
   {
      assign_op(BuildUnary<operations::neg>());
      return this->top();
   }

   friend
   _Matrix& inv_sign(GenericMatrix& me)
   {
      return me.negate();
   }

   friend
   _Matrix& inv_sign(const GenericMatrix& me)
   {
      return non_const(me.top()).negate();
   }

   template <typename _Matrix2>
   _Matrix& operator+= (const GenericMatrix<_Matrix2>& M)
   {
      DIMENSION_CHECK((this->rows() != M.rows() ||
		       this->cols() != M.cols()),  "operator+= - matrix dimension mismatch");
      assign_op(M, BuildBinary<operations::add>());
      return this->top();
   }

   template <typename _Matrix2>
   _Matrix& operator-= (const GenericMatrix<_Matrix2,E>& M)
   {
      DIMENSION_CHECK((this->rows() != M.rows() ||
		       this->cols() != M.cols()),  "operator-= - matrix dimension mismatch");
      assign_op(M, BuildBinary<operations::sub>());
      return this->top();
   }

   template <typename Right>
   _Matrix& operator*= (const Right& r)
   {
      if (r)
	 assign_op(operations::fix2<Right, BuildBinary<operations::mul> >(r));
      else
	 fill(0);
      return this->top();
   }

   template <typename Right>
   _Matrix& operator/= (const Right& r)
   {
      assign_op(operations::fix2<Right, BuildBinary<operations::div> >(r));
      return this->top();
   }

   /// Compare with an implicit zero matrix
   bool operator! () const
   {
      return _zero_check((bool2type<_Matrix::prefer_concat>*)0);
   }

   operator bool() const
   {
      return ! operator!();
   }

   /// @param i==0: main diagonal; i>0: i-th diagonal below the main; i<0: (-i)-th above the main
   IndexedSlice<ConcatRows<_Matrix>&, series>
   diagonal(int i=0)
   {
      const int
	 c=this->cols(),
	 start( i>0 ? i*c : -i ),
	 size( i>0 ? std::min(this->rows()-i, c) : std::min(this->rows(), c+i) );
      return IndexedSlice<ConcatRows<_Matrix>&, series> (concat_rows(*this), series(start,size,c+1));
   }

   const IndexedSlice<const ConcatRows<_Matrix>&, series>
   diagonal(int i=0) const
   {
      const int
	 c=this->cols(),
	 start( i>0 ? i*c : -i ),
	 size( i>0 ? std::min(this->rows()-i, c) : std::min(this->rows(), c+i) );
      return IndexedSlice<const ConcatRows<_Matrix>&, series> (concat_rows(*this), series(start,size,c+1));
   }

   /// @param i==0: main anti-diagonal; i>0: i-th diagonal below the main; i<0: (-i)-th above the main
   IndexedSlice<ConcatRows<_Matrix>&, series>
   anti_diagonal(int i=0)
   {
      const int
	 c=this->cols(),
	 start( i>0 ? (i+1)*c-1 : c+i-1 ),
	 size( i>0 ? std::min(this->rows()-i, c) : std::min(this->rows(), c+i) );
      return IndexedSlice<ConcatRows<_Matrix>&, series> (concat_rows(*this), series(start,size,c-1));
   }

   const IndexedSlice<const ConcatRows<_Matrix>&, series>
   anti_diagonal(int i=0) const
   {
      const int
	 c=this->cols(),
	 start( i>0 ? (i+1)*c-1 : c+i-1 ),
	 size( i>0 ? std::min(this->rows()-i, c) : std::min(this->rows(), c+i) );
      return IndexedSlice<const ConcatRows<_Matrix>&, series> (concat_rows(*this), series(start,size,c-1));
   }

protected:
   template <typename Line, typename E2>
   static void _multiply(Line& l_i, Line& l_j,
			 const E2& a_ii, const E2& a_ij, const E2& a_ji, const E2& a_jj, False*)
   {
      typename Line::iterator e_j=l_j.begin();
      for (typename Entire<Line>::iterator e_i=entire(l_i); !e_i.at_end(); ++e_i, ++e_j) {
	 const E x_i= (*e_i)*a_ii + (*e_j)*a_ij;
	 *e_j =       (*e_i)*a_ji + (*e_j)*a_jj;
	 *e_i =  x_i;
      }
   }

   template <typename Line, typename E2>
   static void _multiply(Line& l_i, Line& l_j,
			 const E2& a_ii, const E2& a_ij, const E2& a_ji, const E2& a_jj, True*)
   {
      typename Line::iterator e_i=l_i.begin(), e_j=l_j.begin();
      int state=zipper_both;
      if (e_i.at_end()) state>>=3;
      if (e_j.at_end()) state>>=6;
      while (state) {
	 if (state >= zipper_both) {
	    state &= ~zipper_cmp;
	    state += 1 << sign(e_i.index() - e_j.index())+1;
	 }
	 if (state & zipper_lt) {
	    if (a_ji)
	       l_j.insert(e_j,e_i.index(),(*e_i)*a_ji);
	    if (a_ii)
	       *e_i++ *= a_ii;
	    else
	       l_i.erase(e_i++);
	    if (e_i.at_end())
	       state>>=3;
	 } else if (state & zipper_gt) {
	    if (a_ij)
	       l_i.insert(e_i,e_j.index(),(*e_j)*a_ij);
	    if (a_jj)
	       *e_j++ *= a_jj;
	    else
	       l_j.erase(e_j++);
	    if (e_j.at_end())
	       state>>=6;
	 } else {
	    const E x_i= (*e_i)*a_ii + (*e_j)*a_ij;
	          *e_j = (*e_i)*a_ji + (*e_j)*a_jj;
	    if (x_i)
	       *e_i++ = x_i;
	    else
	       l_i.erase(e_i++);
	    if (e_i.at_end())
	       state>>=3;
	    if (*e_j)
	       ++e_j;
	    else
	       l_j.erase(e_j++);
	    if (e_j.at_end())
	       state>>=6;
	 }
      }
   }

public:
   void multiply_from_left(const SparseMatrix2x2<E>& U)
   {
      DIMENSION_CHECK((U.i<0 || U.j<0 || U.i >= this->rows() || U.j >= this->rows()),
		      "multiply_from_left - dimension mismatch");
      _multiply(non_const(this->row(U.i)), non_const(this->row(U.j)),
		U.a_ii, U.a_ij, U.a_ji, U.a_jj, (bool2type<is_sparse>*)0);
   }

   void multiply_from_left(const Transposed< SparseMatrix2x2<E> >& U)
   {
      DIMENSION_CHECK((U.i<0 || U.j<0 || U.i >= this->rows() || U.j >= this->rows()),
		      "multiply_from_left - dimension mismatch");
      _multiply(non_const(this->row(U.i)), non_const(this->row(U.j)),
		U.a_ii, U.a_ji, U.a_ij, U.a_jj, (bool2type<is_sparse>*)0);
   }

   void multiply_from_right(const SparseMatrix2x2<E>& U)
   {
      DIMENSION_CHECK((U.i<0 || U.j<0 || U.i >= this->cols() || U.j >= this->cols()),
		      "multiply_from_right - dimension mismatch");
      _multiply(non_const(this->col(U.i)), non_const(this->col(U.j)),
		U.a_ii, U.a_ji, U.a_ij, U.a_jj, (bool2type<is_sparse>*)0);
   }

   void multiply_from_right(const Transposed< SparseMatrix2x2<E> >& U)
   {
      DIMENSION_CHECK((U.i<0 || U.j<0 || U.i >= this->cols() || U.j >= this->cols()),
		      "multiply_from_right - dimension mismatch");
      _multiply(non_const(this->col(U.i)), non_const(this->col(U.j)),
		U.a_ii, U.a_ij, U.a_ji, U.a_jj, (bool2type<is_sparse>*)0);
   }

   template <typename Result>
   struct rebind_generic {
      typedef GenericMatrix<Result, E> type;
   };
};

struct is_matrix;

template <typename Matrix, typename E>
struct spec_object_traits< GenericMatrix<Matrix,E> >
   : spec_or_model_traits<Matrix,is_container> {
   static const int dimension=2, is_resizeable=2;
   static const bool allow_sparse=true;
   typedef is_matrix generic_tag;

   typedef Rows<Matrix> serialized;
   static serialized& serialize(Matrix& M) { return rows(M); }
   static const serialized& serialize(const Matrix& M) { return rows(M); }
};

/* --------------------------------------------
 *  LazyMatrix1
 * lazy evaluation of an unary matrix operator
 * -------------------------------------------- */

template <typename MatrixRef, typename Operation>
struct lazy1_traits {
   typedef unary_op_builder<Operation, void, typename deref<MatrixRef>::type::const_reference>
      op_builder;
   typedef typename op_builder::operation::result_type reference;
   typedef typename deref<reference>::type value_type;
   typedef typename object_traits<value_type>::persistent_type element_type;
};

template <typename MatrixRef, typename Operation>
class LazyMatrix1
   : public modified_container_base<MatrixRef, Operation>,
     public GenericMatrix< LazyMatrix1<MatrixRef,Operation>,
			   typename lazy1_traits<MatrixRef,Operation>::element_type> {
   typedef modified_container_base<MatrixRef, Operation> _base;
public:
   typedef typename lazy1_traits<MatrixRef,Operation>::value_type value_type;
   typedef typename lazy1_traits<MatrixRef,Operation>::reference reference;
   typedef reference const_reference;
   typedef typename deref<MatrixRef>::type matrix_type;
   static const bool prefer_concat=matrix_type::prefer_concat;

   LazyMatrix1(typename _base::arg_type src_arg, const Operation& op_arg=Operation())
      : _base(src_arg, op_arg) { }
};

template <typename MatrixRef, typename Operation>
struct spec_object_traits< LazyMatrix1<MatrixRef, Operation> >
   : spec_object_traits<is_container> {
   static const bool is_lazy=true, is_temporary=true, is_always_const=true;
};

template <typename MatrixRef, typename Operation>
struct check_container_feature<LazyMatrix1<MatrixRef, Operation>, sparse>
   : check_container_ref_feature<MatrixRef, sparse> { };

template <typename MatrixRef, typename Operation>
struct check_container_feature<LazyMatrix1<MatrixRef, Operation>, pure_sparse>
   : check_container_ref_feature<MatrixRef, pure_sparse> { };

template <typename MatrixRef, typename Operation>
class matrix_random_access_methods< LazyMatrix1<MatrixRef, Operation> > {
   typedef LazyMatrix1<MatrixRef,Operation> master;
public:
   typename lazy1_traits<MatrixRef,Operation>::reference
   operator() (int i, int j) const
   {
      const master& me=static_cast<const master&>(*this);
      typedef typename lazy1_traits<MatrixRef,Operation>::op_builder opb;
      return opb::create(me.get_accessor()) (me.get_container()(i,j));
   }
};

template <typename MatrixRef, typename Operation>
class ConcatRows< LazyMatrix1<MatrixRef, Operation> >
   : public LazyVector1<masquerade<ConcatRows, MatrixRef>, Operation> {
protected:
   ConcatRows();
   ~ConcatRows();
};

template <typename MatrixRef, typename Operation>
class Rows< LazyMatrix1<MatrixRef,Operation> >
   : public modified_container_impl< Rows< LazyMatrix1<MatrixRef,Operation> >,
               cons< Container< masquerade<pm::Rows,MatrixRef> >,
	       cons< Accessor< operations::construct_unary2<LazyVector1,Operation> >,
                     Hidden< LazyMatrix1<MatrixRef,Operation> > > >
            > {
   typedef modified_container_impl<Rows> _super;
protected:
   ~Rows();
public:
   const typename _super::container& get_container() const
   {
      return rows(this->hidden().get_container());
   }
   typename _super::accessor get_accessor() const
   {
      return this->hidden().get_accessor();
   }
};

template <typename MatrixRef, typename Operation>
class Cols< LazyMatrix1<MatrixRef,Operation> >
   : public modified_container_impl< Cols< LazyMatrix1<MatrixRef,Operation> >,
               cons< Container< masquerade<pm::Cols,MatrixRef> >,
	       cons< Accessor< operations::construct_unary2<LazyVector1,Operation> >,
                     Hidden< LazyMatrix1<MatrixRef,Operation> > > >
            > {
   typedef modified_container_impl<Cols> _super;
protected:
   ~Cols();
public:
   const typename _super::container& get_container() const
   {
      return cols(this->hidden().get_container());
   }
   typename _super::accessor get_accessor() const
   {
      return this->hidden().get_accessor();
   }
};

/* --------------------------------------------
 *  LazyMatrix2
 * lazy evaluation of a binary matrix operator
 * -------------------------------------------- */

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
struct lazy2_traits {
   typedef binary_op_builder<Operation, void, void,
			     typename deref<MatrixRef1>::type::const_reference,
			     typename deref<MatrixRef2>::type::const_reference>
      op_builder;
   typedef typename op_builder::operation::result_type reference;
   typedef typename deref<reference>::type value_type;
   typedef typename object_traits<value_type>::persistent_type element_type;
};

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
class LazyMatrix2
   : public modified_container_pair_base<MatrixRef1, MatrixRef2, Operation>,
     public GenericMatrix< LazyMatrix2<MatrixRef1,MatrixRef2,Operation>,
                           typename lazy2_traits<MatrixRef1,MatrixRef2,Operation>::element_type > {
   typedef modified_container_pair_base<MatrixRef1, MatrixRef2, Operation> _base;
public:
   typedef typename lazy2_traits<MatrixRef1,MatrixRef2,Operation>::value_type value_type;
   typedef typename lazy2_traits<MatrixRef1,MatrixRef2,Operation>::reference reference;
   typedef reference const_reference;
   typedef typename deref<MatrixRef1>::type matrix1_type;
   typedef typename deref<MatrixRef2>::type matrix2_type;
   static const bool prefer_concat=matrix1_type::prefer_concat && matrix2_type::prefer_concat;

   LazyMatrix2(typename _base::first_arg_type src1_arg, typename _base::second_arg_type src2_arg, const Operation& op_arg=Operation())
      : _base(src1_arg, src2_arg, op_arg) { }
};

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
struct spec_object_traits< LazyMatrix2<MatrixRef1, MatrixRef2, Operation> >
   : spec_object_traits<is_container> {
   static const bool is_lazy=true, is_temporary=true, is_always_const=true;
};

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
struct check_container_feature<LazyMatrix2<MatrixRef1, MatrixRef2, Operation>, sparse> {
   static const bool answer=TransformedContainerPair_helper1<MatrixRef1, MatrixRef2>::sparse_result;
};

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
struct check_container_feature<LazyMatrix2<MatrixRef1, MatrixRef2, Operation>, dense> {
   static const bool answer=TransformedContainerPair_helper1<MatrixRef1, MatrixRef2>::dense_result;
};

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
class matrix_random_access_methods< LazyMatrix2<MatrixRef1, MatrixRef2, Operation> > {
   typedef LazyMatrix2<MatrixRef1,MatrixRef2,Operation> master;
public:
   typename lazy2_traits<MatrixRef1,MatrixRef2,Operation>::reference
   operator() (int i, int j) const
   {
      const master& me=static_cast<const master&>(*this);
      typedef typename lazy2_traits<MatrixRef1,MatrixRef2,Operation>::op_builder opb;
      return opb::create(me.get_accessor()) (me.get_container1()(i,j), me.get_container2()(i,j));
   }
};

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
class ConcatRows< LazyMatrix2<MatrixRef1, MatrixRef2, Operation> >
   : public LazyVector2<masquerade<ConcatRows,MatrixRef1>, masquerade<ConcatRows,MatrixRef2>, Operation> {
protected:
   ConcatRows();
   ~ConcatRows();
};

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
class Rows< LazyMatrix2<MatrixRef1, MatrixRef2, Operation> >
   : public modified_container_pair_impl< Rows< LazyMatrix2<MatrixRef1,MatrixRef2,Operation> >,
               cons< Container1< masquerade<pm::Rows,MatrixRef1> >,
               cons< Container2< masquerade<pm::Rows,MatrixRef2> >,
	       cons< Accessor< operations::construct_binary2<LazyVector2,Operation> >,
		     Hidden< LazyMatrix2<MatrixRef1,MatrixRef2,Operation> > > > >
            > {
   typedef modified_container_pair_impl<Rows> _super;
protected:
   ~Rows();
public:
   const typename _super::container1& get_container1() const
   {
      return rows(this->hidden().get_container1());
   }
   const typename _super::container2& get_container2() const
   {
      return rows(this->hidden().get_container2());
   }
   typename _super::accessor get_accessor() const
   {
      return this->hidden().get_accessor();
   }
};

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
class Cols< LazyMatrix2<MatrixRef1, MatrixRef2, Operation> >
   : public modified_container_pair_impl< Cols< LazyMatrix2<MatrixRef1,MatrixRef2,Operation> >,
               cons< Container1< masquerade<pm::Cols,MatrixRef1> >,
               cons< Container2< masquerade<pm::Cols,MatrixRef2> >,
               cons< Accessor< operations::construct_binary2<LazyVector2,Operation> >,
		     Hidden< LazyMatrix2<MatrixRef1,MatrixRef2,Operation> > > > >
            > {
   typedef modified_container_pair_impl<Cols> _super;
protected:
   ~Cols();
public:
   const typename _super::container1& get_container1() const
   {
      return cols(this->hidden().get_container1());
   }
   const typename _super::container2& get_container2() const
   {
      return cols(this->hidden().get_container2());
   }
   typename _super::accessor get_accessor() const
   {
      return this->hidden().get_accessor();
   }
};

template <typename MatrixRef, typename Operation>
class ConcatRows< Transposed< LazyMatrix1<MatrixRef, Operation> > >
   : public ConcatRows< LazyMatrix1<masquerade<Transposed,MatrixRef>, Operation> > { };

template <typename MatrixRef1, typename MatrixRef2, typename Operation>
class ConcatRows< Transposed< LazyMatrix2<MatrixRef1, MatrixRef2, Operation> > >
   : public ConcatRows< LazyMatrix2<masquerade<Transposed,MatrixRef1>,
                                    masquerade<Transposed,MatrixRef2>, Operation> > { };

/// explicit conversion of matrix elements to another type
template <typename TargetType, typename _Matrix> inline
const LazyMatrix1<const _Matrix&, conv<typename _Matrix::element_type, TargetType> >
convert_to(const GenericMatrix<_Matrix>& M)
{
   return M.top();
}

/* --------------------------
 *  MatrixMinor concatenated
 * -------------------------- */

template <typename MatrixRef, typename RowIndexSetRef, typename ColIndexSetRef,
          bool contiguous=derived_from<typename deref<RowIndexSetRef>::type, sequence>::answer>
class MatrixMinorConcatRows
   : public ConcatRows_default< MatrixMinor<MatrixRef, RowIndexSetRef, ColIndexSetRef> > { };

template <typename MatrixRef, typename RowIndexSetRef>
class MatrixMinorConcatRows<MatrixRef, RowIndexSetRef, const all_selector&, true>
   : public indexed_subset_impl< MatrixMinorConcatRows<MatrixRef, RowIndexSetRef, const all_selector&, true>,
               cons<Container1< masquerade<ConcatRows,MatrixRef> >,
               cons<Container2< sequence >,
                    Hidden< MatrixMinor<MatrixRef, RowIndexSetRef, const all_selector&> > > >
            > {
   typedef indexed_subset_impl<MatrixMinorConcatRows> _super;
public:
   typename _super::container1& get_container1()
   {
      return concat_rows(this->hidden().get_matrix());
   }
   const typename _super::container1& get_container1() const
   {
      return concat_rows(this->hidden().get_matrix());
   }
   sequence get_container2() const
   {
      return sequence(this->hidden().get_subset((int2type<1>*)0).front() * this->hidden().cols(), dim());
   }

   int dim() const { return this->hidden().rows() * this->hidden().cols(); }
};

template <typename MatrixRef, typename RowIndexSetRef, typename ColIndexSetRef>
class ConcatRows< MatrixMinor<MatrixRef, RowIndexSetRef, ColIndexSetRef> >
   : public MatrixMinorConcatRows<MatrixRef,RowIndexSetRef,ColIndexSetRef>,
     public GenericVector< ConcatRows< MatrixMinor<MatrixRef,RowIndexSetRef,ColIndexSetRef> >,
                           typename deref<MatrixRef>::type::element_type > {
protected:
   ~ConcatRows();
public:
   int dim() const
   {
      return this->hidden().rows() * this->hidden().cols();
   }

   ConcatRows& operator= (const ConcatRows& s)
   {
      return ConcatRows::generic_type::operator=(s);
   }
   using ConcatRows::generic_type::operator=;
};

/* -----------------------------
 *  Block matrices concatenated
 * ----------------------------- */

template <typename MatrixRef1, typename MatrixRef2>
class ConcatRows< RowChain<MatrixRef1, MatrixRef2> >
   : public container_chain_impl< ConcatRows< RowChain<MatrixRef1,MatrixRef2> >,
               cons< Container1< masquerade<pm::ConcatRows,MatrixRef1> >,
               cons< Container2< masquerade<pm::ConcatRows,MatrixRef2> >,
                    Hidden< RowChain<MatrixRef1,MatrixRef2> > > >
            >,
     public GenericVector< ConcatRows< RowChain<MatrixRef1,MatrixRef2> >,
                                       typename deref<MatrixRef1>::type::element_type> {
   typedef container_chain_impl<ConcatRows> _super;
protected:
   ~ConcatRows();
public:
   typename ConcatRows::container1& get_container1()
   {
      return concat_rows(this->hidden().get_container1());
   }
   typename ConcatRows::container2& get_container2()
   {
      return concat_rows(this->hidden().get_container2());
   }
   const typename ConcatRows::container1& get_container1() const
   {
      return concat_rows(this->hidden().get_container1());
   }
   const typename ConcatRows::container2& get_container2() const
   {
      return concat_rows(this->hidden().get_container2());
   }

   ConcatRows& operator= (const ConcatRows& s)
   {
      return ConcatRows::generic_type::operator=(s);
   }
   using ConcatRows::generic_type::operator=;
   using _super::dim;
};

template <typename MatrixRef1, typename MatrixRef2>
class ConcatRows< Transposed< RowChain<MatrixRef1,MatrixRef2> > >
   : public ConcatRows< ColChain<masquerade<Transposed,MatrixRef1>, masquerade<Transposed,MatrixRef2> > > { };

template <typename MatrixRef1, typename MatrixRef2>
class ConcatRows< Transposed< ColChain<MatrixRef1,MatrixRef2> > >
   : public ConcatRows< RowChain<masquerade<Transposed,MatrixRef1>, masquerade<Transposed,MatrixRef2> > > { };

/* ----------------------
 *  SingleRow, SingleCol
 * ---------------------- */

template <typename _Vector>
class single_line_matrix {
protected:
   ~single_line_matrix() { }
public:
   typedef typename _Vector::value_type value_type;
   typedef typename _Vector::reference reference;
   typedef typename _Vector::const_reference const_reference;

   _Vector& hidden() { return *reinterpret_cast<_Vector*>(this); }
   const _Vector& hidden() const { return *reinterpret_cast<const _Vector*>(this); }
};

template <typename _Vector>
class SingleRow
   : public single_line_matrix<_Vector>,
     public GenericMatrix< SingleRow<_Vector>, typename _Vector::element_type> {
public:
   static const bool prefer_concat=true;
};

template <typename _Vector>
class SingleCol
   : public single_line_matrix<_Vector>,
     public GenericMatrix< SingleCol<_Vector>, typename _Vector::element_type> {
public:
   static const bool prefer_concat=true;
};

template <typename _Vector>
class matrix_random_access_methods< SingleRow<_Vector> > {
public:
   typename _Vector::reference
   operator() (int i, int j)
   {
      SingleRow<_Vector>& me=*static_cast<SingleRow<_Vector>*>(this);
      INDEX_CHECK((i!=0), "operator() - row index out of range");
      INDEX_CHECK((j<0 || j>=me.hidden().dim()), "operator() - column index out of range");
      return me.hidden()[j];
   }

   typename _Vector::const_reference
   operator() (int i, int j) const
   {
      const SingleRow<_Vector>& me=*static_cast<const SingleRow<_Vector>*>(this);
      INDEX_CHECK((i!=0), "operator() - row index out of range");
      INDEX_CHECK((j<0 || j>=me.hidden().dim()), "operator() - column index out of range");
      return me.hidden()[j];
   }
};

template <typename _Vector>
class matrix_random_access_methods< SingleCol<_Vector> > {
public:
   typename _Vector::reference
   operator() (int i, int j)
   {
      SingleCol<_Vector>& me=*static_cast<SingleCol<_Vector>*>(this);
      INDEX_CHECK((i<0 || i>=me.get_vector().dim()), "operator() - row index out of range");
      INDEX_CHECK((j!=0), "operator() - column index out of range");
      return me.hidden()[i];
   }

   typename _Vector::const_reference
   operator() (int i, int j) const
   {
      const SingleCol<_Vector>& me=*static_cast<const SingleCol<_Vector>*>(this);
      INDEX_CHECK((i<0 || i>=me.get_vector().dim()), "operator() - row index out of range");
      INDEX_CHECK((j!=0), "operator() - column index out of range");
      return me.hidden()[i];
   }
};

template <typename _Vector>
class ConcatRows< SingleRow<_Vector> > : public _Vector {
protected:
   ConcatRows();
   ~ConcatRows();
};

template <typename _Vector>
class ConcatRows< SingleCol<_Vector> > : public _Vector {
protected:
   ConcatRows();
   ~ConcatRows();
};

template <typename _Vector>
struct check_container_feature< SingleRow<_Vector>, sparse>
   : check_container_feature<_Vector, sparse> { };

template <typename _Vector>
struct check_container_feature< SingleRow<_Vector>, pure_sparse>
   : check_container_feature<_Vector, pure_sparse> { };

template <typename _Vector>
struct check_container_feature< SingleCol<_Vector>, sparse>
   : check_container_feature<_Vector, sparse> { };

template <typename _Vector>
struct check_container_feature< SingleCol<_Vector>, pure_sparse>
   : check_container_feature<_Vector, pure_sparse> { };

template <typename _Vector, bool _sparse=check_container_feature<_Vector,sparse>::answer>
class single_line_across
   : public modified_container_impl< single_line_across<_Vector,_sparse>,
               cons< Accessor< operations::construct_or_masquerade<SingleElementVector> >,
		     Hidden< _Vector > >
            > {
public:
   int size() const { return this->hidden().size(); }
};

template <typename _Vector>
class single_line_across<_Vector, true>
   : public modified_container_pair_impl< single_line_across<_Vector,true>,
               cons< Container1< _Vector >,
               cons< Container2< sequence >,
               cons< Coupler< zipping_coupler< operations::cmp<>, set_union_zipper, true, false> >,
               cons< Accessor< BuildBinary<SingleElementSparseVector_factory> >,
                     Hidden< _Vector > > > > >
            > {
public:
   sequence get_container2() const { return sequence(0, size()); }
   int size() const { return this->hidden().dim(); }
};

template <typename _Vector>
class Rows< SingleRow<_Vector> > : public single_elem_container<_Vector&> {
protected:
   ~Rows();
};
template <typename _Vector>
class Cols< SingleRow<_Vector> > : public single_line_across<_Vector> {
protected:
   ~Cols();
};
template <typename _Vector>
class Rows< SingleCol<_Vector> > : public single_line_across<_Vector> {
protected:
   ~Rows();
};
template <typename _Vector>
class Cols< SingleCol<_Vector> > : public single_elem_container<_Vector&> {
protected:
   ~Cols();
};

template <typename _Vector>
class ConcatRows< Transposed< SingleRow<_Vector> > > : public _Vector {
protected:
   ConcatRows();
   ~ConcatRows();
};

template <typename _Vector>
class ConcatRows< Transposed< SingleCol<_Vector> > > : public _Vector {
protected:
   ConcatRows();
   ~ConcatRows();
};

template <typename _Vector> inline
SingleRow<_Vector>&
vector2row(GenericVector<_Vector>& v)
{
   return reinterpret_cast<SingleRow<_Vector>&>(v.top());
}

template <typename _Vector> inline
const SingleRow<_Vector>&
vector2row(const GenericVector<_Vector>& v)
{
   return reinterpret_cast<const SingleRow<_Vector>&>(v.top());
}

template <typename _Vector> inline
SingleCol<_Vector>&
vector2col(GenericVector<_Vector>& v)
{
   return reinterpret_cast<SingleCol<_Vector>&>(v.top());
}

template <typename _Vector> inline
const SingleCol<_Vector>&
vector2col(const GenericVector<_Vector>& v)
{
   return reinterpret_cast<const SingleCol<_Vector>&>(v.top());
}

/* --------------------------
 *  RepeatedRow, RepeatedCol
 * -------------------------- */

template <typename VectorRef>
class repeated_line_matrix {
protected:
   alias<VectorRef> vector;
   int cnt;
public:
   typedef typename deref<VectorRef>::type::value_type value_type;
   typedef typename deref<VectorRef>::type::const_reference const_reference;
   typedef const_reference reference;

   typedef typename alias<VectorRef>::arg_type arg_type;
   repeated_line_matrix(arg_type vector_arg, int cnt_arg) : vector(vector_arg), cnt(cnt_arg) { }

   typename alias<VectorRef>::const_reference get_vector() const { return vector; }
   int get_count() const { return cnt; }
};

template <typename VectorRef,
          bool _same_elem=derived_from_instance<typename deref<VectorRef>::type, repeated_item_container>::answer>
class RepeatedRow
   : public repeated_line_matrix<VectorRef>,
     public GenericMatrix< RepeatedRow<VectorRef,_same_elem>,
			   typename deref<VectorRef>::type::element_type> {
   typedef repeated_line_matrix<VectorRef> _base;
public:
   static const bool prefer_concat=_same_elem;

   RepeatedRow(typename _base::arg_type vector_arg, int cnt_arg)
      : _base(vector_arg, cnt_arg) { }
};

template <typename VectorRef,
          bool _same_elem=derived_from_instance<typename deref<VectorRef>::type, repeated_item_container>::answer>
class RepeatedCol
   : public repeated_line_matrix<VectorRef>,
     public GenericMatrix< RepeatedCol<VectorRef,_same_elem>,
			   typename deref<VectorRef>::type::element_type> {
   typedef repeated_line_matrix<VectorRef> _base;
public:
   static const bool prefer_concat=_same_elem;

   RepeatedCol(typename _base::arg_type vector_arg, int cnt_arg)
      : _base(vector_arg, cnt_arg) { }
};

template <typename VectorRef, bool _same_elem>
struct spec_object_traits< RepeatedRow<VectorRef, _same_elem> >
   : spec_object_traits<is_container> {
   static const bool is_temporary=true, is_always_const=true;
};

template <typename VectorRef, bool _same_elem>
struct spec_object_traits< RepeatedCol<VectorRef, _same_elem> >
   : spec_object_traits<is_container> {
   static const bool is_temporary=true, is_always_const=true;
};

template <typename VectorRef, bool _same_elem>
class matrix_random_access_methods< RepeatedRow<VectorRef,_same_elem> > {
   typedef RepeatedRow<VectorRef,_same_elem> master;
public:
   typename container_traits<VectorRef>::const_reference
   operator() (int i, int j) const
   {
      const master& me=static_cast<const master&>(*this);
      INDEX_CHECK((i<0 || i>=me.rows()), "operator() - row index out of range");
      INDEX_CHECK((j<0 || j>=me.cols()), "operator() - column index out of range");
      return me.get_vector()[j];
   }
};

template <typename VectorRef, bool _same_elem>
class matrix_random_access_methods< RepeatedCol<VectorRef,_same_elem> > {
   typedef RepeatedCol<VectorRef,_same_elem> master;
public:
   typename container_traits<VectorRef>::const_reference
   operator() (int i, int j) const
   {
      const master& me=static_cast<const master&>(*this);
      INDEX_CHECK((i<0 || i>=me.rows()), "operator() - row index out of range");
      INDEX_CHECK((j<0 || j>=me.cols()), "operator() - column index out of range");
      return me.get_vector()[i];
   }
};

template <typename VectorRef>
class ConcatRows< RepeatedRow<VectorRef, true> >
   : public constant_container< ConcatRows< RepeatedRow<VectorRef,true> >,
				typename container_traits<VectorRef>::value_type >,
     public GenericVector< ConcatRows< RepeatedRow<VectorRef,true> >,
                           typename deref<VectorRef>::type::element_type > {
protected:
   typedef repeated_line_matrix<VectorRef> hidden_type;
   const hidden_type& hidden() const { return reinterpret_cast<const hidden_type&>(*this); }
public:
   const typename container_traits<VectorRef>::value_type& front() const
   {
      return hidden().get_vector().front();
   }

   int size() const
   {
      return hidden().get_count() * hidden().get_vector().dim();
   }
};

template <typename VectorRef, bool _same_elem>
class Rows< RepeatedRow<VectorRef,_same_elem> >
   : public constant_container< Rows< RepeatedRow<VectorRef,_same_elem> >,
				typename deref<VectorRef>::type > {
protected:
   ~Rows();
   typedef repeated_line_matrix<VectorRef> hidden_type;
   const hidden_type& hidden() const { return reinterpret_cast<const hidden_type&>(*this); }
public:
   const typename deref<VectorRef>::type& front() const { return hidden().get_vector(); }
   int size() const { return hidden().get_count(); }
};

template <typename VectorRef, bool _same_elem>
class Rows< RepeatedCol<VectorRef,_same_elem> >
   : public modified_container_impl< Rows< RepeatedCol<VectorRef,_same_elem> >,
               cons< Container< typename attrib<VectorRef>::plus_const >,
               cons< Accessor< operations::construct_unary_with_arg<SameElementVector,int> >,
		     Hidden< repeated_line_matrix<VectorRef> > > >
            > {
   typedef modified_container_impl<Rows> _super;
protected:
   ~Rows();
public:
   const typename _super::container& get_container() const { return this->hidden().get_vector(); }
   typename _super::accessor get_accessor() const { return this->hidden().get_count(); }
   int size() const { return get_container().dim(); }
};

template <typename VectorRef, bool _same_elem>
class Cols< RepeatedRow<VectorRef,_same_elem> >
   : public Rows< RepeatedCol<VectorRef,_same_elem> > { };

template <typename VectorRef, bool _same_elem>
class Cols< RepeatedCol<VectorRef,_same_elem> >
   : public Rows< RepeatedRow<VectorRef,_same_elem> > { };

template <typename VectorRef>
class ConcatRows< RepeatedCol<VectorRef, true> >
   : public ConcatRows< RepeatedRow<VectorRef, true> > { };

template <typename VectorRef, bool _same_elem>
class ConcatRows< Transposed< RepeatedRow<VectorRef,_same_elem> > >
   : public ConcatRows< RepeatedCol<VectorRef,_same_elem> > { };

template <typename VectorRef, bool _same_elem>
class ConcatRows< Transposed< RepeatedCol<VectorRef,_same_elem> > >
   : public ConcatRows< RepeatedRow<VectorRef,_same_elem> > { };

template <typename _Vector> inline
const RepeatedRow<typename prevent_lazy<const _Vector&>::type>
repeat_row(const GenericVector<_Vector>& v, int n)
{
   return RepeatedRow<typename prevent_lazy<const _Vector&>::type> (v.top(), n);
}

template <typename _Vector> inline
const RepeatedCol<typename prevent_lazy<const _Vector&>::type>
repeat_col(const GenericVector<_Vector>& v, int n)
{
   return RepeatedCol<typename prevent_lazy<const _Vector&>::type> (v.top(), n);
}

template <typename E> inline
const RepeatedRow< SameElementVector<const E&> >
same_element_matrix(const E& x, int m, int n)
{
   return RepeatedRow< SameElementVector<const E&> > (same_element_vector(x,n), m);
}

template <typename E> inline
const RepeatedRow< SameElementVector<E> >
ones_matrix(int m, int n)
{
   return RepeatedRow< SameElementVector<E> > (ones_vector<E>(n), m);
}

template <typename E> inline
const RepeatedRow< SameElementVector<E> >
zero_matrix(int m, int n)
{
   return RepeatedRow< SameElementVector<E> > (zero_vector<E>(n), m);
}

template <typename E>
class unlimited_constant_matrix : public unlimited_constant_container<E> {
   typedef unlimited_constant_container<E> super;
public:
   static const bool prefer_concat=true;
   unlimited_constant_matrix(typename function_argument<E>::type arg) : super(arg) { }

   typename super::reference operator() (int, int) const { return this->front(); }
};

template <typename E>
class unlimited_constant_matrix<E&> : public unlimited_constant_container<E&> {
   typedef unlimited_constant_container<E&> super;
protected:
   unlimited_constant_matrix();
   ~unlimited_constant_matrix();
public:
   static const bool prefer_concat=true;

   typename super::reference operator() (int, int) const { return this->front(); }
};

template <typename E>
class ConcatRows< unlimited_constant_matrix<E> >
   : public unlimited_constant_container<E> {
protected:
   ConcatRows();
   ~ConcatRows();
};

template <typename E>
class Rows< unlimited_constant_matrix<E> >
   : public unlimited_constant_container< const unlimited_constant_container<E>& > {
protected:
   Rows();
   ~Rows();
};

template <typename E>
class Cols< unlimited_constant_matrix<E> >
   : public unlimited_constant_container< const unlimited_constant_container<E>& > {
protected:
   Cols();
   ~Cols();
};

/* --------------------------------------
 *  diagonal and block-diagonal matrices
 * -------------------------------------- */

template <typename VectorRef>
class DiagMatrix
   : public GenericMatrix< DiagMatrix<VectorRef>, typename deref<VectorRef>::type::element_type > {
protected:
   alias<VectorRef> vector;
public:
   typedef typename deref<VectorRef>::type::value_type value_type;
   typedef typename deref<VectorRef>::type::const_reference const_reference;
   typedef const_reference reference;

   typedef typename alias<VectorRef>::arg_type arg_type;
   DiagMatrix(arg_type vector_arg) : vector(vector_arg) { }

   typename alias<VectorRef>::const_reference get_vector() const { return vector; }
};

template <typename VectorRef>
struct spec_object_traits< DiagMatrix<VectorRef> >
   : spec_object_traits<is_container> {
   static const bool is_temporary=true, is_always_const=true;
};

template <typename VectorRef>
class matrix_random_access_methods< DiagMatrix<VectorRef> > {
   typedef DiagMatrix<VectorRef> master;
public:
   typename container_traits<VectorRef>::const_reference
   operator() (int i, int j) const
   {
      const master& me=*static_cast<const master*>(this);
      INDEX_CHECK((i<0 || i>=me.rows()), "operator() - row index out of range");
      INDEX_CHECK((j<0 || j>=me.cols()), "operator() - column index out of range");
      if (i==j) return me.get_vector()[i];
      operations::clear<typename deref<VectorRef>::type::value_type> zero;
      return zero();
   }
};

template <typename VectorRef>
class ConcatRows< DiagMatrix<VectorRef> >
   : public modified_container_impl< ConcatRows< DiagMatrix<VectorRef> >,
               cons< Container< masquerade_add_features<VectorRef, cons<indexed, end_sensitive> > >,
               cons< Accessor< pair<nothing, operations::fix2<int, operations::composed12< BuildUnaryIt<operations::index2element>, void, BuildBinary<operations::mul> > > > >,
                     Hidden< DiagMatrix<VectorRef> > > >
            >,
     public GenericVector< ConcatRows< DiagMatrix<VectorRef> >,
			   typename deref<VectorRef>::type::element_type > {
   typedef modified_container_impl<ConcatRows> _super;
public:
   const typename _super::container& get_container() const
   {
      return reinterpret_cast<const typename _super::container&>(this->hidden().get_vector());
   }

   typename _super::accessor::second_type get_accessor() const
   {
      const int n=this->hidden().get_vector().dim();
      return n+1;
   }

   int dim() const
   {
      const int n=this->hidden().get_vector().dim();
      return n*n;
   }
};

template <typename VectorRef>
struct check_container_feature< DiagMatrix<VectorRef>, sparse> : True { };

template <typename VectorRef>
struct check_container_feature< DiagMatrix<VectorRef>, pure_sparse>
   : check_container_ref_feature<VectorRef, pure_sparse> { };

template <typename VectorRef>
class Rows< DiagMatrix<VectorRef> >
   : public modified_container_pair_impl< Rows< DiagMatrix<VectorRef> >,
               cons< Container1< sequence >,
               cons< Container2< masquerade_add_features<VectorRef, dense> >,
               cons< Accessor< UnitVector_factory<> >,
		     Hidden< DiagMatrix<VectorRef> > > > >
            > {
   typedef modified_container_pair_impl<Rows> _super;
protected:
   ~Rows();
public:
   sequence get_container1() const
   {
      return sequence(0,size());
   }
   const typename _super::container2& get_container2() const
   {
      return reinterpret_cast<const typename _super::container2&>(this->hidden().get_vector());
   }

   typename _super::accessor get_accessor() const { return size(); }
   int size() const { return this->hidden().get_vector().dim(); }
};

template <typename VectorRef>
class Cols< DiagMatrix<VectorRef> > : public Rows< DiagMatrix<VectorRef> > { };

template <typename Vector> inline
const DiagMatrix<const Vector&>
diag(const GenericVector<Vector>& v)
{
   return v.top();
}

template <typename E> inline
const DiagMatrix< SameElementVector<E> >
unit_matrix(int dim)
{
   return ones_vector<E>(dim);
}

template <typename MatrixRef1, typename MatrixRef2>
class BlockDiagMatrix
   : public container_pair_base<typename attrib<MatrixRef1>::plus_const, typename attrib<MatrixRef2>::plus_const>,
     public GenericMatrix< BlockDiagMatrix<MatrixRef1,MatrixRef2>,
			   typename deref<MatrixRef1>::type::element_type > {
   typedef container_pair_base<typename attrib<MatrixRef1>::plus_const, typename attrib<MatrixRef2>::plus_const> _base;
public:
   typedef typename identical<typename container_elem_traits<MatrixRef1>::value_type,
			      typename container_elem_traits<MatrixRef2>::value_type>::type
      value_type;
   typedef typename compatible<typename container_elem_traits<MatrixRef1>::const_reference,
			       typename container_elem_traits<MatrixRef2>::const_reference>::type
      const_reference;
   typedef const_reference reference;

   BlockDiagMatrix(typename _base::first_arg_type matrix1_arg, typename _base::second_arg_type matrix2_arg)
      : _base(matrix1_arg, matrix2_arg) { }

   int rows() const { return this->get_container1().rows() + this->get_container2().rows(); }
   int cols() const { return this->get_container1().cols() + this->get_container2().cols(); }
};

template <typename MatrixRef1, typename MatrixRef2>
struct spec_object_traits< BlockDiagMatrix<MatrixRef1, MatrixRef2> >
   : spec_object_traits<is_container> {
   static const bool is_temporary=true, is_always_const=true;
};

template <typename MatrixRef1, typename MatrixRef2>
struct check_container_feature<BlockDiagMatrix<MatrixRef1, MatrixRef2>, sparse> : True { };

template <typename MatrixRef1, typename MatrixRef2>
struct check_container_feature<BlockDiagMatrix<MatrixRef1, MatrixRef2>, pure_sparse> {
   static const bool answer=check_container_ref_feature<MatrixRef1, pure_sparse>::answer &&
                            check_container_ref_feature<MatrixRef2, pure_sparse>::answer;
};

template <typename MatrixRef1, typename MatrixRef2>
class matrix_random_access_methods< BlockDiagMatrix<MatrixRef1, MatrixRef2> > {
   typedef BlockDiagMatrix<MatrixRef1, MatrixRef2> master;
public:
   typename compatible<typename container_elem_traits<MatrixRef1>::const_reference,
		       typename container_elem_traits<MatrixRef2>::const_reference>::type
   operator() (int i, int j) const
   {
      master& me=static_cast<master&>(*this);
      INDEX_CHECK((i<0 || i>=me.rows() || j<0 || j>=me.cols()), "BlockDiagMatrix - indices out of range");
      const int r1=me.get_container1().rows(), c1=me.get_container1().cols();
      if (i<r1 && j<c1) return me.get_container1()(i,j);
      if (i>=r1 && j>=c1) return me.get_container2()(i-r1,j-c1);
      operations::clear<typename master::element_type> zero;
      return zero();
   }
};

template <typename MatrixRef1, typename MatrixRef2>
class block_diag_rows1 :
   public modified_container_impl< block_diag_rows1<MatrixRef1,MatrixRef2>,
	     cons< Container< masquerade<Rows, MatrixRef1> >,
	     cons< Accessor< ExpandedVector_factory<> >,
		   Hidden< BlockDiagMatrix<MatrixRef1,MatrixRef2> > > > > {
   typedef modified_container_impl<block_diag_rows1> _super;
public:
   const typename _super::container& get_container() const
   {
      return rows(this->hidden().get_container1());
   }
   typename _super::accessor get_accessor() const
   {
      return typename _super::accessor(0, this->hidden().cols());
   }
};

template <typename MatrixRef1, typename MatrixRef2>
class block_diag_rows2 :
   public modified_container_impl< block_diag_rows2<MatrixRef1,MatrixRef2>,
	     cons< Container< masquerade<Rows, MatrixRef2> >,
	     cons< Accessor< ExpandedVector_factory<> >,
		   Hidden< BlockDiagMatrix<MatrixRef1,MatrixRef2> > > > > {
   typedef modified_container_impl<block_diag_rows2> _super;
public:
   const typename _super::container& get_container() const
   {
      return rows(this->hidden().get_container2());
   }
   typename _super::accessor get_accessor() const
   {
      return typename _super::accessor(this->hidden().get_container1().cols(), this->hidden().cols());
   }
};

template <typename MatrixRef1, typename MatrixRef2>
class block_diag_cols1 :
   public modified_container_impl< block_diag_cols1<MatrixRef1,MatrixRef2>,
	     cons< Container< masquerade<Cols, MatrixRef1> >,
	     cons< Accessor< ExpandedVector_factory<> >,
		   Hidden< BlockDiagMatrix<MatrixRef1,MatrixRef2> > > > > {
   typedef modified_container_impl<block_diag_cols1> _super;
public:
   const typename _super::container& get_container() const
   {
      return cols(this->hidden().get_container1());
   }
   typename _super::accessor get_accessor() const
   {
      return typename _super::accessor(0, this->hidden().rows());
   }
};

template <typename MatrixRef1, typename MatrixRef2>
class block_diag_cols2 :
   public modified_container_impl< block_diag_cols2<MatrixRef1,MatrixRef2>,
	     cons< Container< masquerade<Cols, MatrixRef2> >,
	     cons< Accessor< ExpandedVector_factory<> >,
		   Hidden< BlockDiagMatrix<MatrixRef1,MatrixRef2> > > > > {
   typedef modified_container_impl<block_diag_cols2> _super;
public:
   const typename _super::container& get_container() const
   {
      return cols(this->hidden().get_container2());
   }
   typename _super::accessor get_accessor() const
   {
      return typename _super::accessor(this->hidden().get_container1().rows(), this->hidden().rows());
   }
};

template <typename MatrixRef1, typename MatrixRef2>
class Rows< BlockDiagMatrix<MatrixRef1,MatrixRef2> >
   : public container_chain_impl< Rows< BlockDiagMatrix<MatrixRef1,MatrixRef2> >,
	       cons< Container1< block_diag_rows1<MatrixRef1,MatrixRef2> >,
	       cons< Container2< block_diag_rows2<MatrixRef1,MatrixRef2> >,
		     Hidden< BlockDiagMatrix<MatrixRef1,MatrixRef2> > > > > {
   typedef container_chain_impl<Rows> _super;
protected:
   ~Rows();
};

template <typename MatrixRef1, typename MatrixRef2>
class Cols< BlockDiagMatrix<MatrixRef1,MatrixRef2> >
   : public container_chain_impl< Cols< BlockDiagMatrix<MatrixRef1,MatrixRef2> >,
	       cons< Container1< block_diag_cols1<MatrixRef1,MatrixRef2> >,
	       cons< Container2< block_diag_cols2<MatrixRef1,MatrixRef2> >,
		     Hidden< BlockDiagMatrix<MatrixRef1,MatrixRef2> > > > > {
   typedef container_chain_impl<Cols> _super;
protected:
   ~Cols();
};

template <typename E, typename Matrix1, typename Matrix2> inline
const BlockDiagMatrix<const Matrix1&, const Matrix2&>
diag(const GenericMatrix<Matrix1,E>& m1, const GenericMatrix<Matrix2,E>& m2)
{
   return BlockDiagMatrix<const Matrix1&, const Matrix2&>(m1.top(), m2.top());
}

template <typename E, typename Vector1, typename Matrix2> inline
const BlockDiagMatrix<DiagMatrix<const Vector1&>, const Matrix2&>
diag(const GenericVector<Vector1,E>& v1, const GenericMatrix<Matrix2,E>& m2)
{
   return BlockDiagMatrix<DiagMatrix<const Vector1&>, const Matrix2&>(diag(v1), m2.top());
}

template <typename E, typename Matrix1, typename Vector2> inline
const BlockDiagMatrix<const Matrix1&, DiagMatrix<const Vector2&> >
diag(const GenericMatrix<Matrix1,E>& m1, const GenericVector<Vector2,E>& v2)
{
   return BlockDiagMatrix<const Matrix1&, DiagMatrix<const Vector2&> >(m1.top(), diag(v2));
}

/* -----------------
 *  SparseMatrix2x2
 * ----------------- */

/* This is a matrix of arbitrary dimensions.
   Only 4 elements have distinct values: a_ii, a_ij, a_ji, a_jj.
   The rest elements are implicitly 1 if on the diagonal and 0 otherwise.
*/
template <typename E>
class SparseMatrix2x2 {
public:
   typedef E value_type;
   typedef E& reference;
   typedef const E& const_reference;

   int i, j;	// row & col
   E a_ii, a_ij, a_ji, a_jj;

   SparseMatrix2x2() { }

   SparseMatrix2x2(int i_arg, int j_arg,
		   typename function_argument<E>::type a_ii_arg, typename function_argument<E>::type a_ij_arg,
		   typename function_argument<E>::type a_ji_arg, typename function_argument<E>::type a_jj_arg)
      : i(i_arg), j(j_arg), a_ii(a_ii_arg), a_ij(a_ij_arg), a_ji(a_ji_arg), a_jj(a_jj_arg) { }
};

template <typename E>
class Transposed< SparseMatrix2x2<E> > : public SparseMatrix2x2<E> {
protected:
   ~Transposed();
};

template <typename E> inline
E det(const SparseMatrix2x2<E>& U)
{
   return U.a_ii*U.a_jj-U.a_ij*U.a_ji;
}

/* ------------------------------------
 *  MatrixProduct
 * lazy evaluation of a matrix product
 * ------------------------------------ */

template <typename MatrixRef1, typename MatrixRef2>
struct lazy_product_traits {
   typedef typename deref<typename operations::mul<typename deref<MatrixRef1>::type::row_type,
						   typename deref<MatrixRef2>::type::col_type>::result_type>::type
      reference;
   typedef typename deref<reference>::type value_type;
   typedef typename object_traits<value_type>::persistent_type element_type;
};

template <typename MatrixRef1, typename MatrixRef2>
class MatrixProduct
   : public container_pair_base<MatrixRef1, MatrixRef2>,
     public GenericMatrix< MatrixProduct<MatrixRef1,MatrixRef2>,
			   typename lazy_product_traits<MatrixRef1,MatrixRef2>::element_type > {
   typedef container_pair_base<MatrixRef1, MatrixRef2> _base;
public:
   typedef typename lazy_product_traits<MatrixRef1,MatrixRef2>::reference reference;
   typedef typename lazy_product_traits<MatrixRef1,MatrixRef2>::value_type value_type;
   typedef reference const_reference;
   typedef typename lazy_product_traits<MatrixRef1,MatrixRef2>::element_type element_type;

   MatrixProduct(typename _base::first_arg_type src1_arg, typename _base::second_arg_type src2_arg)
      : _base(src1_arg,src2_arg) { }

   static const bool is_sparse=check_container_ref_feature<MatrixRef1, sparse>::answer &&
                               check_container_ref_feature<MatrixRef2, sparse>::answer;
   typedef typename if_else<is_sparse, SparseMatrix<element_type>, Matrix<element_type> >::result
      persistent_type;
};

template <typename MatrixRef1, typename MatrixRef2>
struct spec_object_traits< MatrixProduct<MatrixRef1, MatrixRef2> >
   : spec_object_traits<is_container> {
   static const bool is_lazy=true, is_temporary=true, is_always_const=true;
};

template <typename MatrixRef1, typename MatrixRef2>
class matrix_random_access_methods< MatrixProduct<MatrixRef1, MatrixRef2> > {
   typedef MatrixProduct<MatrixRef1,MatrixRef2> master;
public:
   typename operations::mul<typename container_elem_traits<MatrixRef1>::const_reference,
			    typename container_elem_traits<MatrixRef2>::const_reference>::result_type
   operator() (int i, int j) const
   {
      const master& me=static_cast<const master&>(*this);
      INDEX_CHECK((i<0 || i>=me.get_container1().rows()), "operator() - row index out of range");
      INDEX_CHECK((j<0 || j>=me.get_container2().cols()), "operator() - column index out of range");
      return me.get_container1().row(i) * me.get_container2().col(j);
   }
};

template <typename MatrixRef1, typename MatrixRef2>
class ConcatRows< MatrixProduct<MatrixRef1,MatrixRef2> >
   : public container_product_impl< ConcatRows< MatrixProduct<MatrixRef1,MatrixRef2> >,
	       cons< Container1< masquerade<Rows,MatrixRef1> >,
	       cons< Container2< masquerade<Cols,MatrixRef2> >,
	       cons< Accessor< BuildBinary<operations::mul> >,
		     Hidden< MatrixProduct<MatrixRef1,MatrixRef2> > > > >
            >,
     public GenericVector< ConcatRows< MatrixProduct<MatrixRef1,MatrixRef2> >,
			   typename lazy_product_traits<MatrixRef1,MatrixRef2>::element_type > {
   typedef container_product_impl<ConcatRows> _super;
public:
   const typename _super::container1& get_container1() const
   {
      return rows(this->hidden().get_container1());
   }
   const typename _super::container2& get_container2() const
   {
      return cols(this->hidden().get_container2());
   }
};

template <typename MatrixRef1, typename MatrixRef2>
class Rows< MatrixProduct<MatrixRef1,MatrixRef2> >
   : public modified_container_pair_impl< Rows< MatrixProduct<MatrixRef1,MatrixRef2> >,
	       cons< Container1< masquerade<pm::Rows,MatrixRef1> >,
	       cons< Container2< unlimited_constant_container<typename attrib<MatrixRef2>::plus_const_ref> >,
	       cons< Accessor< BuildBinary<operations::mul> >,
		     Hidden< MatrixProduct<MatrixRef1,MatrixRef2> > > > >
            > {
   typedef modified_container_pair_impl<Rows> _super;
protected:
   ~Rows();
public:
   const typename _super::container1& get_container1() const
   {
      return rows(this->hidden().get_container1());
   }
   const typename _super::container2& get_container2() const
   {
      return constant(this->hidden().get_container2());
   }
};

template <typename MatrixRef1, typename MatrixRef2>
class Cols< MatrixProduct<MatrixRef1,MatrixRef2> >
   : public modified_container_pair_impl< Cols< MatrixProduct<MatrixRef1,MatrixRef2> >,
	       cons< Container1< unlimited_constant_container<typename attrib<MatrixRef1>::plus_const_ref> >,
	       cons< Container2< masquerade<pm::Cols,MatrixRef2> >,
	       cons< Accessor< BuildBinary<operations::mul> >,
		     Hidden< MatrixProduct<MatrixRef1,MatrixRef2> > > > >
            > {
   typedef modified_container_pair_impl<Cols> _super;
protected:
   ~Cols();
public:
   const typename _super::container1& get_container1() const
   {
      return constant(this->hidden().get_container1());
   }
   const typename _super::container2& get_container2() const
   {
      return cols(this->hidden().get_container2());
   }
};

/* ---------------------------------------
 *  matrix arithmetic operations
 * --------------------------------------- */
namespace operations {

template <typename OpRef>
struct neg_impl<OpRef, is_matrix> {
   typedef OpRef argument_type;
   typedef const LazyMatrix1<typename attrib<OpRef>::plus_const, BuildUnary<neg> > result_type;

   result_type operator() (typename function_argument<OpRef>::const_type x) const { return x; }

   typename function_argument<OpRef>::lvalue
   assign(typename function_argument<OpRef>::lvalue x) const
   {
      return static_cast<typename function_argument<OpRef>::lvalue> (inv_sign(x));
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct add_impl<LeftRef, RightRef, check, cons<is_matrix, is_matrix> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyMatrix2<typename attrib<LeftRef>::plus_const,
			     typename attrib<RightRef>::plus_const, BuildBinary<add> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.rows() != r.rows() ||
				  l.cols() != r.cols()), "operator+ - matrix dimension mismatch");
      return result_type(l,r);
   }

   template <typename Iterator2>
   typename function_argument<LeftRef>::const_type
   operator() (partial_left, typename function_argument<LeftRef>::const_type l, const Iterator2&) const
   {
      return l;
   }

   template <typename Iterator1>
   typename function_argument<RightRef>::const_type
   operator() (partial_right, const Iterator1&, typename function_argument<RightRef>::const_type r) const
   {
      return r;
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l+=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct sub_impl<LeftRef, RightRef, check, cons<is_matrix, is_matrix> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyMatrix2<typename attrib<LeftRef>::plus_const,
			     typename attrib<RightRef>::plus_const, BuildBinary<sub> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.rows() != r.rows() ||
				  l.cols() != r.cols()), "operator- - matrix dimension mismatch");
      return result_type(l,r);
   }

   template <typename Iterator2>
   typename function_argument<LeftRef>::const_type
   operator() (partial_left, typename function_argument<LeftRef>::const_type l, const Iterator2&) const
   {
      return l;
   }

   template <typename Iterator1>
   typename neg<RightRef>::result_type
   operator() (partial_right, const Iterator1&, typename function_argument<RightRef>::const_type r) const
   {
      neg<RightRef> n;
      return n(r);
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l-=r);
   }
};

template <typename LeftRef, typename RightRef, bool check,
	  bool left_diag=derived_from_instance<typename deref<LeftRef>::type, DiagMatrix>::answer,
	  bool right_diag=derived_from_instance<typename deref<RightRef>::type, DiagMatrix>::answer>
struct matrix_prod_chooser {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const MatrixProduct<typename prevent_lazy<LeftRef>::type,
			       typename prevent_lazy<RightRef>::type> result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.cols() != r.rows()), "operator* - matrix dimension mismatch");
      return result_type(l,r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct matrix_prod_chooser<LeftRef, RightRef, check, true, false> {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef RepeatedCol<typename extract_template_type_param<typename deref<LeftRef>::type,0>::type> Left;
   typedef const LazyMatrix2<Left, RightRef, BuildBinary<mul> > result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.cols() != r.rows()), "operator* - matrix dimension mismatch");
      return result_type(Left(l.get_vector(),r.cols()), r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct matrix_prod_chooser<LeftRef, RightRef, check, false, true> {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef RepeatedRow<typename extract_template_type_param<typename deref<RightRef>::type,0>::type> Right;
   typedef const LazyMatrix2<LeftRef, Right, BuildBinary<mul> > result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.cols() != r.rows()), "operator* - matrix dimension mismatch");
      return result_type(l, Right(r.get_vector(),l.rows()));
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct matrix_prod_chooser<LeftRef, RightRef, check, true, true> {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef LazyVector2<typename extract_template_type_param<typename deref<LeftRef>::type,0>::type,
		       typename extract_template_type_param<typename deref<RightRef>::type,0>::type,
		       BuildBinary<mul> >
      prod_diag;
   typedef const DiagMatrix<prod_diag> result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.cols() != r.rows()), "operator* - matrix dimension mismatch");
      return result_type(prod_diag(l.get_vector(), r.get_vector()));
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct mul_impl<LeftRef, RightRef, check, cons<is_matrix, is_matrix> >
   : matrix_prod_chooser<LeftRef,RightRef,check> { };

template <typename LeftRef, typename RightRef, bool check>
struct mul_impl<LeftRef, RightRef, check, cons<is_matrix, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyVector2<pm::masquerade<Rows, typename attrib<LeftRef>::plus_const>,
			     pm::masquerade_keep_ref<unlimited_constant_container, typename prevent_lazy<RightRef>::type>,
                             BuildBinary<mul> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.cols() != r.dim()), "operator* - matrix dimension mismatch");
      return result_type(l,r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct mul_impl<LeftRef, RightRef, check, cons<is_vector, is_matrix> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyVector2<pm::masquerade_keep_ref<unlimited_constant_container, typename prevent_lazy<LeftRef>::type>,
			     pm::masquerade<Cols, typename attrib<RightRef>::plus_const>,
                             BuildBinary<mul> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      if (check) DIMENSION_CHECK((l.dim() != r.rows()), "operator* - matrix dimension mismatch");
      return result_type(l,r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct mul_impl<LeftRef, RightRef, check, cons<is_matrix, is_scalar> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyMatrix2<typename attrib<LeftRef>::plus_const,
			     pm::masquerade_keep_ref<unlimited_constant_matrix, typename prevent_lazy<RightRef>::type>,
			     BuildBinary<mul> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      return result_type(l,r);
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l*=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct mul_impl<LeftRef, RightRef, check, cons<is_scalar, is_matrix> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyMatrix2<pm::masquerade_keep_ref<unlimited_constant_matrix, typename prevent_lazy<LeftRef>::type>,
			     typename attrib<RightRef>::plus_const,
			     BuildBinary<mul> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      return result_type(l,r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct div_impl<LeftRef, RightRef, check, cons<is_matrix, is_scalar> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef const LazyMatrix2<typename attrib<LeftRef>::plus_const, 
			     pm::masquerade_keep_ref<unlimited_constant_matrix, typename prevent_lazy<RightRef>::type>,
			     BuildBinary<div> >
      result_type;

   result_type operator() (typename function_argument<LeftRef>::const_type l,
			   typename function_argument<RightRef>::const_type r) const
   {
      return result_type(l,r);
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l/=r);
   }
};

/* ------------------------------------
 *  operations building block matrices
 * ------------------------------------ */

template <typename LeftRef, typename RightRef, bool check>
struct div_impl<LeftRef, RightRef, check, cons<is_matrix, is_matrix> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef RowChain<LeftRef, RightRef> chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;
   
   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      if (check) DIMENSION_CHECK((l.cols() && r.cols() && l.cols() != r.cols()),
				 "concat_rows - column dimensions mismatch");
      return chain(l,r);
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l/=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct div_impl<LeftRef, RightRef, check, cons<is_matrix, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef RowChain<LeftRef, pm::masquerade<SingleRow, RightRef> > chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;

   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      if (check) DIMENSION_CHECK((l.cols() && r.dim() && l.cols() != r.dim()),
				 "concat_rows - column dimensions mismatch");
      return chain(l,r);
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l/=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct div_impl<LeftRef, RightRef, check, cons<is_vector, is_matrix> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef RowChain<pm::masquerade<SingleRow, LeftRef>, RightRef> chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;

   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      if (check) DIMENSION_CHECK((l.dim() && r.cols() && l.dim() != r.cols()),
				 "concat_rows - column dimensions mismatch");
      return chain(l,r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct div_impl<LeftRef, RightRef, check, cons<is_vector, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef RowChain<pm::masquerade<SingleRow, LeftRef>, pm::masquerade<SingleRow, RightRef> > chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;

   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      if (check) DIMENSION_CHECK((l.dim() != r.dim()),
				 "concat_rows - column dimensions mismatch");
      return chain(l,r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct bitwise_or_impl<LeftRef, RightRef, check, cons<is_matrix, is_matrix> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef ColChain<LeftRef, RightRef> chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;

   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      if (check) DIMENSION_CHECK((l.rows() && r.rows() && l.rows() != r.rows()),
				 "concat_cols - row dimensions mismatch");
      return chain(l,r);
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l|=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct bitwise_or_impl<LeftRef, RightRef, check, cons<is_matrix, is_vector> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef ColChain<LeftRef, pm::masquerade<SingleCol, RightRef> > chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;

   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      if (check) DIMENSION_CHECK((l.rows() && r.dim() && l.rows() != r.dim()),
				 "concat_cols - row dimensions mismatch");
      return chain(l,r);
   }

   typename function_argument<LeftRef>::lvalue
   assign(typename function_argument<LeftRef>::lvalue l,
	  typename function_argument<RightRef>::const_type r) const
   {
      return static_cast<typename function_argument<LeftRef>::lvalue> (l|=r);
   }
};

template <typename LeftRef, typename RightRef, bool check>
struct bitwise_or_impl<LeftRef, RightRef, check, cons<is_vector, is_matrix> > {
   typedef LeftRef first_argument_type;
   typedef RightRef second_argument_type;
   typedef ColChain<pm::masquerade<SingleCol, LeftRef>, RightRef> chain;
   typedef typename assign_const<chain, object_traits<chain>::is_always_const>::type result_type;

   result_type operator() (typename function_argument<LeftRef>::type l,
			   typename function_argument<RightRef>::type r) const
   {
      if (check) DIMENSION_CHECK((l.dim() && r.rows() && l.dim() != r.rows()),
				 "concat_cols - row dimensions mismatch");
      return chain(l,r);
   }
};

} // end namespace operations
namespace operators {

template <typename Matrix1, typename Matrix2> inline
typename operations::add_impl<const Matrix1&, const Matrix2&, true>::result_type
operator+ (const GenericMatrix<Matrix1>& l, const GenericMatrix<Matrix2>& r)
{
   operations::add_impl<const Matrix1&, const Matrix2&, true> op;
   return op(l.top(), r.top());
}

template <typename Matrix1, typename Matrix2> inline
typename operations::sub_impl<const Matrix1&, const Matrix2&, true>::result_type
operator- (const GenericMatrix<Matrix1>& l, const GenericMatrix<Matrix2>& r)
{
   operations::sub_impl<const Matrix1&, const Matrix2&, true> op;
   return op(l.top(), r.top());
}

template <typename Matrix1, typename Matrix2> inline
bool operator== (const GenericMatrix<Matrix1>& l, const GenericMatrix<Matrix2>& r)
{
   DIMENSION_CHECK((l.rows() != r.rows() ||
		    l.cols() != r.cols()),  "operator== - matrix dimension mismatch");
   operations::eq<const Matrix1&, const Matrix2&> op;
   return op(l.top(), r.top());
}

template <typename Matrix1, typename Matrix2> inline
bool operator< (const GenericMatrix<Matrix1>& l, const GenericMatrix<Matrix2>& r)
{
   DIMENSION_CHECK((l.rows() != r.rows() ||
		    l.cols() != r.cols()),  "operator< - matrix dimension mismatch");
   operations::lt<const Matrix1&, const Matrix2&> op;
   return op(l.top(), r.top());
}

} // end namespace operators
} // end namespace pm

namespace polymake {
   using pm::GenericMatrix;
   using pm::ConcatRows;
   using pm::same_element_matrix;
   using pm::convert_to;
   using pm::zero_matrix;
   using pm::ones_matrix;
   using pm::unit_matrix;
   using pm::SparseMatrix2x2;
}

namespace std {
   template <typename Matrix1, typename Matrix2, typename E> inline
   void swap(pm::GenericMatrix<Matrix1,E>& M1, pm::GenericMatrix<Matrix2,E>& M2)
   {
      M1.top().swap(M2.top());
   }

   template <typename Matrix1, typename Matrix2, typename E> inline
   void swap(const pm::GenericMatrix<Matrix1,E>& M1, const pm::GenericMatrix<Matrix2,E>& M2)
   {
      pm::non_const(M1.top()).swap(pm::non_const(M2.top()));
   }

   // due to silly overloading rules
   template <typename Matrix>
   void swap(pm::ConcatRows<Matrix>& M1, pm::ConcatRows<Matrix>& M2)
   {
      M1.swap(M2);
   }
}

#ifdef _POLYMAKE_GENERIC_INCIDENCE_MATRIX_H
#  include <Incidence_and_SparseMatrix.h>
#endif
#endif // _POLYMAKE_GENERIC_MATRIX_H

// Local Variables:
// mode:C++
// End:
