include/algebramix/matrix_bareiss.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : matrix_bareiss.hpp
* DESCRIPTION: Bareiss pivoting for dense matrices
* COPYRIGHT  : (C) 2008 Bernard Mourrain, 2012 Gregoire 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_BAREISS__HPP
#define __MMX__MATRIX_BAREISS__HPP
#include <algebramix/matrix.hpp>

namespace mmx {
#define TMPL template<typename C>

/******************************************************************************
* Matrices for rings with exact division
******************************************************************************/

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

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

/******************************************************************************
* Determinant via Bareiss algorithm
******************************************************************************/

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

TMPL static void
bareiss_pivoting (C& d, C* m, nat r, nat c) {
  nat l= 0, i;
  int s= 1;
  set_as (d, 1);
  for (nat k= 0; k < min (r, c) && l < r; k++) {
    for (i= l; i < r &&  Mat::entry (m, i, k, r, c) == 0; i++) {}
    if (i < r) {
      if (i != l) {
        s *= -1;
        Mat::row_swap (m, i, l, r, c); 
      }
      for (i = l+1; i < r; i++) {
        C si= Mat::entry (m, i, k, r, c);
        C sl= Mat::entry (m, l, k, r, c);
        Mat::row_combine_sub (m, i, sl, l, si, r, c);
        Mat::row_div (m, d, i, r, c);
      }
      d= Mat::entry (m, l, k, r, c);
      l++;
    }
  }
  d= Mat::entry (m, r-1, c-1, r, c);
  if (s == -1) d= -d; }

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;
  Mat::copy (c, m, n, n);
  bareiss_pivoting (r, c, n, n);
  mmx_delete<C> (c, len_c); }
};

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

---