/Users/mourrain/Devel/mmx/algebramix/include/algebramix/matrix.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix.hpp
* DESCRIPTION: Dense matrices
* COPYRIGHT  : (C) 2004--2007  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_MATRIX_HPP
#define __MMX_MATRIX_HPP
#include <basix/iterator.hpp>
#include <basix/vector.hpp>
#include <algebramix/permutation.hpp>
#include <algebramix/matrix_naive.hpp>

namespace mmx {
#define Matrix_variant(C) matrix_variant_helper<C>::MV
#define TMPL_DEF template<typename C, typename V= typename Matrix_variant(C)>
#define TMPL template<typename C, typename V>
#define TMPLK template<typename C, typename V, typename K>
#define Format format<C>
#define Matrix matrix<C,V>
#define Matrix_rep matrix_rep<C,V>
#define Abs_matrix matrix<Abs_type(C),V>
#define Real_matrix matrix<Real_type(C),V>
#define Center_matrix matrix<Center_type(C),V>
#define Radius_matrix matrix<Radius_type(C),V>
TMPL class matrix_rep;
TMPL class matrix;
TMPL inline nat cols (const Matrix& m);
TMPL inline nat rows (const Matrix& m);
TMPL inline bool is_a_scalar   (const Matrix& m);
TMPL inline bool is_non_scalar (const Matrix& m);
TMPL inline C* tab (Matrix& m);
TMPL inline const C* tab (const Matrix& m);

/******************************************************************************
* Dense matrices of arbitrary sizes
******************************************************************************/

TMPL_DEF
class matrix_rep REP_STRUCT_1(C) {
  C*   a;
  nat  nr;
  nat  nc;
  bool scalar_flag;
public:
  inline matrix_rep (C* a2, nat nr2, nat nc2, bool flag, const Format& fm):
    Format (fm), a (a2), nr (nr2), nc (nc2), scalar_flag (flag) {}
  inline ~matrix_rep () { mmx_delete<C> (a, aligned_size<C,V> (nr * nc)); }
  friend class Matrix;
  friend nat cols LESSGTR (const Matrix& m);
  friend nat rows LESSGTR (const Matrix& m);
  friend C* tab LESSGTR (Matrix& m);
  friend const C* tab LESSGTR (const Matrix& m);
  friend bool is_a_scalar   LESSGTR (const Matrix& m);
  friend bool is_non_scalar LESSGTR (const Matrix& m);
};

TMPL_DEF
class matrix {
INDIRECT_PROTO_2 (matrix, matrix_rep, C, V)
public:
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<matrix_linear,V> Mat;
private:
  static C zero;
public:
  matrix (const Format& fm= Format (no_format ())) {
    nat nc = Mat::def_rows;
    nat nr = Mat::def_cols;
    nat l  = aligned_size<C,V> (nr * nc);
    C*  tab= mmx_formatted_new<C> (l, fm);
    rep= new Matrix_rep (tab, nr, nc, false, fm);
  }

  template<typename T>
  matrix (const T& c) {
    Format fm= get_format (as<C> (c));
    nat nr = Mat::init_rows;
    nat nc = Mat::init_cols;
    nat l  = aligned_size<C,V> (nr * nc);
    C*  tab= mmx_formatted_new<C> (l, fm);
    rep= new Matrix_rep (tab, nr, nc, true, fm);
    Mat::set (rep->a, C(c), nr, nc);
  }

  template<typename T>
  matrix (const T& c, const Format& fm) {
    nat nr = Mat::init_rows;
    nat nc = Mat::init_cols;
    nat l  = aligned_size<C,V> (nr * nc);
    C*  tab= mmx_formatted_new<C> (l, fm);
    rep= new Matrix_rep (tab, nr, nc, true, fm);
    Mat::set (rep->a, C(c), nr, nc);
  }

  template<typename T, typename W>
  matrix (const matrix<T,W>& m) {
    Format fm= as<Format > (CF(m));
    nat l= aligned_size<C,V> (rows(m) * cols(m));
    C*  t= mmx_formatted_new<C> (l, fm);
    rep= new Matrix_rep (t, rows(m), cols(m), is_a_scalar (m), fm);
    Vec::cast (rep->a, tab(m), rows(m) * cols(m));
  }

  template<typename T, typename W>
  matrix (const matrix<T,W>& m, const Format& fm) {
    nat l= aligned_size<C,V> (rows(m) * cols(m));
    C*  t= mmx_formatted_new<C> (l, fm);
    rep= new Matrix_rep (t, rows(m), cols(m), is_a_scalar (m), fm);
    Vec::cast (rep->a, tab(m), rows(m) * cols(m));
  }

  matrix (const C& x) {
    Format fm= get_format (x);
    nat nr = Mat::init_rows;
    nat nc = Mat::init_cols;
    nat l  = aligned_size<C,V> (nr * nc);
    C*  tab= mmx_formatted_new<C> (l, fm);
    rep= new Matrix_rep (tab, nr, nc, true, fm);
    Mat::set (rep->a, x, nr, nc);
  }

  matrix (const C& x, nat nr, nat nc) {
    Format fm= get_format (x);
    nat l  = aligned_size<C,V> (nr * nc);
    C*  tab= mmx_formatted_new<C> (l, fm);
    rep= new Matrix_rep (tab, nr, nc, false, fm);
    Mat::set (rep->a, x, nr, nc);
  }

  matrix (C* tab, nat nr, nat nc, const Format& fm) {
    rep= new Matrix_rep (tab, nr, nc, false, fm);
  }

  matrix (C* tab, nat nr, nat nc, bool flag, const Format& fm) {
    rep= new Matrix_rep (tab, nr, nc, flag, fm);
  }

  matrix (const iterator<C>& it) {
    Format fm= CF(it);
    nat i, j, nr= Mat::def_rows, nc= Mat::def_cols;
    if (nr==0) nr=1;
    nat l  = aligned_size<C,V> (nr * nc);
    C*  tab= mmx_formatted_new<C> (l, fm);
    for (i=0; i<nr; i++, ++it)
      for (j=0; j<nc; j++, ++it)
        tab[Mat::index (i, j, nr, nc)]= *it;
    rep= new Matrix_rep (tab, nc, nr, false, fm);
  }

  inline const C& operator () (nat i, nat j) const {
    VERIFY (is_non_scalar (*this), "non-scalar matrix expected");
    VERIFY (i<rep->nr, "row out of range");
    VERIFY (j<rep->nc, "column out of range");
    return rep->a[Mat::index (i, j, rows (*this), cols(*this))]; }
  inline C& operator () (nat i, nat j) {
    VERIFY (is_non_scalar (*this), "non-scalar matrix expected");
    VERIFY (i<rep->nr, "row out of range");
    VERIFY (j<rep->nc, "column out of range");
    return rep->a[Mat::index (i, j, rows (*this), cols(*this))]; }
  inline C scalar () const {
    VERIFY (is_a_scalar (*this), "scalar matrix expected");
    return rep->a[0]; }
  inline C& scalar () {
    VERIFY (is_a_scalar (*this), "scalar matrix expected");
    return rep->a[0]; }
};
INDIRECT_IMPL_2 (matrix, matrix_rep, typename C, C, typename V, V)
DEFINE_UNARY_FORMAT_2 (matrix)

TMPL C Matrix::zero= C(0);
TMPL inline nat rows (const Matrix& m) { return m->nr; }
TMPL inline nat cols (const Matrix& m) { return m->nc; }
TMPL inline nat nbcol (const Matrix& m) {return cols(m);}
TMPL inline nat nbrow (const Matrix& m) {return rows(m);} 
TMPL inline nat N (const Matrix& m) { return rows (m) * cols (m); } 
TMPL inline Format CF (const Matrix& m) { return m->tfm (); }
TMPL inline bool is_a_scalar (const Matrix& m) { return m->scalar_flag; }
TMPL inline bool is_non_scalar (const Matrix& m) { return !m->scalar_flag; }
TMPL inline bool is_square_matrix (const Matrix& m) {
  return rows (m) == cols (m); }
TMPL inline const C& read (const Matrix& m, nat i, nat j) { return m(i,j); }
TMPL inline const C* tab (const Matrix& m) { return m->a; }
TMPL inline C* tab (Matrix& m) {
  VERIFY (is_non_scalar (m), "non-scalar matrix expected");
  m.secure(); return m->a; }
template<typename C, typename V, typename C2, typename V2> inline Matrix
extend (const Matrix& m, const matrix<C2,V2>& n) {
  VERIFY (is_a_scalar (m), "scalar matrix expected");
  VERIFY (is_non_scalar (n), "non-scalar matrix expected");
  return Matrix (m.scalar(), rows (n), cols (n)); }

STYPE_TO_TYPE(TMPL,scalar_type,Matrix,C);
UNARY_RETURN_TYPE(TMPL,abs_op,Matrix,Abs_matrix);
UNARY_RETURN_TYPE(TMPL,Re_op,Matrix,Real_matrix);
UNARY_RETURN_TYPE(TMPL,center_op,Matrix,Center_matrix);
UNARY_RETURN_TYPE(TMPL,radius_op,Matrix,Radius_matrix);

TMPL
struct fast_helper<Matrix > {
  typedef Fast_type(C) FC;
  typedef matrix<FC> fast_type;
  typedef typename Matrix_variant(FC) FV;
  static inline fast_type dd (const Matrix& m) {
    typedef implementation<vector_linear,FV> FVec;
    format<FC> fm= fast (CF(m));
    nat l= aligned_size<FC,FV> (rows (m) * cols (m));
    FC* r= mmx_formatted_new<FC> (l, fm);
    FVec::template vec_unary<fast_op> (r, tab (m), rows (m) * cols (m));
    return fast_type (r, rows (m), cols (m), is_a_scalar (m), fm); }
  static inline Matrix uu (const fast_type& m) {
    typedef implementation<vector_linear,V> Vec;
    Format fm= slow<C> (CF(m));
    nat l= aligned_size<C,V> (rows (m) * cols (m));
    C*  r= mmx_formatted_new<C> (l, fm);
    Vec::template vec_unary<slow_op> (r, tab (m), rows (m) * cols (m));
    return Matrix (r, rows (m), cols (m), is_a_scalar (m), fm); }
};

template<typename T, typename F, typename TV, typename FV>
struct as_helper<matrix<T,TV>,matrix<F,FV> > {
  static inline matrix<T,TV>
  cv (const matrix<F,FV>& m) {
    format<T> fm= as<format<T> > (CF(m));
    nat n= rows(m) * cols(m);
    nat l= aligned_size<T,TV> (rows (m) * cols (m));
    const F* sm= tab (m);
    T* c= mmx_formatted_new<T> (l, fm);
    for (nat i=0; i<n; i++) c[i]= as<T> (sm[i]);
    return matrix<T,TV> (c, rows (m), cols (m), is_a_scalar (m), fm); }
};

template<typename T, typename F, typename TV, typename FV> inline void
set_as (matrix<T,TV>& r, const matrix<F,FV>& m) {
  r= matrix<T,TV> (m, CF(r));
}

template<typename C, typename V, typename T> inline void
set_as (Matrix& r, const T& x) {
  r= Matrix (x, CF(r));
}

/******************************************************************************
* Fixed size dense matrices
******************************************************************************/

template<typename V, typename RS, typename CS>
struct matrix_fixed: public V {
  typedef implementation<vector_linear,V> Vec;
  typedef matrix_fixed<typename V::Naive,RS,CS> Naive;
  typedef matrix_fixed<typename V::Positive,RS,CS> Positive;
  typedef matrix_fixed<typename V::No_simd,RS,CS> No_simd;
  typedef matrix_fixed<typename V::No_thread,RS,CS> No_thread;
  typedef matrix_fixed<typename V::No_scaled,RS,CS> No_scaled;
};

template<typename F, typename V, typename RS, typename CS>
struct implementation<F,matrix_fixed<V,RS,CS> >:
  public implementation<F,V> {};

template<typename V, typename RS, typename CS>
struct implementation<matrix_defaults,matrix_fixed<V,RS,CS> > {
  static const nat def_rows = RS::val;
  static const nat def_cols = CS::val;
  static const nat init_rows= RS::val;
  static const nat init_cols= CS::val;
};

#define FMatrix matrix<C,matrix_fixed<V,RS,CS> >
#define FMatrix_rep matrix_rep<C,matrix_fixed<V,RS,CS> >

template<typename C, typename V, typename RS, typename CS>
class FMatrix_rep REP_STRUCT_1(C) {
  C* a;
public:
  inline matrix_rep (C* a2, nat, nat, bool, const Format& fm):
    Format (fm), a (a2) {}
  inline virtual ~matrix_rep () { mmx_delete<C> (a, RS::val * CS::val); }
  friend class FMatrix;
  friend nat cols LESSGTR (const FMatrix& m);
  friend nat rows LESSGTR (const FMatrix& m);
  friend C*  tab  LESSGTR (const FMatrix& m);
  friend bool is_a_scalar   LESSGTR (const Matrix& m);
  friend bool is_non_scalar LESSGTR (const Matrix& m);
};

template<typename C, typename V, typename RS, typename CS> inline nat
cols (const FMatrix& m) { return CS::val; }
template<typename C, typename V, typename RS, typename CS> inline nat
rows (const FMatrix& m) { return RS::val; }
template<typename C, typename V, typename RS, typename CS> inline bool
is_a_scalar (const FMatrix& m) { (void) m; return false; }
template<typename C, typename V, typename RS, typename CS> inline bool
is_non_scalar (const FMatrix& m) { (void) m; return true; }

#undef FMatrix
#undef FMatrix_rep

/******************************************************************************
* Matrix iteration
******************************************************************************/

TMPL_DEF
class matrix_iterator_rep: public iterator_rep<C> {
public:
  typedef implementation<matrix_linear,V> Mat;
private:
  const Matrix m;  // the matrix to traverse
  nat i;           // current row
  nat j;           // current column
  nat nc;
protected:
  bool is_busy () { return i < rows(m); }
  void advance () { j++; if (j >= cols(m)) { j=0; i++; } }
  C current () {return m (i, j); }
public:
  matrix_iterator_rep (const Matrix& m2, nat c = Mat::def_cols):
    iterator_rep<C> (CF(m)), m(m2), i(0), j(0), nc(c) {}
};

TMPL iterator<C>
iterate (const Matrix& m) {
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  return iterator<C> (new matrix_iterator_rep<C,V> (m));
}

/******************************************************************************
* Hashing and generic output routines for matrices
******************************************************************************/

TMPL syntactic
flatten (const Matrix& m) {
  if (is_a_scalar (m)) return flatten (m.scalar());
  int i, j, nr= rows(m), nc= cols(m);
  vector<syntactic> v;
  for (i=0; i<nr; i++) {
    vector<syntactic> h;
    for (j=0; j<nc; j++)
      h << flatten (m (i, j));
    v << apply (GEN_ROW, h);
  }
  return apply (GEN_SQTUPLE, v);
}

TMPL
struct binary_helper<Matrix >: public void_binary_helper<Matrix > {
  static inline string short_type_name () {
    return "M" * Short_type_name (C); }
  static inline generic full_type_name () {
    return gen ("Matrix", Full_type_name (C)); }
  static generic disassemble (const Matrix& m) {
    if (is_a_scalar (m)) return as<generic> (m.scalar ());
    vector<generic> v= vec (as<generic> (rows (m)), as<generic> (cols (m)));
    for (nat i=0; i<rows(m)*cols(m); i++)
      v << as<generic> (tab(m)[i]);
    return as<generic> (v); }
  static Matrix assemble (const generic& x) {
    // FIXME: what to do with the format?
    if (!is<vector<generic> > (x)) return Matrix (as<C> (x));
    vector<generic> v= as<vector<generic> > (x);
    nat rs= as<nat> (v[0]);
    nat cs= as<nat> (v[1]);
    nat l = aligned_size<C,V> (rs*cs);
    C*  r = mmx_new<C> (l);
    for (nat i=0; i<rs*cs; i++)
      r[i]= as<C> (v[i+2]);
    return Matrix (r, rs, cs, Format ()); }
  static void write (const port& out, const Matrix& m) {
    binary_write<Format > (out, CF (m));
    binary_write<bool> (out, is_a_scalar (m));
    if (is_a_scalar (m)) binary_write<C> (out, m.scalar ());
    else {
      binary_write<nat> (out, rows(m));
      binary_write<nat> (out, cols(m));
      for (nat i=0; i<rows(m)*cols(m); i++)
        binary_write<C> (out, tab(m)[i]); } }
  static Matrix read (const port& in) {
    Format fm= binary_read<Format > (in);
    bool sc= binary_read<bool> (in);
    if (sc) return Matrix (binary_read<C> (in), fm);
    else {
      nat rs= binary_read<nat> (in);
      nat cs= binary_read<nat> (in);
      nat l = aligned_size<C,V> (rs*cs);
      C*  r = mmx_formatted_new<C> (l, fm);
      for (nat i=0; i<rs*cs; i++)
        r[i]= binary_read<C> (in);
      return Matrix (r, rs, cs, fm); } }
};

/******************************************************************************
* Abstract routines on dense matrices
******************************************************************************/

template<typename Op, typename C, typename V> nat
unary_hash (const matrix<C,V>& m) {
  register nat i, j, h= 214365, nr= rows (m), nc= cols (m);
  for (i=0; i<nr; i++)
    for (j=0; j<nc; j++)
      h= (h<<1) ^ (h<<5) ^ (h>>27) ^ Op::op (m (i, j));
  return h;
}

template<typename Op, typename C, typename V>
matrix<Unary_return_type(Op,C),V>
unary_map (const matrix<C,V>& m) {
  typedef implementation<vector_linear,V> Vec;
  typedef Unary_return_type(Op,C) T;
  format<T> fm= unary_map<Op> (CF(m));
  if (is_a_scalar (m)) return matrix<T,V> (Op::op (m.scalar()));
  nat nrows= rows (m);
  nat ncols= cols (m);
  nat l= aligned_size<T,V> (nrows * ncols);
  T*  r= mmx_formatted_new<T> (l, fm);
  Vec::template vec_unary<Op> (r, tab (m), nrows*ncols);
  return matrix<T,V> (r, nrows, ncols, fm);
}

template<typename Op, typename C1, typename C2, typename V>
matrix<Binary_return_type(Op,C1,C2),V>
binary_map (const matrix<C1,V>& m, const matrix<C2,V>& n) {
  typedef implementation<vector_linear,V> Vec;
  typedef Binary_return_type(Op,C1,C2) T;
  format<T> fm= binary_map<Op> (CF(m), CF(n));
  if (is_a_scalar (m) || is_a_scalar (n)) {
    if (is_non_scalar (m)) return binary_map<Op> (m, extend (n, m));
    if (is_non_scalar (n)) return binary_map<Op> (extend (m, n), n);
    return matrix<T,V> (Op::op (m.scalar(), n.scalar()));
  }
  nat nrows= rows (m);
  nat ncols= cols (m);
  ASSERT (rows (n) == nrows, "unequal number of rows");
  ASSERT (cols (n) == ncols, "unequal number of columns");
  nat l= aligned_size<T,V> (nrows * ncols);
  T*  r= mmx_formatted_new<T> (l, fm);
  Vec::template vec_binary<Op> (r, tab (m), tab (n), nrows*ncols);
  return matrix<T,V> (r, nrows, ncols, fm);
}

template<typename Op, typename C, typename V, typename X>
matrix<Binary_return_type(Op,C,X),V>
binary_map_scalar (const matrix<C,V>& m, const X& x) {
  typedef implementation<vector_linear,V> Vec;
  typedef Binary_return_type(Op,C,X) T;
  format<T> fm= binary_map_scalar<C> (CF(m), x);
  if (is_a_scalar (m)) return matrix<T,V> (Op::op (m.scalar(), x));
  nat nrows= rows (m);
  nat ncols= cols (m);
  nat l= aligned_size<T,V> (nrows * ncols);
  T*  r= mmx_formatted_new<T> (l, fm);
  Vec::template vec_binary_scalar<Op> (r, tab (m), x, nrows*ncols);
  return matrix<T,V> (r, nrows, ncols, fm);
}

template<typename Op, typename C, typename V> matrix<C,V>&
nullary_set (matrix<C,V>& m) {
  typedef implementation<vector_linear,V> Vec;
  Vec::template vec_nullary<Op> (tab (m), rows (m) * cols (m));
  return m;
}

template<typename Op, typename T, typename C, typename V> matrix<T,V>&
unary_set (matrix<T,V>& m, const matrix<C,V>& n) {
  typedef implementation<vector_linear,V> Vec;
  if (is_a_scalar (m) || is_a_scalar (n)) {
    if (is_non_scalar (m))
      return unary_set<Op> (m, extend (n, m));
    else if (is_non_scalar (n))
      m= extend (m, n);
    else {
      Op::set_op (m.scalar(), n.scalar());
      return m;
    }
  }
  nat nrows= rows (m);
  nat ncols= cols (m);
  ASSERT (rows (n) == nrows, "unequal number of rows");
  ASSERT (cols (n) == ncols, "unequal number of columns");
  Vec::template vec_unary<Op> (tab (m), tab (n), nrows*ncols);
  return m;
}

template<typename Op, typename T, typename V, typename X> matrix<T,V>&
unary_set_scalar (matrix<T,V>& m, const X& x) {
  typedef implementation<vector_linear,V> Vec;
  if (is_a_scalar (m)) {
    Op::set_op (m.scalar(), x);
    return m;
  }
  nat nrows= rows (m);
  nat ncols= cols (m);
  Vec::template vec_unary_scalar<Op> (tab (m), x, nrows*ncols);
  return m;
}

template<typename Op, typename C1, typename C2, typename V> bool
binary_test (const matrix<C1,V>& m, const matrix<C2,V>& n) {
  typedef implementation<vector_linear,V> Vec;
  if (is_a_scalar (m) || is_a_scalar (n)) {
    if (is_non_scalar (m)) return binary_test<Op> (m, extend (n, m));
    if (is_non_scalar (n)) return binary_test<Op> (extend (m, n), n);
    return Op::op (m.scalar(), n.scalar());
  }
  nat nrows= rows (m);
  nat ncols= cols (m);
  if (rows (n) != nrows || cols (n) != ncols) return false;
  return Vec::template vec_binary_test<Op> (tab (m), tab (n), nrows*ncols);
}

template<typename Op, typename C, typename V> Unary_return_type(Op,C)
big (const matrix<C,V>& m) {
  typedef implementation<vector_linear,V> Vec;
  ASSERT (is_non_scalar (m), "non scalar matrix expected");
  return Vec::template vec_unary_big<Op> (tab (m), rows (m) * cols (m));
}

/******************************************************************************
* Dynamic mappers
******************************************************************************/

template<typename S1, typename D> matrix<D>
map (const function_1<D,Argument(S1) >& fun,
     const matrix<S1>& m, const format<D>& fm)
{
  nat n= rows(m) * cols(m);
  nat l= default_aligned_size<D> (n);
  D* r= mmx_formatted_new<D> (l, fm);
  const S1* src= tab (m);
  for (nat i=0; i<n; i++) r[i]= fun (src[i]);
  return matrix<D> (r, rows(m), cols(m), fm);
}

template<typename S1, typename D> inline matrix<D>
map (D (*fun) (const S1&),
     const matrix<S1>& x1, const format<D>& fm) {
  return map (function_1<D,Argument(S1) > (fun), x1, fm);
}

/******************************************************************************
* Constants
******************************************************************************/

TMPL void set_pi (Matrix& m) { nullary_set<pi_as_op> (m); }
TMPL void set_log2 (Matrix& m) { nullary_set<log2_as_op> (m); }
TMPL void set_euler (Matrix& m) { nullary_set<euler_as_op> (m); }
TMPL void set_catalan (Matrix& m) { nullary_set<catalan_as_op> (m); }
TMPL void set_imaginary (Matrix& m) { nullary_set<imaginary_as_op> (m); }
TMPL void set_nan (Matrix& m) { nullary_set<nan_as_op> (m); }
TMPL void set_fuzz (Matrix& m) { nullary_set<fuzz_as_op> (m); }
TMPL void set_smallest (Matrix& m) { nullary_set<smallest_as_op> (m); }
TMPL void set_largest (Matrix& m) { nullary_set<largest_as_op> (m); }
TMPL void set_accuracy (Matrix& m) { nullary_set<accuracy_as_op> (m); }
TMPL void set_infinity (Matrix& m) { nullary_set<infinity_as_op> (m); }
TMPL void set_maximal (Matrix& m) { nullary_set<maximal_as_op> (m); }
TMPL void set_minimal (Matrix& m) { nullary_set<minimal_as_op> (m); }

/******************************************************************************
* Instantiations of abstract routines on matrices
******************************************************************************/

TMPL Matrix copy (const Matrix& m) {
  return unary_map<id_op> (m); }
TMPL Matrix duplicate (const Matrix& m) {
  return unary_map<duplicate_op> (m); }

TMPL Matrix operator - (const Matrix& m) {
  return unary_map<neg_op> (m); }
TMPL Matrix operator + (const Matrix& m, const Matrix& n) {
  return binary_map<add_op> (m, n); }
TMPL Matrix operator + (const Matrix& m, const C& n) {
  return binary_map<add_op> (m, Matrix (n)); }
TMPL Matrix operator + (const C& m, const Matrix& n) {
  return binary_map<add_op> (Matrix (m), n); }
TMPL Matrix operator - (const Matrix& m, const Matrix& n) {
  return binary_map<sub_op> (m, n); }
TMPL Matrix operator - (const Matrix& m, const C& n) {
  return binary_map<sub_op> (m, Matrix (n)); }
TMPL Matrix operator - (const C& m, const Matrix& n) {
  return binary_map<sub_op> (Matrix (m), n); }
TMPL Matrix operator * (const Matrix& m, const C& c) {
  return binary_map_scalar<rmul_op> (m, c); }
TMPL Matrix operator * (const C& c, const Matrix& m) {
  return binary_map_scalar<lmul_op> (m, c); }
TMPL Matrix operator / (const Matrix& m, const C& c) {
  return binary_map_scalar<rdiv_op> (m, c); }
TMPL Matrix quo (const Matrix& m, const C& c) {
  return binary_map_scalar<rquo_op> (m, c); }
TMPL Matrix rem (const Matrix& m, const C& c) {
  return binary_map_scalar<rrem_op> (m, c); }
TMPL Matrix& operator += (Matrix& m, const Matrix& n) {
  return unary_set<add_op> (m, n); }
TMPL Matrix& operator -= (Matrix& m, const Matrix& n) {
  return unary_set<sub_op> (m, n); }
TMPL Matrix& operator *= (Matrix& m, const C& x) {
  return unary_set_scalar<rmul_op> (m, x); }
TMPL Matrix& operator /= (Matrix& m, const C& x) {
  return unary_set_scalar<rdiv_op> (m, x); }
TMPL bool operator <= (const Matrix& m, const Matrix& n) {
  return binary_test<lesseq_op> (m, n); }
TMPL bool operator >= (const Matrix& m, const Matrix& n) {
  return binary_test<gtreq_op> (m, n); }
template<typename C, typename V, typename K> bool
operator == (const Matrix& m, const K& c) {
  return m == Matrix (c, rows (m), cols (m)); }
template<typename C, typename V, typename K> bool
operator != (const Matrix& m, const K& c) {
  return m != Matrix (c, rows (m), cols (m)); }
template<typename C, typename V, typename K> bool
operator <= (const Matrix& m, const C& c) {
  return m <= Matrix (c, rows (m), cols (m)); }
template<typename C, typename V, typename K> bool
operator >= (const Matrix& m, const C& c) {
  return m >= Matrix (c, rows (m), cols (m)); }

TRUE_IDENTITY_OP_SUGAR(TMPL,Matrix)
EXACT_IDENTITY_OP_SUGAR(TMPL,Matrix)
ARITH_SCALAR_INT_SUGAR(TMPL,Matrix);
STRICT_COMPARE_SUGAR(TMPL,Matrix)

TMPL Matrix derive (const Matrix& m) {
  return unary_map<derive_op> (m); }
TMPL Matrix integrate (const Matrix& m) {
  return unary_map<integrate_op> (m); }

TMPL inline C big_add (const Matrix& m) { return big<add_op> (m); }
TMPL inline C big_sup (const Matrix& m) { return big<sup_op> (m); }
TMPL inline C big_max (const Matrix& m) { return big<max_op> (m); }

/******************************************************************************
* Floating point related functions
******************************************************************************/

TMPL inline bool is_finite (const Matrix& m) {
  if (is_a_scalar (m)) return is_finite (m.scalar());
  return big<and_is_finite_op> (m); }
TMPL inline bool is_nan (const Matrix& m) {
  if (is_a_scalar (m)) return is_nan (m.scalar());
  return big<or_is_nan_op> (m); }
TMPL inline bool is_infinite (const Matrix& m) {
  if (is_a_scalar (m)) return is_infinite (m.scalar());
  return !is_nan (m) && big<or_is_infinite_op> (m); }
TMPL inline bool is_fuzz (const Matrix& m) {
  if (is_a_scalar (m)) return is_fuzz (m.scalar());
  return !is_nan (m) && big<or_is_fuzz_op> (m); }
TMPL inline bool is_reliable (const Matrix& m) {
  if (is_a_scalar (m)) return is_reliable (m.scalar());
  return is_reliable (C (0)); }

TMPL inline Matrix change_precision (const Matrix& m, xnat p) {
  return binary_map_scalar<change_precision_op> (m, p); }
TMPL inline xnat precision (const Matrix& m) {
  return big<min_precision_op> (m); }

TMPL inline xint exponent (const Matrix& m) {
  return big<max_exponent_op> (m); }
TMPL inline double magnitude (const Matrix& m) {
  return big<max_magnitude_op> (m); }

/******************************************************************************
* Complex and ball matrices
******************************************************************************/

TMPL Abs_matrix abs (const Matrix& m) { return unary_map<abs_op> (m); }
TMPL Real_matrix Re (const Matrix& m) { return unary_map<Re_op> (m); }
TMPL Real_matrix Im (const Matrix& m) { return unary_map<Im_op> (m); }
TMPL Matrix conj (const Matrix& m) { return unary_map<conj_op> (m); }

TMPL Center_matrix center (const Matrix& m) {
  return unary_map<center_op> (m); }
TMPL Radius_matrix radius (const Matrix& m) {
  return unary_map<radius_op> (m); }
TMPL Center_matrix lower (const Matrix& m) {
  return unary_map<lower_op> (m); }
TMPL Center_matrix upper (const Matrix& m) {
  return unary_map<upper_op> (m); }
TMPL inline Matrix sharpen (const Matrix& m) {
  return unary_map<sharpen_op> (m); }
TMPLK inline Matrix blur (const Matrix& m, const K& x) {
  return binary_map_scalar<blur_op> (m, x); }
template<typename C, typename V, typename C2, typename V2> inline Matrix
blur (const Matrix& m, const matrix<C2,V2>& n) {
  return binary_map<blur_op> (m, n); }

/******************************************************************************
* Explicit matrices
******************************************************************************/

template<typename C> matrix<C>
identity_matrix (const int& n, const format<C>& fm= format<C> ()) {
  return matrix<C> (promote (1, fm), n, n);
}

template<typename C> matrix<C>
fill_matrix (const C& x, const int& r, const int& c) {
  matrix<C> m (promote (0, x), (nat) r, (nat) c);
  for (nat i=0; i<((nat) r); i++)
    for (nat j=0; j<((nat) c); j++)
      m (i, j)= x;
  return m;
}

template<typename C> matrix<C>
tensor_matrix (const vector<C>& v, const vector<C>& w) {
  ASSERT (is_non_scalar (v), "non-scalar vector expected");
  ASSERT (is_non_scalar (w), "non-scalar vector expected");
  matrix<C> m (promote (0, CF(v)), N(v), N(w));
  for (nat i=0; i<N(v); i++)
    for (nat j=0; j<N(w); j++)
      m (i, j)= v[i] * w[j];
  return m;
}

template<typename C> matrix<C>
hilbert_matrix (const int& n, const format<C>& fm= format<C> ()) {
  matrix<C> m (promote (0, fm), n, n);
  for (nat i=0; i<(nat) n; i++)
    for (nat j=0; j<(nat) n; j++)
      m (i, j)= 1 / promote (1 + i + j, fm);
  return m;
}

template<typename C> matrix<C>
vandermonde (const vector<C>& v) {
  ASSERT (is_non_scalar (v), "non-scalar vector expected");
  matrix<C> m (promote (0, CF(v)), N(v), N(v));
  for (nat i=0; i<N(v); i++) {
    C p= promote (1, CF(v));
    for (nat j=0; j<N(v); j++, p *= v[i])
      m (i, j)= p;
  }
  return m;
}

template<typename C> matrix<C>
jordan_matrix (const C& c, const int& n) {
  matrix<C> m (c, n, n);
  for (nat i=0; i<((nat) (n-1)); i++)
    m (i, i+1)= promote (1, c);
  return m;
}

template<typename C> matrix<C>
toeplitz_matrix (const vector<C>& v) {
  ASSERT ((N(v)&1) == 1, "odd dimension expected");
  nat n= (N(v) >> 1) + 1;
  matrix<C> m (promote (0, CF(v)), n, n);
  for (nat i=0; i<n; i++)
    for (nat j=0; j<n; j++)
      m (i, j)= v [n-1 + i-j];
  return m;
}

template<typename C> matrix<C>
hankel_matrix (const vector<C>& v) {
  ASSERT ((N(v)&1) == 1, "odd dimension expected");
  nat n= (N(v) >> 1) + 1;
  matrix<C> m (promote (0, CF(v)), n, n);
  for (nat i=0; i<n; i++)
    for (nat j=0; j<n; j++)
      m (i, j)= v [i+j];
  return m;
}

/******************************************************************************
* Further routines on matrices
******************************************************************************/

TMPL Matrix
horizontal_join (const Matrix& m1, const Matrix& m2) {
  ASSERT (is_non_scalar (m1) || is_non_scalar (m2),
          "non-scalar matrix expected");
  if (!is_non_scalar (m1))
    return horizontal_join (Matrix (m1.scalar(), rows (m2), rows (m2)), m2);
  if (!is_non_scalar (m2))
    return horizontal_join (m1, Matrix (m2.scalar(), rows (m1), rows (m1)));
  ASSERT (rows (m1) == rows (m2), "unequal number of rows");
  Matrix r (promote (0, CF(m1)), rows (m1), cols (m1) + cols (m2));
  for (nat i=0; i<rows(m1); i++) {
    for (nat j=0; j<cols(m1); j++)
      r(i,j)= m1(i,j);
    for (nat j=0; j<cols(m2); j++)
      r(i,j+cols(m1))= m2(i,j);
  }
  return r;
}

TMPL Matrix
vertical_join (const Matrix& m1, const Matrix& m2) {
  ASSERT (is_non_scalar (m1) || is_non_scalar (m2),
          "non-scalar matrix expected");
  if (!is_non_scalar (m1))
    return vertical_join (Matrix (m1.scalar(), cols (m2), cols (m2)), m2);
  if (!is_non_scalar (m2))
    return vertical_join (m1, Matrix (m2.scalar(), cols (m1), cols (m1)));
  ASSERT (cols (m1) == cols (m2), "unequal number of columns");
  Matrix r (promote (0, CF(m1)), rows (m1) + rows (m2), cols (m1));
  for (nat j=0; j<cols(m1); j++) {
    for (nat i=0; i<rows(m1); i++)
      r(i,j)= m1(i,j);
    for (nat i=0; i<rows(m2); i++)
      r(i+rows(m1),j)= m2(i,j);
  }
  return r;
}

TMPL Matrix
transpose (const Matrix & m) {
  typedef implementation<matrix_linear,V> Mat;
  if (is_a_scalar (m)) return m;
  nat nrows= rows (m), ncols= cols (m);
  nat l= aligned_size<C,V> (nrows * ncols);
  C*  r= mmx_formatted_new<C> (l, CF(m));
  Mat::transpose (r, tab (m), nrows, ncols);
  return Matrix (r, ncols, nrows, CF(m));
}

TMPL inline Matrix
range (const Matrix& m, nat r1, nat c1, nat r2, nat c2) {
  typedef implementation<matrix_linear,V> Mat;
  if (is_a_scalar (m)) return Matrix (m.scalar(), r2-r1, c2-c1);
  nat nrows= rows (m), ncols= cols (m);
  nat l= aligned_size<C,V> ((r2-r1) * (c2-c1));
  C*  r= mmx_formatted_new<C> (l, CF(m));
  Mat::get_range (r, tab (m), r1, c1, r2, c2, nrows, ncols);
  return Matrix (r, r2-r1, c2-c1, CF(m));
}

TMPL inline Matrix
delete_row (const Matrix& m, nat r) {
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  ASSERT (r < rows (m), "index out of range");
  Matrix d (promote (0, CF(m)), rows (m) - 1, cols (m));
  for (nat i= 0; i < r; i++)
    for (nat j= 0; j < cols (m); j++) d(i,j)= m(i,j);
  for (nat i= r+1; i < rows (m); i++)
    for (nat j= 0; j < cols (m); j++) d(i-1,j)= m(i,j);
  return d;
}
  
TMPL Matrix
delete_col (const Matrix& m, nat c) {
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  ASSERT (c < cols (m), "index out of range");
  Matrix d (promote (0, CF(m)), rows (m), cols (m) - 1);
  for (nat j= 0; j < c; j++)
    for (nat i= 0; i < rows (m); i++) d(i,j)= m(i,j);
  for (nat j= c+1; j < cols (m); j++)
    for (nat i= 0; i < rows (m); i++) d(i,j-1)= m(i,j);
  return d;
}
  
TMPL Matrix
operator * (const Matrix& m, const Matrix& n) {
  typedef implementation<matrix_multiply,V> Mat;
  if (is_a_scalar (m) || is_a_scalar (n)) {
    if (is_non_scalar (m)) return m * n.scalar();
    if (is_non_scalar (n)) return m.scalar() * n;
    return Matrix (m.scalar() * n.scalar());
  }
  nat mrows= rows (m), mcols= cols (m), nrows= rows(n), ncols= cols(n);
  ASSERT (nrows == mcols, "numbers of rows and columns don't match");
  nat l= aligned_size<C,V> (mrows * ncols);
  C*  r= mmx_formatted_new<C> (l, CF(m));
  Mat::mul (r, tab (m), tab (n), mrows, mcols, ncols);
  return Matrix (r, mrows, ncols, CF(m));
}

template<typename C, typename V, typename W> vector<C,W>
operator * (const matrix<C,V>& m, const vector<C,W>& v) {
  typedef implementation<vector_linear,W> Vec;
  typedef implementation<matrix_linear,V> Mat;
  if (is_a_scalar (m)) return m.scalar() * v;
  ASSERT (is_non_scalar (v), "non-scalar vector expected");
  nat rr= rows (m), cc= cols (m);
  ASSERT (cc == N(v), "sizes don't match");
  nat l= aligned_size<C,W> (rr);
  C*  a= mmx_formatted_new<C> (l, CF(m));
  for (nat i=0; i<rr; i++)
    a[i]= Vec::template vec_binary_big_stride<mul_add_op>
      (tab (m) + Mat::index (i, 0, rr, cc), Mat::index (0, 1, rr, cc),
       seg (v), 1, cc);
  return vector<C,W> (a, rr, l, CF(m));
}

template<typename C, typename V, typename W> vector<C,W>
operator * (const vector<C,W>& v, const matrix<C,V>& m) {
  typedef implementation<vector_linear,W> Vec;
  typedef implementation<matrix_linear,V> Mat;
  if (is_a_scalar (m)) return v * m.scalar();
  ASSERT (is_non_scalar (v), "non-scalar vector expected");
  nat rr= rows (m), cc= cols (m);
  ASSERT (rr == N(v), "sizes don't match");
  nat l= aligned_size<C,W> (cc);
  C*  a= mmx_formatted_new<C> (l, CF(m));
  for (nat i=0; i<cc; i++)
    a[i]= Vec::template vec_binary_big_stride<mul_add_op>
      (seg (v), 1,
       tab (m) + Mat::index (0, i, rr, cc), Mat::index (1, 0, rr, cc), rr);
  return vector<C,W> (a, cc, l, CF(m));
}

/******************************************************************************
* Operations on rows or columns
******************************************************************************/

TMPL vector<C>
row (const Matrix& m, nat i) {
  nat n= cols (m);
  nat l= default_aligned_size<C> (n);
  C*  a= mmx_formatted_new<C> (l, CF(m));
  for (nat j=0; j<n; j++) a[j]= m(i,j);
  return vector<C> (a, n, l, CF(m));
}

TMPL vector<C>
column (const Matrix& m, nat j) {
  nat n= rows (m);
  nat l= default_aligned_size<C> (n);
  C*  a= mmx_formatted_new<C> (l, CF(m));
  for (nat i=0; i<n; i++) a[i]= m(i,j);
  return vector<C> (a, n, l, CF(m));
}

template<typename C> vector<vector<C> >
row_vectors (const matrix<C>& m) {
  vector<vector<C> > v=  fill<vector<C> > (rows (m));
  for (nat i=0; i<rows(m); i++)
    v[i]= row (m, i);
  return v;
}

template<typename C> matrix<C>
row_matrix (const vector<vector<C> >& v) {
  if (N(v) == 0) return matrix<C> (0);
  matrix<C> r (promote (0, get_format1 (CF(v))), N(v), N(v[0]));
  for (nat i=0; N(v)>i; i++) {
    ASSERT (N(v[i]) == N(v[0]), "unequal row lengths");
    for (nat j=0; j<N(v[0]); j++)
      r (i, j)= v[i][j];
  }
  return r;
}

TMPL void
swap_row (Matrix& m, nat i, nat j) {
  typedef implementation<matrix_linear,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  nat mrows= rows (m), mcols= cols (m);
  ASSERT (i < mrows && j < mrows, "out of range");
  Mat::row_swap (tab (m), i, j, mrows, mcols);
}

TMPL void
swap_col (Matrix& m, nat i, nat j) {
  typedef implementation<matrix_linear,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  nat mrows= rows (m), mcols= cols (m);
  ASSERT (i < mcols && j < mcols, "out of range");
  Mat::col_swap (tab (m), i, j, mrows, mcols);
}

TMPL void
reverse_cols (Matrix& m) {
  typedef implementation<matrix_linear,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Mat::col_reverse (tab (m), rows (m), cols (m));
}

#define INLINE_SCAL_LINE_OP(name1,name2)                        \
TMPL void name1 (Matrix& m, C c, nat i) {                       \
  typedef implementation<matrix_linear,V> Mat;                  \
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");     \
  nat mrows=rows(m), mcols=cols(m);                             \
  Mat::name2 (tab(m), c, i, mrows, mcols);                      \
}

#define INLINE_LINES_OP(name1, name2)                           \
TMPL void name1 (Matrix& m, nat i, C ci, nat j, C cj) {         \
  typedef implementation<matrix_linear,V> Mat;                  \
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");     \
  nat mrows=rows(m), mcols=cols(m);                             \
  Mat::name2 (tab(m), i, ci, j, cj, mrows, mcols);              \
}

INLINE_SCAL_LINE_OP(row_mul, row_mul);
INLINE_SCAL_LINE_OP(row_div, row_div);
INLINE_SCAL_LINE_OP(col_mul, col_mul);
INLINE_SCAL_LINE_OP(col_div, col_div);
INLINE_LINES_OP(rows_linsub,row_combine_sub);
INLINE_LINES_OP(cols_linsub,col_combine_sub);

/******************************************************************************
* Permutations
******************************************************************************/

TMPL matrix<C>
as_matrix (const permutation& p, const Format& fm= Format ()) {
  nat n= N(p);
  matrix<C> m (promote (0, fm), n, n);
  for (nat i=0; i<n; i++)
    for (nat j=0; j<n; j++)
      m (i, j) = promote (i == p (j) ? 1 : 0, fm);
  return m;
}

TMPL void
permute_columns (Matrix& m, const permutation& p) {
  // replace m by m * as_matrix (p)
  typedef implementation<matrix_permute,V> Mat;
  if (is_a_scalar (m)) return;
  nat rs= rows (m), cs= cols (m);
  Mat::col_permute (tab (m), seg (*p), rs, cs);
}

TMPL void
permute_rows (Matrix& m, const permutation& p) {
  // replace m by transpose (as_matrix (p)) * m
  typedef implementation<matrix_permute,V> Mat;
  if (is_a_scalar (m)) return;
  nat rs= rows (m), cs= cols (m);
  Mat::row_permute (tab (m), seg (*p), rs, cs);
}

TMPL Matrix
operator * (const Matrix& m, const permutation& p) {
  // returns m * as_matrix (p)
  typedef implementation<matrix_permute,V> Mat;
  if (is_a_scalar (m)) return m;
  nat rs= rows (m), cs= cols (m);
  Matrix dest (promote (0, CF(m)), rs, cs);
  Mat::col_permute (tab (dest), tab (m), seg (*p), rs, cs);
  return dest;
}

TMPL Matrix
operator * (const permutation& p, const Matrix& m) {
  // returns transpose (as_matrix (p)) * m
  typedef implementation<matrix_permute,V> Mat;
  if (is_a_scalar (m)) return m;
  nat rs= rows (m), cs= cols (m);
  Matrix dest (promote (0, CF(m)), rs, cs);
  Mat::row_permute (tab (dest), tab (m), seg (*p), rs, cs);
  return dest;
}

/******************************************************************************
* Echelon forms
******************************************************************************/

TMPL Matrix
column_echelon (const Matrix& m, bool reduced= false) {
  typedef implementation<matrix_echelon,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  C* k= NULL;
  Mat::col_echelon (tab(c), k, rows(m), cols(m), reduced);
  return c;
}

TMPL Matrix
column_reduced_echelon (const Matrix& m) {
  return column_echelon (m, true);
}

TMPL Matrix
column_reduced_echelon (const Matrix& m, permutation& permut) {
  typedef implementation<matrix_echelon,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  C* k= NULL;
  nat* buf= mmx_new<nat> (default_aligned_size<nat> (cols(m)));
  Mat::col_echelon (tab(c), k, rows(m), cols(m), true, buf);
  permut= permutation (vector<nat> (buf, cols(m), format<nat> ()));
  return c;
}

TMPL Matrix
column_echelon (const Matrix& m, Matrix& k, bool reduced= false) {
  typedef implementation<matrix_echelon,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  k= Matrix (promote (0, CF(m)), cols(m), cols(m));
  Mat::col_echelon (tab(c), tab(k), rows(m), cols(m), false);
  return c;
}

TMPL Matrix
column_reduced_echelon (const Matrix& m, Matrix& k) {
  return column_echelon (m, k, true);
}

TMPL inline Matrix
row_echelon (const Matrix& m, bool reduced= false) {
  return transpose (column_echelon (transpose (m), reduced));
}

TMPL Matrix
row_reduced_echelon (const Matrix& m) {
  return row_echelon (m, true);
}

TMPL inline Matrix
row_echelon (const Matrix& m, Matrix& k, bool reduced= false) {
  Matrix c= column_echelon (transpose (m), k, reduced);
  k= transpose (k);
  return transpose (c);
}

TMPL Matrix
row_reduced_echelon (const Matrix& m, Matrix& k) {
  return row_echelon (m, k, true);
}

/******************************************************************************
* System solving
******************************************************************************/

TMPL inline C
det (const Matrix& m) {
  typedef implementation<matrix_determinant,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  ASSERT (cols(m) == rows(m), "square matrix expected");
  C r= promote (1, CF(m));
  Mat::det (r, tab(m), rows(m));
  return r;
}

TMPL inline C
first_minor (const Matrix& m, nat i, nat j) {
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  ASSERT (cols(m) == rows(m), "square matrix expected");
  ASSERT (i < rows(m), "index out of range");
  ASSERT (j < cols(m), "index out of range");
  return det (delete_col (delete_row (m, i), j));
}

TMPL inline C
cofactor (const Matrix& m, nat i, nat j) {
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  ASSERT (cols(m) == rows(m), "square matrix expected");
  ASSERT (i < rows(m), "index out of range");
  ASSERT (j < cols(m), "index out of range");
  C c= first_minor (m, i, j);
  return ((i + j) & 1) ? -c : c;
}

TMPL Matrix
invert (const Matrix& m) {
  typedef implementation<matrix_invert,V> Mat;
  if (is_a_scalar (m)) return Matrix (invert (m.scalar()));
  ASSERT (cols(m) == rows(m), "square matrix expected");
  Matrix a (promote (1, CF(m)), rows(m), cols(m));
  Mat::invert (tab(a), tab(m), rows(m));
  return a;
}

TMPL Matrix operator / (const C& c, const Matrix& m) {
  return c * invert (m); }

/******************************************************************************
* Kernel, image and rank
******************************************************************************/

TMPL Matrix
kernel (const Matrix& m) {
  typedef implementation<matrix_kernel,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix k (promote (0, CF(m)), cols(m), cols(m));
  nat dim= Mat::kernel (tab(k), tab(m), rows(m), cols(m));
  return range (k, 0, 0, cols(m), dim);
}

TMPL Matrix
image (const Matrix& m) {
  typedef implementation<matrix_image,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix k (promote (0, CF(m)), rows(m), rows(m));
  nat dim= Mat::image (tab(k), tab(m), rows(m), cols(m));
  return range (k, 0, 0, rows(m), dim);
}

TMPL nat
rank (const Matrix& m) {
  typedef implementation<matrix_image,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  return Mat::rank (tab(m), rows(m), cols(m));
}

TMPL Matrix
krylov (const Matrix& m, const Matrix& v) {
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  ASSERT (is_square_matrix (m), "square matrix expected");
  Matrix r= image (v);
  Matrix p= m;
  while (true) {
    nat rk= rows (r);
    r= image (vertical_join (r, r*p));
    if (rows (r) <= rk) return r;
    p= p*p;
  }
}

/******************************************************************************
* Gram Schmidt and related decompositions
******************************************************************************/

TMPL inline Matrix
row_orthogonalization (const Matrix& m) {
  typedef implementation<matrix_orthogonalization,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  vector<C> n (promote (0, CF(m)), rows(m));
  Mat::row_orthogonalize (tab(c), rows(m), cols(m), seg(n));
  return c;
}  

TMPL inline Matrix
column_orthogonalization (const Matrix& m) {
  typedef implementation<matrix_orthogonalization,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  vector<C> n (promote (0, CF(m)), cols(m));
  Mat::col_orthogonalize (tab(c), rows(m), cols(m), seg(n));
  return c;
}  

TMPL inline Matrix
row_orthogonalization (const Matrix& m, Matrix& l) {
  typedef implementation<matrix_orthogonalization,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  vector<C> n (promote (0, CF(m)), rows(m));
  l= Matrix (promote (0, CF(m)), rows(m), rows(m));
  Mat::row_orthogonalize (tab(c), rows(m), cols(m), tab(l), seg(n));
  return c;
}  

TMPL inline Matrix
column_orthogonalization (const Matrix& m, Matrix& l) {
  typedef implementation<matrix_orthogonalization,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  vector<C> n (promote (0, CF(m)), cols(m));
  l= Matrix (promote (0, CF(m)), cols(m), cols(m));
  Mat::col_orthogonalize (tab(c), rows(m), cols(m), tab(l), seg(n));
  return c;
}  

TMPL inline Matrix
row_orthonormalization (const Matrix& m) {
  typedef implementation<matrix_orthogonalization,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  Mat::row_orthonormalize (tab(c), rows(m), cols(m));
  return c;
}  

TMPL inline Matrix
column_orthonormalization (const Matrix& m) {
  typedef implementation<matrix_orthogonalization,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  Mat::col_orthonormalize (tab(c), rows(m), cols(m));
  return c;
}  

TMPL inline Matrix
row_orthonormalization (const Matrix& m, Matrix& l) {
  typedef implementation<matrix_orthogonalization,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  l= Matrix (promote (0, CF(m)), rows(m), rows(m));
  Mat::row_orthonormalize (tab(c), rows(m), cols(m), tab(l));
  return c;
}  

TMPL inline Matrix
column_orthonormalization (const Matrix& m, Matrix& l) {
  typedef implementation<matrix_orthogonalization,V> Mat;
  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c= copy (m);
  l= Matrix (promote (0, CF(m)), cols(m), cols(m));
  Mat::col_orthonormalize (tab(c), rows(m), cols(m), tab(l));
  return c;
}  

#undef TMPL_DEF
#undef TMPL
#undef TMPLK
#undef Format
#undef Matrix
#undef Matrix_rep
#undef Abs_matrix
#undef Real_matrix
#undef Center_matrix
#undef Radius_matrix
} // namespace mmx
#endif // __MMX_MATRIX_HPP