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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_naive.hpp
* DESCRIPTION: Naive matrix operations
* COPYRIGHT  : (C) 2007  J. der Hoeven and G. Lecerf and O. Ruatta
*******************************************************************************
* 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_NAIVE__HPP
#define __MMX__MATRIX_NAIVE__HPP
#include <basix/vector.hpp>

namespace mmx {
#define TMPL template<typename C>

/******************************************************************************
* Naive variant for matrices
******************************************************************************/

struct matrix_naive {
  typedef vector_naive Vec;        // Matrices as vectors
  typedef matrix_naive Naive;      // Naive variant
  typedef matrix_naive Positive;   // Multiplication without subtractions
  typedef matrix_naive No_aligned; // Multiplication without memory alignement
  typedef matrix_naive No_simd;    // Variant without SIMD instructions
  typedef matrix_naive No_thread;  // Variant without threads
  typedef matrix_naive No_scaled;  // Variant without preconditioning
};

template<typename C>
struct matrix_variant_helper {
  typedef matrix_naive MV;
};

/******************************************************************************
* Matrix defaults
******************************************************************************/

struct matrix_defaults {};

template<typename V>
struct implementation<matrix_defaults,V,matrix_naive> {
  static const nat def_rows = 0;
  static const nat def_cols = 0;
  static const nat init_rows= 1;
  static const nat init_cols= 1;
}; // implementation<matrix_defaults,V,matrix_naive>

/******************************************************************************
* Inheritance from vectors
******************************************************************************/

template<typename V>
struct implementation<vector_defaults,V,matrix_naive>:
  public implementation<vector_defaults,V,typename V::Vec> {};

template<typename V>
struct implementation<vector_allocate,V,matrix_naive>:
  public implementation<vector_allocate,V,typename V::Vec> {};

template<typename V>
struct implementation<vector_abstractions,V,matrix_naive>:
  public implementation<vector_abstractions,V,typename V::Vec> {};

template<typename V>
struct implementation<vector_abstractions_stride,V,matrix_naive>:
  public implementation<vector_abstractions_stride,V,typename V::Vec> {};

template<typename V>
struct implementation<vector_linear,V,matrix_naive>:
  public implementation<vector_linear,V,typename V::Vec> {};

template<typename C, typename V> inline nat
aligned_size (nat r, nat c) {
  return aligned_size<C,V> (r * c); }

template<typename C> inline nat
default_aligned_size (nat r, nat c) {
  return default_aligned_size<C> (r * c); }

/******************************************************************************
* Vectorial operations
******************************************************************************/

struct matrix_vectorial {};

template<typename V>
struct implementation<matrix_vectorial,V,matrix_naive>:
  public implementation<matrix_defaults,V>
{
  typedef implementation<vector_linear,V> Vec;

static inline nat
index (nat row, nat col, nat rows, nat cols) {
  (void) cols;
  return col * rows + row; }

TMPL
static inline const C&
entry (const C* m, nat row, nat col, nat rows, nat cols) {
  return * (m + index (row, col, rows, cols)); }

/*
static inline nat
index (nat row, nat col, nat rows, nat cols) {
  (void) rows;
  return row * cols + col; }
*/

template<typename Op, typename T, typename C> static inline void
mat_unary_stride (T* dest, nat dest_rs, nat dest_cs,
                  const C* src, nat src_rs, nat src_cs,
                  nat rows, nat cols) {
  for (; cols != 0; dest += dest_cs, src += src_cs, cols--)
    Vec::template vec_unary_stride<Op> (dest, dest_rs, src, src_rs, rows); }

template<typename Op, typename T, typename C1, typename C2> static inline void
mat_binary_stride (T* dest, nat dest_rs, nat dest_cs,
                   const C1* src1, nat src1_rs, nat src1_cs,
                   const C2* src2, nat src2_rs, nat src2_cs,
                   nat rows, nat cols) {
  for (; cols != 0; dest += dest_cs, src1 += src1_cs, src2 += src2_cs, cols--)
    Vec::template vec_binary_stride<Op>
      (dest, dest_rs, src1, src1_rs, src2, src2_rs, rows); }

TMPL static inline void
clear (C* dest, nat r, nat c) {
  Vec::clear (dest, r * c);
}

TMPL static void
set (C* dest, const C& c, nat rows, nat cols) {
  clear (dest, rows, cols);
  nat m= min (rows, cols), plus= index (1, 1, rows, cols);
  while (m != 0) {
    *dest= c;
    dest += plus; m--; } }

TMPL static void
get_range (C* d, const C* s, nat r1, nat c1, nat r2, nat c2, nat r, nat c) {
  nat rr= r2 - r1, cc = c2 - c1;
  mat_unary_stride<id_op>
    (d, index (1, 0, rr, cc), index (0, 1, rr, cc),
     s + index (r1, c1, r, c), index (1, 0, r, c), index (0, 1, r, c),
     rr, cc);
}

TMPL static void
clear_range (C* d, nat r1, nat c1, nat r2, nat c2, nat r, nat c) {
  d += index (r1, c1, r, c);
  C* dd= d;
  nat rp= index (1, 0, r, c), cp= index (0, 1, r, c);
  for (nat j=0; j<(c2-c1); j++, dd += cp, d = dd)
    for (nat i=0; i<(r2-r1); i++, d += rp)
      clear_op::set_op (*d);
}

TMPL static void
set_range (C* d, const C& s, nat r1, nat c1, nat r2, nat c2, nat r, nat c) {
  d += index (r1, c1, r, c);
  C* dd= d;
  nat rp= index (1, 0, r, c), cp= index (0, 1, r, c);
  for (nat j=0; j<(c2-c1); j++, dd += cp, d = dd)
    for (nat i=0; i<(r2-r1); i++, d += rp)
      if (i == j) *d = s;
      else clear_op::set_op (*d);
}

TMPL static void
copy (C* dest, const C* s, nat r, nat c) {
  Vec::copy (dest, s, r * c); }

TMPL static void
transpose (C* dest, const C* src, nat rows, nat cols) {
  mat_unary_stride<id_op>
    (dest, index (0, 1, cols, rows), index (1, 0, cols, rows),
     src,  index (1, 0, rows, cols), index (0, 1, rows, cols),
     rows, cols); }

TMPL static inline void
print (const port& out, const C* s, nat rows, nat cols) {
  out << "[ ";
  for (nat i=0; i < rows; i++) {
    for (nat j=0; j < cols; j++) {
      out << * (s + index (i, j, rows, cols));
      if (j + 1 != cols) out << ", ";
    }
    if (i + 1 != rows) out << ";" << lf;
  }
  out << " ]";
}

}; // implementation<matrix_vectorial,V,matrix_naive>

/******************************************************************************
* Operations on rows and columns
******************************************************************************/

struct matrix_linear {};

template<typename V>
struct implementation<matrix_linear,V,matrix_naive>:
  public implementation<matrix_vectorial,V>
{
  typedef implementation<matrix_vectorial,V> Mat;
  typedef implementation<vector_linear,V> Vec;

template<typename Op, typename C, typename X> static inline void
col_unary_scalar (C* dest, const X& x, nat c1, nat r, nat rr, nat cc) {
  Vec::template vec_unary_scalar_stride<Op>
    (dest + Mat::index (0, c1, rr, cc), Mat::index (1, 0, rr, cc), x, r);
}

template<typename Op, typename C, typename X> static inline void
col_binary_scalar (C* dest, const X& x, nat c1, nat c2, nat r,
                   nat rr, nat cc) {
  Vec::template vec_binary_scalar_stride<Op>
    (dest + Mat::index (0, c1, rr, cc), Mat::index (1, 0, rr, cc),
     dest + Mat::index (0, c2, rr, cc), Mat::index (1, 0, rr, cc), x, r);
}

template<typename Op, typename C> static inline void
col_binary_combine (C* dest, nat c1, nat c2, nat r, nat rr, nat cc) {
  Vec::template vec_binary_combine_stride<Op>
    (dest + Mat::index (0, c1, rr, cc), Mat::index (1, 0, rr, cc),
     dest + Mat::index (0, c2, rr, cc), Mat::index (1, 0, rr, cc), r);
}

template<typename Op, typename C, typename X> static inline void
row_unary_scalar (C* dest, const X& x, nat r1, nat c, nat rr, nat cc) {
  Vec::template vec_unary_scalar_stride<Op>
    (dest + Mat::index (r1, 0, rr, cc), Mat::index (0, 1, rr, cc), x, c);
}

template<typename Op, typename C, typename X> static inline void
row_binary_scalar (C* dest, const X& x, nat r1, nat r2, nat c,
                   nat rr, nat cc) {
  Vec::template vec_binary_scalar_stride<Op>
    (dest + Mat::index (r1, 0, rr, cc), Mat::index (0, 1, rr, cc),
     dest + Mat::index (r2, 0, rr, cc), Mat::index (0, 1, rr, cc), x, c);
}

template<typename Op, typename C> static inline void
row_binary_combine (C* dest, nat r1, nat r2, nat c, nat rr, nat cc) {
  Vec::template vec_binary_combine_stride<Op>
    (dest + Mat::index (r1, 0, rr, cc), Mat::index (0, 1, rr, cc),
     dest + Mat::index (r2, 0, rr, cc), Mat::index (0, 1, rr, cc), c);
}

TMPL static void col_mul (C* s, const C& sc, nat i, nat r, nat c) {
  ASSERT (i<c, "out of range");
  col_unary_scalar<rmul_op> (s, sc, i, r, r, c); }
TMPL static void col_div (C* s, const C& sc, nat i, nat r, nat c) {
  ASSERT (i<c, "out of range");
  col_unary_scalar<rdiv_op> (s, sc, i, r, r, c); }
TMPL static void row_mul (C* s, const C& sc, nat i, nat r, nat c) {
  ASSERT (i<r, "out of range");
  row_unary_scalar<rmul_op> (s, sc, i, c, r, c); }
TMPL static void row_div (C* s, const C& sc, nat i, nat r, nat c) {
  ASSERT (i<r, "out of range");
  row_unary_scalar<rdiv_op> (s, sc, i, c, r, c); }

TMPL static void col_sweep (C* s, nat i, nat j, const C& sc, nat r, nat c) {
  ASSERT (i<c && j<c, "out of range");
  col_binary_scalar<mul_sub_op> (s, sc, i, j, r, r, c); }
TMPL static void row_sweep (C* s, nat i, nat j, const C& sc, nat r, nat c) {
  ASSERT (i<r && j<r, "out of range");
  row_binary_scalar<mul_sub_op> (s, sc, i, j, c, r, c); }
TMPL static void col_sweep (C* s, nat i, nat j, nat r, const C& sc,
                            nat rr, nat cc) {
  col_binary_scalar<mul_sub_op> (s, sc, i, j, r, rr, cc); }
TMPL static void row_sweep (C* s, nat i, nat j, nat c, const C& sc,
                            nat rr, nat cc) {
  row_binary_scalar<mul_sub_op> (s, sc, i, j, c, rr, cc); }

TMPL static void col_swap (C* s, nat i, nat j, nat r, nat c) {
  ASSERT (i<c && j<c, "out of range");
  col_binary_combine<swap_op> (s, i, j, r, r, c); }
TMPL static void row_swap (C* s, nat i, nat j, nat r, nat c) {
  ASSERT (i<r && j<r, "out of range");
  row_binary_combine<swap_op> (s, i, j, c, r, c); }
TMPL static void col_swap (C* s, nat i, nat j, nat r, nat rr, nat cc) {
  col_binary_combine<swap_op> (s, i, j, r, rr, cc); }
TMPL static void row_swap (C* s, nat i, nat j, nat c, nat rr, nat cc) {
  row_binary_combine<swap_op> (s, i, j, c, rr, cc); }

TMPL static void
row_combine_sub (C* s, nat i, const C& si, nat j, const C& sj, nat r, nat c) {
  ASSERT (i<r && j<r, "out of range");
  C* iti= s + Mat::index (i, 0, r, c);
  C* itj= s + Mat::index (j, 0, r, c);
  nat plus= Mat::index (0, 1, r, c);
  while (c>0) {
    *iti= si * (*iti) - sj * (*itj);
    iti += plus; itj += plus; c--; } }

TMPL static void
col_combine_sub (C* s, nat i, const C& si, nat j, const C& sj, nat r, nat c) {
  ASSERT (i<c && j<c, "out of range");
  C* iti= s + Mat::index (0, i, r, c);
  C* itj= s + Mat::index (0, j, r, c);
  nat plus= Mat::index (1, 0, r, c);
  while (r>0) {
    *iti= si * (*iti) - sj * (*itj);
    iti += plus; itj += plus; r--; } }

TMPL static bool
row_is_zero (const C* s, nat i, nat r, nat c) {
  return Vec::template vec_binary_test_scalar_stride<equal_op>
    (s + Mat::index (i, 0, r, c), Mat::index (0, 1, r, c), C (0), c); }

TMPL static bool
col_is_zero (const C* s, nat j, nat r, nat c) {
  return Vec::template vec_binary_test_scalar_stride<equal_op>
    (s + Mat::index (0, j, r, c), Mat::index (1, 0, r, c), C (0), r); }

}; // implementation<matrix_linear,V,matrix_naive>

/******************************************************************************
* Multiplication
******************************************************************************/

template<typename V>
struct matrix_multiply_threshold {};

struct matrix_multiply_base {};
struct matrix_multiply {};

template<typename V>
struct implementation<matrix_multiply_base,V,matrix_naive>:
  public implementation<matrix_linear,V>
{
  typedef implementation<matrix_linear,V> Mat;

template<typename Op, typename C>
struct clear_helper {
  static inline void op (C* d, nat r, nat rr, nat c, nat cc) {
    clear_range (d, 0, 0, r, c, rr, cc); }
};

template<typename C>
struct clear_helper<mul_add_op,C> {
  static inline void op (C* d, nat r, nat rr, nat c, nat cc) {}
};

template<typename Op, typename C> static inline void
clr (C* d, nat r, nat rr, nat c, nat cc) {
  clear_helper<Op,C>::op (d, r, rr, c, cc);
}

template<typename Op, typename D, typename S1, typename S2>
static inline void
mul (D* d, const S1* s1, const S2* s2,
     nat r, nat rr, nat l, nat ll, nat c, nat cc) {
  typedef typename Op::acc_op Acc;
  if (l == 0) clr<Op> (d, r, rr, c, cc);
  else {
    nat ri, ci, li;
    D *dr, *dc;
    const S1 *s1r, *s1c;
    const S2 *s2r, *s2c;
    nat drp = Mat::index (1, 0, rr, cc), dcp = Mat::index (0, 1, rr, cc);
    nat s1rp= Mat::index (1, 0, rr, ll), s1cp= Mat::index (0, 1, rr, ll);
    nat s2rp= Mat::index (1, 0, ll, cc), s2cp= Mat::index (0, 1, ll, cc);
    for (ri= r, dr= d, s1r= s1; ri!=0; ri--, dr += drp, s1r += s1rp)
      for (ci= c, dc= dr, s2c= s2; ci!=0; ci--, dc += dcp, s2c += s2cp) {
        D tmp= *dc;
        for (li= l, s1c= s1r, s2r= s2c; li==l; li--, s1c += s1cp, s2r += s2rp)
          Op ::set_op (tmp, *s1c, *s2r);
        for (                         ; li!=0; li--, s1c += s1cp, s2r += s2rp)
          Acc::set_op (tmp, *s1c, *s2r);
        *dc= tmp;
      }
  }
}

}; // implementation<matrix_multiply_base,V,matrix_naive>

template<typename V>
struct implementation<matrix_multiply,V,matrix_naive>:
  public implementation<matrix_multiply_base,V>
{
  typedef implementation<matrix_multiply_base,V> Mat;

template<typename Op, typename D, typename S1, typename S2>
static inline void
mul (D* d, const S1* s1, const S2* s2,
     nat r, nat rr, nat l, nat ll, nat c, nat cc) {
  Mat::template mul<Op> (d, s1, s2, r, rr, l, ll, c, cc);
}

template<typename D, typename S1, typename S2> static inline void
mul (D* dest, const S1* m1, const S2* m2, nat r, nat l, nat c) {
  Mat::template mul<mul_op> (dest, m1, m2, r, r, l, l, c, c);
}

}; // implementation<matrix_multiply,V,matrix_naive>

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

struct matrix_iterate {};

template<typename V>
struct implementation<matrix_iterate,V,matrix_naive>:
  public implementation<matrix_multiply,V>
{
  typedef implementation<matrix_multiply,V> Mat;

TMPL static inline void
iterate_mul (C** d, const C* s, const C* x,
             nat rs, nat cs, nat rx, nat cx, nat n) {
  // return s^0 x, s^1 x,..., s^{n-1} x in d
  // d must have allocated size n * rs * cx
  VERIFY (cs == rx && rs == cs, "sizes do not match");
  if (rs == 0) return;
  Mat::get_range (d[0], x, 0, 0, rx, cx, rx, cx);
  for (nat i = 1; i < n; i++)
    Mat::mul (d[i], s, d[i-1], rs, cs, cx);
}

TMPL static inline void
project_iterate_mul (C** d, const C* y, const C* s, const C* x,
                     nat ry, nat cy, nat rs, nat cs, nat rx, nat cx, nat n) {
  // return y s^0 x, y s^1 x,..., y s^{n-1} x in d
  VERIFY (cy == rs && rs == cs && cs == rx, "sizes do not match");
  if (n == 0) return;
  nat l= aligned_size<C,V> (rx * cx);
  C* buf= mmx_new<C> (l << 1); 
  C* t1= buf, * t2= t1 + l;
  Mat::get_range (t1, x, 0, 0, rx, cx, rx, cx);
  Mat::mul (d[0], y, t1, ry, cy, cx);
  for (nat i = 1; i < n; i++) {
    Mat::mul (t2, s, t1, rs, cs, cx);
    Mat::mul (d[i], y, t2, ry, cy, cx);
    swap (t1, t2);
  }
  mmx_delete<C> (buf, l << 1);
}

}; // implementation<matrix_iterate,V,matrix_naive>

/******************************************************************************
* Low level action of permutations on matrices
******************************************************************************/

struct matrix_permute {};

template<typename V>
struct implementation<matrix_permute,V,matrix_naive>:
  public implementation<matrix_multiply,V>
{
  typedef implementation<matrix_multiply,V> Mat;
  typedef implementation<vector_linear,V> Vec;

  TMPL static void
  col_permute (C* dest, const C* m, const nat* p, nat r, nat c) {
    // Permute columns of r x c matrix according to length c permutation p
    nat rs= Mat::index (1, 0, r, c);
    for (nat i=0; i<c; i++)
      Vec::template vec_unary_stride<id_op>
        (dest + Mat::index (0, i, r, c), rs,
         m + Mat::index (0, p[i], r, c), rs, r);
  }

  TMPL static void
  row_permute (C* dest, const C* m, const nat* p, nat r, nat c) {
    // Permute rows of r x c matrix according to length c permutation p
    nat cs= Mat::index (0, 1, r, c);
    for (nat i=0; i<r; i++)
      Vec::template vec_unary_stride<id_op>
        (dest + Mat::index (i, 0, r, c), cs,
         m + Mat::index (p[i], 0, r, c), cs, c);
  }

  TMPL static void
  col_permute (C* m, const nat* p, nat r, nat c, bool inv= false) {
    // Permute columns of r x c matrix according to length c permutation p
    // Apply the inverse permutation if inv is true
    nat q[c];
    if (inv) for (nat i=0; i<c; i++) q[i]= p[i];
    else for (nat i=0; i<c; i++) q[p[i]]= i;
    for (nat i=0; i<c; ) {
      if (q[i] == i) i++;
      else {
        nat j= q[i], k= q[j];
        Mat::col_swap (m, i, j, r, c);
        q[i]= k; q[j]= j;
      }
    }
  }

  TMPL static void
  row_permute (C* m, const nat* p, nat r, nat c, bool inv= false) {
    // Permute rows of r x c matrix according to length c permutation p
    // Apply the inverse permutation if inv is true
    nat q[r];
    if (inv) for (nat i=0; i<c; i++) q[i]= p[i];
    else for (nat i=0; i<c; i++) q[p[i]]= i;
    for (nat i=0; i<r; ) {
      if (q[i] == i) i++;
      else {
        nat j= q[i], k= q[j];
        Mat::row_swap (m, i, q[i], r, c);
        q[i]= k; q[j]= j;
      }
    }
  }

}; // implementation<matrix_permute,V,matrix_naive>

/******************************************************************************
* Helper for pivot quality
******************************************************************************/

template<typename C>
struct pivot_helper {
  static inline bool
  better (const C& x1, const C& x2) { return false; }
};

template<typename C> inline bool
better_pivot (const C& x1, const C& x2) {
  return pivot_helper<C>::better (x1, x2);
}

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

struct matrix_echelon {};

template<typename V>
struct implementation<matrix_echelon,V,matrix_naive>:
  public implementation<matrix_multiply,V>
{
  typedef implementation<matrix_multiply,V> Mat;

private:

TMPL static void
col_echelon (C* m, C* k, nat ri, nat ci, nat rf, nat cf, nat rm, nat cm, 
             C& num, C& den, bool reduced, nat* permut)
{
  // Perform column pivoting from (ri, ci) to (rf, cf).
  // In return m is echelonized and k is the transformation matrix,
  // so that the computed m is the product of m * k.
  // If k is NULL then it is not computed. Otherwise k must have
  // allocated size at least cm * cm. If the input matrix is regular,
  // then num/den yields the determinant. We have den=1 if the number
  // of row swaps is even and den=-1 otherwise. Note that the first
  // entries of the echelonized rows need not to be one.
  //
  // if not NULL, permut returns the permutation performed on the columns
  VERIFY ((ri < rm) && (ci < cm) && (rf <= rm) && (cf <= cm), "out of range");
  format<C> fm= get_format (m[0]);
  num= promote (1, fm); den= promote (1, fm);
  nat cb= ci;
  if (permut != NULL) for (nat i= 0; i < cm; i++) permut[i]= i;
  while (ri < rf && ci < cf) {
    // search pivot
    nat i= ri, j= ci;
    nat best_index;
    C   best_val= promote (0, fm);
    for (i= ri; i<rf; i++) {
      best_index= cf;
      for (j= ci; j<cf; j++) {
        C next_val= Mat::entry (m, i, j, rm, cm);
        if (next_val != promote (0, fm))
          if (best_index == cf || better_pivot (next_val, best_val)) {
            best_index= j;
            best_val= next_val;
          }
      }
      if (best_index != cf) { j= best_index; break; }
    }
    
    if (i<rf && j<cf) {
      // put pivot on the left
      if (j != ci) {
        if (k != NULL) Mat::col_swap (k, ci, j, cm, cm);
        Mat::col_swap (m, ci, j, rm, cm);
        den= -den;
        if (permut != NULL) swap (permut[j], permut[ci]); 
      }
      ri= i;
      // column reduction
      C p= Mat::entry (m, ri, ci, rm, cm);
      num *= p;
      for (nat index= reduced ? cb : ci + 1; index < cf; index++) {
        if (index == ci) continue;
        C t= Mat::entry (m, ri, index, rm, cm);
        if (k != NULL) Mat::col_sweep (k, index, ci, t/p, cm, cm);
        Mat::col_sweep (m, index, ci, t/p, rm, cm);
      }
      ri++; ci++;
    }
    else  break;
  }
}

public:

TMPL static inline void
col_echelon (C* m, C* k, nat rm, nat cm, C& num, C& den,
             bool reduced= false, nat* permut= NULL) {
  if (k != NULL && cm != 0) Mat::set (k, promote (1, k[0]), cm, cm);
  col_echelon (m, k, 0, 0, rm, cm, rm, cm, num, den, reduced, permut);
}

TMPL static inline void
col_echelon (C* m, C* k, nat rm, nat cm,
             bool reduced= false, nat* permut= NULL) {
  if (k != NULL && cm != 0) Mat::set (k, promote (1, k[0]), cm, cm);
  C num, den;
  col_echelon (m, k, rm, cm, num, den, reduced, permut);
}

}; // implementation<matrix_echelon,V,matrix_naive>

/******************************************************************************
* Determinant and characteristic polynomial with exact division
******************************************************************************/

struct matrix_determinant {};

template<typename V>
struct implementation<matrix_determinant,V,matrix_naive>:
  public implementation<matrix_echelon,V>
{
  typedef implementation<matrix_echelon,V> Mat;

TMPL static void
det (C& r, const C* m, nat n) {
  if (n == 0) { set_as (r, 1); return; }
  nat len_c= aligned_size<C,V> (n * n);
  C* c= mmx_new<C> (len_c), * k= NULL;
  C num, den;
  Mat::copy (c, m, n, n);
  Mat::col_echelon (c, k, n, n, num, den);
  if (col_is_zero (c, n-1, n, n)) set_as (r, 0);
  else r= num / den;
  mmx_delete<C> (c, len_c);
}

}; // implementation<matrix_determinant,V,matrix_naive>

/******************************************************************************
* Kernel
******************************************************************************/

struct matrix_kernel {};

template<typename V>
struct implementation<matrix_kernel,V,matrix_naive>:
  public implementation<matrix_echelon,V>
{
  typedef implementation<matrix_echelon,V> Mat;

TMPL static nat
kernel (C* k, const C* m, nat rm, nat cm) {
  // return the dimension of the kernel
  VERIFY (rm > 0 && cm > 0, "unexpected empty matrix");
  nat i, j, dim, len_c= aligned_size<C,V> (rm * cm);
  C* c= mmx_new<C> (len_c);
  Mat::copy (c, m, rm, cm);
  Mat::col_echelon (c, k, rm, cm);
  for (i=0; i<cm; i++)
    if (Mat::col_is_zero (c, i, rm, cm)) break;
  mmx_delete<C> (c, len_c);
  dim= cm - i;
  if (dim < cm)
    for (j=0; j < dim; j++)
      Mat::col_swap (k, j, i+j, cm, cm);
    // Note that the next columns are not filled with zeros
  C x;
  for (j=0; j < dim; j++)
    for (i=0; i < cm; i++)
      if ((x= *(k + Mat::index (i, j, cm, cm))) != 0)
        Mat::col_div (k, x, j, cm, cm);
  return dim;
}

}; // implementation<matrix_kernel,V,matrix_naive>

/******************************************************************************
* Image and rank
******************************************************************************/

struct matrix_image {};

template<typename V>
struct implementation<matrix_image,V,matrix_naive>:
  public implementation<matrix_echelon,V>
{
  typedef implementation<matrix_echelon,V> Mat;

TMPL static nat
image (C* k, const C* m, nat rm, nat cm) {
  // return the dimension of the image
  VERIFY (rm > 0 && cm > 0, "unexpected empty matrix");
  nat i, len_c= aligned_size<C,V> (rm * cm);
  C* c= mmx_new<C> (len_c);
  Mat::copy (c, m, rm, cm);
  Mat::col_echelon (c, (C*) NULL, rm, cm);
  for (i=0; i<cm; i++)
    if (Mat::col_is_zero (c, i, rm, cm)) break;
  Mat::get_range (k, c, 0, 0, rm, i, rm, cm);
  // Note that the next columns are not filled with zeros
  mmx_delete<C> (c, len_c);
  return i;
}

TMPL static nat
rank (const C* m, nat rm, nat cm) {
  VERIFY (rm > 0 && cm > 0, "unexpected empty matrix");
  nat i, len_c= aligned_size<C,V> (rm * cm);
  C* c= mmx_new<C> (len_c);
  Mat::copy (c, m, rm, cm);
  Mat::col_echelon (c, (C*) NULL, rm, cm);
  for (i=0; i<cm; i++)
    if (Mat::col_is_zero (c, i, rm, cm)) break;
  mmx_delete<C> (c, len_c);
  return i;
}

}; // implementation<matrix_image,V,matrix_naive>

/******************************************************************************
* Inversion
******************************************************************************/

struct matrix_invert {};

template<typename V>
struct implementation<matrix_invert,V,matrix_naive>:
  public implementation<matrix_echelon,V>
{
  typedef implementation<matrix_echelon,V> Mat;

TMPL static void
invert_lower_triangular (C* inv, const C* m, nat n) {
  for (nat i=0; i<n; i++) {
    C d= Mat::entry (m, i, i, n, n);
    ASSERT (d != 0, "non-invertible matrix");
    for (nat k=0; k<=i; k++) {
      C sum= promote (0, d);
      for (nat j=k; j<i; j++)
        sum -= Mat::entry (m, i, j, n, n) * Mat::entry (inv, j, k, n, n);
      if (i == k) sum += promote (1, d);
      inv[Mat::index (i, k, n, n)]= sum / d;
    }
    for (nat k=i+1; k<n; k++)
      inv[Mat::index (i, k, n, n)]= promote (0, d);
  }
}

TMPL static void
invert (C* k, const C* m, nat n) {
  nat l= aligned_size<C,V> (n * n);
  C* a= mmx_new<C> (l << 1);
  C* b= a + l;
  Mat::copy (k, m, n, n);
  Mat::col_echelon (k, b, n, n); // m * b = k, k lower triangular
  //mmout << "K\n----------\n"; Mat::print (mmout, k, n, n); mmout << "\n";
  //mmout << "B\n----------\n"; Mat::print (mmout, b, n, n); mmout << "\n";
  invert_lower_triangular (a, k, n);
  //mmout << "A\n----------\n"; Mat::print (mmout, a, n, n); mmout << "\n";
  Mat::mul (k, b, a, n, n, n);
  mmx_delete<C> (a, l << 1);
}

}; // implementation<matrix_invert,V,matrix_naive>

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

struct matrix_orthogonalization {};

template<typename V>
struct implementation<matrix_orthogonalization,V,matrix_naive>:
  public implementation<matrix_multiply,V>
{
  typedef implementation<matrix_linear,V> Mat;
  typedef implementation<vector_linear,V> Vec;

TMPL static inline C
row_inn_prod (C* m, nat i, nat j, nat r, nat c) {
  // dot product of rows i and j
  return Vec::template vec_binary_big_stride<mul_add_op> (
    m + Mat::index (i, 0, r, c), Mat::index (0, 1, r, c),
    m + Mat::index (j, 0, r, c), Mat::index (0, 1, r, c), c);
}

TMPL static inline C
col_inn_prod (C* m, nat i, nat j, nat r, nat c) {
  // dot product of columns i and j
  return Vec::template vec_binary_big_stride<mul_add_op> (
    m + Mat::index (0, i, r, c), Mat::index (1, 0, r, c),
    m + Mat::index (0, j, r, c), Mat::index (1, 0, r, c), r);
}

TMPL static void
row_orthogonalize (C* m, nat r, nat c, C* n) {
// In place Gram Schmidt orthogonalization of the row vectors of
// m. The vector n contains the squares of the norms of the row
// vectors of the computed m, ie n[i]= ||m[i,.]||^2.
  for (nat i= 0; i < r; i++) {
    for (nat j= 0; j < i; j++) {
      if (n[j] == 0) continue;
      C mu= row_inn_prod (m, i, j, r, c) / n[j];
      Mat::row_combine_sub (m, i, C(1), j, mu, r, c);
    }
    n[i]= row_inn_prod (m, i, i, r, c);
  }
}

TMPL static void
col_orthogonalize (C* m, nat r, nat c, C* n) {
// Idem for columns
  for (nat i= 0; i < c; i++) {
    for (nat j= 0; j < i; j++) {
      if (n[j] == 0) continue;
      C mu= col_inn_prod (m, i, j, r, c) / n[j];
      Mat::col_combine_sub (m, i, C(1), j, mu, r, c);
    }
    n[i]= col_inn_prod (m, i, i, r, c);
  }
}

TMPL static void
row_orthogonalize (C* m, nat r, nat c, C* l, C* n) {
// In place Gram Schmidt orthogonalization of the column vectors of
// m. The matrix l is lower triangular. In return m * l equals the
// input matrix m. The vector n contains the squares of the norms of
// the column vectors of the computed m, ie n[i]= ||m[i,.]||^2.
  for (nat i= 0; i < r; i++) {
    for (nat j= 0; j < i; j++) {
      if (n[j] == 0) continue;
      C mu= row_inn_prod (m, i, j, r, c) / n[j];
      Mat::row_combine_sub (m, i, C(1), j, mu, r, c);
      *(l + Mat::index (i, j, r, r))= mu;
    }
    n[i]= row_inn_prod (m, i, i, r, c);
    *(l + Mat::index (i, i, r, r))= C(1);
  }
}

TMPL static void
col_orthogonalize (C* m, nat r, nat c, C* l, C* n) {
// Idem for columns.
  for (nat i= 0; i < c; i++) {
    for (nat j= 0; j < i; j++) {
      if (n[j] == 0) continue;
      C mu= col_inn_prod (m, i, j, r, c) / n[j];
      Mat::col_combine_sub (m, i, C(1), j, mu, r, c);
      *(l + Mat::index (j, i, c, c))= mu;
    }
    n[i]= col_inn_prod (m, i, i, r, c);
    *(l + Mat::index (i, i, c, c))= C(1);
  }
}

TMPL static void
row_orthonormalize (C* m, nat r, nat c) {
// In place Gram Schmidt orthonormalization of the row vectors of m.
  for (nat i= 0; i < r; i++) {
    for (nat j= 0; j < i; j++)
      Mat::row_combine_sub (m, i, C(1), j, row_inn_prod (m, i, j, r, c), r, c);
    C a= row_inn_prod (m, i, i, r, c);
    if (a != 0) Mat::row_div (m, sqrt (a), i, r, c);
  }
}

TMPL static void
col_orthonormalize (C* m, nat r, nat c) {
// Idem for columns
  for (nat i= 0; i < c; i++) {
    for (nat j= 0; j < i; j++)
      Mat::col_combine_sub (m, i, C(1), j, col_inn_prod (m, i, j, r, c), r, c);
    C a= col_inn_prod (m, i, i, r, c);
    if (a != 0) Mat::col_div (m, sqrt (a), i, r, c);
  }
}

TMPL static void
row_orthonormalize (C* m, nat r, nat c, C* l) {
// In place Gram Schmidt orthonormalization of the row vectors of m.
// The matrix l is lower triangular. In return m * l equals the
// input matrix m.
  for (nat i= 0; i < r; i++) {
    for (nat j= 0; j < i; j++) {
      C mu= row_inn_prod (m, i, j, r, c);
      Mat::row_combine_sub (m, i, C(1), j, mu, r, c);
      *(l + Mat::index (i, j, r, r))= mu;
    }
    C a= row_inn_prod (m, i, i, r, c);
    if (a != 0) {
      C b= sqrt (a);
      Mat::row_div (m, b, i, r, c);
      *(l + Mat::index (i, i, r, r))= b;
    }
  }
}

TMPL static void
col_orthonormalize (C* m, nat r, nat c, C* l) {
// Idem for columns.
  for (nat i= 0; i < c; i++) {
    for (nat j= 0; j < i; j++) {
      C mu= col_inn_prod (m, i, j, r, c);
      Mat::col_combine_sub (m, i, C(1), j, mu, r, c);
      *(l + Mat::index (j, i, c, c))= mu;
    }
    C a= col_inn_prod (m, i, i, r, c);
    if (a != 0) {
      C b= sqrt (a);
      Mat::col_div (m, b, i, r, c);
      *(l + Mat::index (i, i, c, c))= b;
    }
  }
}

}; // implementation<matrix_orthogonalization,V,matrix_naive>

#undef TMPL
}//namespace mmx
#endif //__MMX__MATRIX_NAIVE__HPP