/Users/mourrain/Devel/mmx/algebramix/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
*******************************************************************************
* 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 Scalar, typename V>
#define Matrix matrix<Scalar,V>
#define nbrow rows
#define nbcol cols

/*********************************************************************
* Pivoting
**********************************************************************/

TMPL
Scalar bareiss_pivoting(Matrix & A, unsigned & r)
{
  typedef unsigned size_t;
  r=0;
  size_t nr = nbrow(A), nc = nbcol(A);
  size_t i;//,j;
  int s = 1;
  Scalar d=1;// t;
  for(size_t k=0; k< min (nr, nc) && r < nr; k++) 
  {
    for(i=r; i<nr && ( A(i,k) == 0);i++);

    //  for(j=i+1; j<nr; j++) {
    //       if(A(j,k) !=0 && size(A(j, k)) < size(A(i,k))){
    //          i=j;
    //       }
    //   }
    if(i < nr) {
      if( i != r) {
        s *= -1;
        swap_row(A,i,r); 
        // swap_row(A,i,r,k,nc);
        // for(j=k;j<nc;j++) {swap (A(i,j), A(r,j));    }
      }
      for(i = r+1;i<nr;i++) {
        rows_linsub(A,i,A(r,k),r,A(i,k));
        // rows_combine(A,i,A(r,k),r,-A(i,k),0,nr,k+1,nc);
        row_div(A,d,i);
//      for(j=k+1;j<nc;j++) {
//        A(i,j) = (A(i,j)*A(r,k) - A(r,j)*A(i,k));
//        A(i,j) /= d;
//      }
        A(i, k) = 0;
      }
      d = A(r, k);
      r++;
    }
  }
  if(s==-1)
    return -A(nr-1,nc-1);
  else
    return  A(nr-1,nc-1);
}

//--------------------------------------------------------------------
TMPL inline
Scalar bareiss_pivoting(Matrix & A)
{
  unsigned rk=0; return bareiss_pivoting(A,rk);
}

//----------------------------------------------------------------------
TMPL
void bareiss_extended_pivoting(Matrix & A, Matrix& P, unsigned & r)
{
  typedef unsigned size_t;
  r=0;
  size_t nr = nbrow(A), nc = nbcol(A);
  ASSERT(nbrow(A)==nbrow(P),"number of rows of A and P don't match");
  size_t i;//,j;
  int s = 1;
  Scalar d=1;// t;
  for(size_t k=0; k< min (nr, nc) && r < nr; k++) {
    for(i=r; i<nr && (A(i,k) == 0);i++);

    // for(j=i+1; j<nr; j++) {
//       if(A(j,k) !=0){// && size(A(j, k)) < size(A(i,k))){
//              i=j;
//       }
//     }
    if(i < nr) {
      if( i != r) {
        s *= -1;
        swap_row(A,i,r);
        //      swap_row(A,i,r,k,nc);
        swap_row(P,i,r);
        //      swap_row(P,i,r,0,nbcol(P));
        // for(j=k;j<nc;j++) { swap (A(i,j), A(r,j));   }
      }
      for(i = r+1;i<nr;i++) {
        rows_linsub(P,i,A(r,k),r,A(i,k));
        //      rows_combine(P,i,A(r,k),r,-A(i,k),0,nbrow(P),0,nbcol(P));
        rows_linsub(A,i,A(r,k),r,A(i,k));
        //      rows_combine(A,i,A(r,k),r,-A(i,k),0,nr,k+1,nc);

        A(i, k) = 0;
        row_div(A,d,i);
        row_div(P,d,i);

      }
      d= A(r, k);
      r++;
    }
  }
}
//--------------------------------------------------------------------
TMPL inline 
Scalar bareiss_extended_pivoting(Matrix & A, Matrix& P) 
{
  unsigned r=0; bareiss_extended_pivoting(A,P,r);
}

/******************************************************************************
* Applications
******************************************************************************/

TMPL unsigned
bareiss_rank (const Matrix& m) {
  //  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c (copy (m));
  unsigned r;
  bareiss_pivoting (c, r); //0, 0, nbrow(c), nbcol(c));
  return r;
}
//--------------------------------------------------------------------
TMPL Matrix
bareiss_triangulate (const Matrix& m) 
{
  //  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  typedef Scalar C;
  Matrix c (copy (m));
  (void) bareiss_pivoting (c);
  return c;
}

//--------------------------------------------------------------------
TMPL Matrix
bareiss_cokernel (const Matrix& m) {
  //  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c (copy (m));
  unsigned nr =nbrow(c);//, nc = nbcol(c);;
  Matrix p =identity_matrix<Scalar>(nr);
  unsigned rk=0;
  bareiss_extended_pivoting(c,p,rk);
  return range (p, rk, 0, nr, nr);
}
//--------------------------------------------------------------------
TMPL Matrix
bareiss_kernel (const Matrix& m) {
  //  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c =transpose(m);
  unsigned nr =nbrow(c);//, nc = nbcol(c);
  Matrix p =identity_matrix<Scalar>(nr);
  unsigned rk=0;
  bareiss_extended_pivoting(c,p,rk);

  p = range (p, rk, 0, nr, nr);
  //  in_place_transpose(p);
  return transpose(p);
}
//--------------------------------------------------------------------
TMPL Matrix
bareiss_image (const Matrix& m) {
  //  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  Matrix c = transpose(m);
  unsigned 
    //    nr=nbrow(c), 
    nc=nbcol(c), rk;

  bareiss_pivoting(c,rk);
  c = range (c, 0, 0, rk, nc);
  return transpose(c);
}
//--------------------------------------------------------------------
TMPL Matrix
bareiss_coimage (const Matrix& m) {
  //  ASSERT (is_non_scalar (m), "non-scalar matrix expected");

  Matrix c (copy (m));
  unsigned 
    //    nr=nbrow(c), 
    nc=nbcol(c), rk;
  bareiss_pivoting(c,rk);
  return range (c, 0, 0, rk, nc);
}

//--------------------------------------------------------------------
TMPL Scalar 
bareiss_det (const Matrix& m) {
  //  ASSERT (is_non_scalar (m), "non-scalar matrix expected");
  ASSERT (nbcol(m) == nbrow(m), "square matrix expected");
  Matrix c (copy (m));
  unsigned r;
  return bareiss_pivoting (c, r);
}
//--------------------------------------------------------------------

TMPL Matrix
bareiss_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= bareiss_image (v);
  Matrix p= m*r;
  while (true) {
    unsigned rk= nbrow (r);
    r= bareiss_image (horizontal_join (r, p));
    if (nbrow (r) <= rk) return r;
    p = m*p;
  }
}
//====================================================================
}//end of namespace mmx
//====================================================================
#undef Scalar
#undef TMPL
#endif //mmx_matrix_bareiss_hpp