include/algebramix/matrix_ring_naive.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_ring_naive.hpp
* DESCRIPTION: Naive matrix operations over any ring
* COPYRIGHT  : (C) 2007  G. Lecerf
*******************************************************************************
* 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_RING_NAIVE__HPP
#define __MMX__MATRIX_RING_NAIVE__HPP
#include <algebramix/matrix_naive.hpp>
#include <algebramix/polynomial.hpp>
#include <algebramix/series.hpp>

namespace mmx {
#define TMPL template<typename C>

/******************************************************************************
* Naive matrices over any ring
******************************************************************************/

template<typename V>
struct matrix_ring_naive: public V {
  typedef typename V::Vec Vec;
  typedef matrix_ring_naive<typename V::Naive> Naive;
  typedef matrix_ring_naive<typename V::Positive> Positive;
  typedef matrix_ring_naive<typename V::No_simd> No_simd;
  typedef matrix_ring_naive<typename V::No_thread> No_thread;
  typedef matrix_ring_naive<typename V::No_scaled> No_scaled;
};

template<typename F, typename V, typename W>
struct implementation<F,V,matrix_ring_naive<W> >:
  public implementation<F,V,W> {};

template<typename C, typename V>
struct matrix_variant_helper<polynomial<C,V> > {
  typedef matrix_ring_naive<matrix_naive> MV;
};

template<typename C, typename V>
struct matrix_variant_helper<series<C,V> > {
  typedef matrix_ring_naive<matrix_naive> MV;
};

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

template<typename V, typename W>
struct implementation<matrix_echelon,V,matrix_ring_naive<W> >:
  public implementation<matrix_iterate,V>
{
  typedef implementation<matrix_iterate,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).
  // Comments are the same as in matrix_naive
  VERIFY ((ri < rm) && (ci < cm) && (rf <= rm) && (cf <= cm), "out of range");
  ASSERT (reduced == false, "Reduced echelon form not available over a ring");
  num= C(1); den= C(1);
  if (permut != NULL) for (nat i= 0; i < cm; i++) permut[i]= i;

  // search pivot
  nat i= ri, j= ci;
  nat best_index;
  C   best_val= 0;
  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 != 0)
        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);
      if (permut != NULL) swap (permut[j], permut[ci]); 
      den= -den;
    }
    ri= i;
    // column reduction
    C p= Mat::entry (m, ri, ci, rm, cm);
    num *= p;
    for (nat index = 1; index < cf-ci; index++) {
      den *= p;
      C t= Mat::entry (m, ri, ci + index, rm, cm);
      if (k != NULL) Mat::col_combine_sub (k, ci + index, p, ci, t, cm, cm);
      Mat::col_combine_sub (m, ci + index, p, ci, t, rm, cm);
    }
    // call the pivoting on the non triangular sub-matrix
    if (ri+1 < rf && ci+1 < cf) {
      C x, y;
      col_echelon (m, k, ri+1, ci+1, rf, cf, rm, cm, x, y, reduced, permut);
      num *= x; den *=y;
    }
  }
}

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) Mat::set (k, C(1), 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) Mat::set (k, C(1), cm, cm);
  C num, den;
  col_echelon (m, k, rm, cm, num, den, reduced, permut);
}

}; // implementation<matrix_echelon,V,matrix_ring_naive<W> >

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

template<typename V, typename W>
struct implementation<matrix_determinant,V,matrix_ring_naive<W> >:
  public implementation<matrix_echelon,V>
{
  typedef implementation<vector_linear,V> Vec;
  typedef implementation<matrix_echelon,V> Mat;

private:

TMPL static inline void
berkowitz (C* d, const C* s, nat n) {
  // x^n + d[n-1] x^[n-1] + ... + d[0] is the characteristic
  // polynomial of s.
  // d must have allocated size at least n + 1, and we set d[n]=1
  // This is an implementation of Berkowitz' algorithm (1984).
  if (n == 0) { d[0]= 1; return; }
  d[1]= 1;
  d[0]= - *(s + Mat::index (0, 0, n, n));
  if (n == 1) return;
  nat len_t= aligned_size<C,V> (n + 1);
  nat len_s= aligned_size<C,V> ((n-1) * (n-1));
  nat l= aligned_size<C,V> (n-1);
  C** t= mmx_new<C*> (len_t);
  C* buf= mmx_new<C> (3 * len_t + len_s + 2 * l);
  for (nat j = 0; j < len_t; j++) t[j]= buf + j;
  C* t1= buf + len_t, * t2= buf + 2 * len_t;
  C* loc_s= buf + 3 * len_t;
  C* x= buf + 3 * len_t + len_s, * y= buf + 3 * len_t + len_s + l;
  for (nat i = 1; i < n; i++) {
    // Apply Samuelson's formula to deduce charpoly (s[0..i, 0..i])
    // from charpoly (s[0..i-1, 0..i-1]).
    Mat::get_range (loc_s, s, 0, 0, i, i, n, n);
    Mat::get_range (y, s, i, 0, i+1, i  , n, n);
    Mat::get_range (x, s, 0, i, i  , i+1, n, n);
    Mat::project_iterate_mul (t, y, loc_s, x, 1, i, i, i, i, 1, i);
    for (nat j = 0; j < i; j++) t1[j] = t[j][0];
    C a= * (s + Mat::index (i, i, n, n));
    Mat::mul (t2, t1, d + 1, 1, i, 1);  
    t2[0] += a * d[0];
    for (nat j = 1; j < i; j++) {
      Mat::mul (t2 + j, t1, d + j + 1, 1, i - j, 1);  
      t2[j] += a * d[j] - d[j-1];
    }
    t2[i]  = a * d[i] - d[i-1];
    t2[i+1]= - d[i];
    Vec::neg (d, t2, i+2);
  }
  mmx_delete<C> (buf, 3 * len_t + len_s + 2 * l);
  mmx_delete<C*> (t, len_t);
}

public:

TMPL static inline void
det (C& r, const C* s, nat n) {
  nat l= default_aligned_size<C> (n + 1);
  C* d= mmx_new<C> (l);
  berkowitz (d, s, n);
  //implementation<vector_linear,vector_naive>
  //  ::print (mmout, d, n+1); mmout << lf;
  r= (n & 1) ? - d[0] : d[0];
  mmx_delete<C> (d, l);
}

TMPL static inline void
characteristic_polynomial (C* d, const C* s, nat n) {
  berkowitz (d, s, n);
}

}; // implementation<matrix_determinant,matrix_ring_naive<W> >

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

template<typename V, typename W>
struct implementation<matrix_kernel,V,matrix_ring_naive<W> >:
  public implementation<matrix_echelon,V>
{

TMPL static nat
kernel (C* k, const C* m, nat rm, nat cm) {
  // return the dimension of the kernel
  ERROR ("kernel not available for matrix over any ring");
}

}; // implementation<matrix_kernel,V,matrix_ring_naive<W> >

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

template<typename V, typename W>
struct implementation<matrix_image,V,matrix_ring_naive<W> >:
  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
  ERROR ("image not available for matrix over any ring");
}

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, 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_ring_naive<W> >

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

template<typename V, typename W>
struct implementation<matrix_invert,V,matrix_ring_naive<W> >:
  public implementation<matrix_echelon,V>
{
typedef implementation<matrix_determinant,V> Mat;
typedef implementation<vector_linear,vector_naive> Vec;

TMPL static void
invert (C* k, const C* s, nat n) {
  if (n == 0) return;
  nat l= default_aligned_size<C> (n + 1);
  C* d= mmx_new<C> (l);
  Mat::characteristic_polynomial (d, s, n);
  // Vec::print (mmout, d, n+1); mmout << lf;
  C p= - mmx::invert (d[0]);
  nat ll= aligned_size<C,V> (n * n);
  C* aux= mmx_new<C> (ll);
  Vec::clear (aux, n * n);
  for (nat j= 0; j < n; j++)
    * (k + Mat::index (j, j, n, n))= d[n];
  for (nat i= n-1; i > 0; i--) {
    Mat::mul (aux, k, s, n, n, n);
    Mat::copy (k, aux, n, n);
    for (nat j= 0; j < n; j++)
      * (k + Mat::index (j, j, n, n)) += d[i];
  }
  Vec::mul (k, p, n * n);
  mmx_delete<C> (d, l);
  mmx_delete<C> (aux, ll);
}

}; // implementation<matrix_invert,V,matrix_ring_naive<W> >

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

---