synaps/linalg/SousMatMax.h

/*************************************************************************
*      This file is part of the source code of SYNAPS.                   *
*      O. Ruatta, GALAAD, INRIA                                          *
**************************************************************************
$Id: SousMatMax.h,v 1.1 2005/07/11 08:33:34 mourrain Exp $
**********************************************************************/
#ifndef SYNAPS_LINALG_SOUSMATMAX_H
#define SYNAPS_LINALG_SOUSMATMAX_H
//----------------------------------------------------------------------
#include <stdio.h>
#include <list>
#include <cmath>
#include <synaps/init.h>
#include <synaps/bignum.h>
#include <synaps/linalg.h>
#include <synaps/mpol.h>
#include <synaps/mpol/bezoutian.h>



__BEGIN_NAMESPACE_SYNAPS

//-----------------------------------------------------------------------

template<class A, class B>
void PermuteLig(A &M, B &IR, const typename A::size_type a, const typename A::size_type b)
{
  typedef typename A::size_type size_type; 
  typedef typename A::value_type value_type;
  M.swaprow(a,b);
  typename B::size_type tmp;
  tmp = IR[a];
  IR[a] = IR[b];
  IR[b] = tmp;
}


template<class A, class B>
void PermuteCol(A &M, B &IC, const typename A::size_type a, const typename A::size_type b)
{
  typedef typename A::size_type size_type; 
  typedef typename A::value_type value_type;
  M.swapcol(a,b);
  typename B::size_type tmp;
  tmp = IC[a];
  IC[a] = IC[b];
  IC[b] = tmp;
}

template<class A, class B> 
void CherchePivot(A & M, B &IR, B &IC, const typename A::size_type & ind, int & TesT)
{
  // Si on rentre dans cette fonction on a M(ind,ind) = 0
  typedef typename A::size_type size_type;
  typedef typename A::value_type value_type;
  size_type indr = ind;
  size_type indc = ind;
  size_type NbRow = M.nbrow();
  size_type NbCol = M.nbcol();
  TesT = 0; //car M(ind,ind) = 0
  while((TesT == 0) && (indc < NbCol)){
    indr = ind;
    while((TesT == 0) && (indr < NbRow)){
      if (M(indr,indc) != 0){
        TesT = 1;
        PermuteLig(M,IR,ind,indr);
        PermuteCol(M,IC,ind,indc);
      }
      indr++;
    }
    indc++;
  }
  //A la fin de cette procedure, TesT = 0 <=> pas de pivot non nul
}



template<class A, class B>
typename A::size_type FieldGaussElim(const A & MM, B &IndRow , B &IndCol)
{
  typedef typename A::size_type size_type; 
  typedef typename A::value_type value_type;
  A M = MM;
  size_type NbRow, NbCol, SizeMin, ind, Rang;
  NbRow = MM.nbrow();
  NbCol = MM.nbcol();
  SizeMin = std::min(NbRow,NbCol);
  value_type tmp, pivot;
  ind = 0;
  int TesT = 1;
  Rang = 0;
  while((TesT != 0) && (ind < SizeMin)){
    // On suppose que le bloc d'indice [0..ind-1,0..ind-1] est deja triangulaire 
    // Si le pivot suppose est nul
    if(M(ind,ind) == 0){
      //Alors on cherche un nouveau pivot
      CherchePivot(M, IndRow, IndCol, ind, TesT);
    }
    // Si il existe un nouveau pivot
    if(TesT != 0){
      //On stocke le nouveau pivot
      pivot = M(ind,ind);
      // On elimine pqr un pivot brutal
      for(size_type i = ind+1; i < NbRow; i++){
        if(M(i,ind) != 0){
          tmp = value_type(-1)*M(i,ind)/pivot;
          M.addrow(i,ind,tmp);
        }
      }
      // Comme on a un nouveau pivot non nul, on incremente le rang
      Rang++;
    }
    //On passe au pivot suivant si il y a lieu de le faire
    ind++;
  }

  //cout<<"Est-ce bien triangulaire ?\n"<<M<<endl;
  // On retourne le rang de la matrice
  return(Rang);
}


template<class A, class B>
B SousMatMax(const A & M, const B & MM)
{
  typedef typename A::size_type size_type; 
  typedef typename A::value_type value_type;

  int NbRow = M.nbrow();
  int NbCol = M.nbcol();
  
  VectDse<int> IndRow(NbRow);
  VectDse<int> IndCol(NbCol);
  
  for(int i = 0; i < NbRow; i++){
    IndRow[i] = i;
  }
  for(int i = 0; i < NbCol; i++){
    IndCol[i] = i;
  }

  size_type Rang = FieldGaussElim(M, IndRow , IndCol);

  //cout<<"le rang -> "<<Rang<<endl;

  int rang = Rang;
  
  B Resultat(rang, rang);
  
  for(int i=0; i < rang; i++){
    for(int j=0; j < rang; j++){
      Resultat(i,j) = MM(IndRow[i],IndCol[j]);
    }
  }

  //cout<<"le resultat de SousMatMax ->"<<endl<<Resultat<<endl;
  
  return Resultat;
}

__END_NAMESPACE_SYNAPS

#endif // SYNAPS_LINALG_SOUSMATMAX_H