synaps/linalg/Det.h

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*      Author(s): B. Mourrain, H. Prieto, GALAAD, INRIA               *
**********************************************************************
$Id: Det.h,v 1.2 2005/09/28 06:40:50 mourrain Exp $
**********************************************************************/
#ifndef SYNAPS_LINALG_DET_H
#define SYNAPS_LINALG_DET_H
//--------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/linalg/MethodName.h>
#include <synaps/linalg/MATRIX.m>
#include <synaps/linalg/Decomp.h>
//--------------------------------------------------------------------
__BEGIN_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
template<class X> bool IsNull(const X & x) {return( x==0);}
//--------------------------------------------------------------------
inline int eps(unsigned int k) 
{
  if(k%2==0) return 1;  else   return -1;
}
/**********************************************************************/
// /*
//    Triangulate a matrix by Bareiss method. The operations are performed 
// inplace. \index{triang, Bareiss} 
// */
// template <class R>
// typename R::value_type triang(const Bareiss & mth, R & M)
// {
//   typedef typename R::value_type value_type;
//   typedef typename R::size_type  size_type;
//   value_type d=1;
//   int eps=1;

//   for(size_type k=0; k<M.nbcol();k++){
//     size_type l=k;
//     while((l<M.nbrow()) && (M(l,k)==0) ) l++;
//     if(l!=M.nbrow()){
//       if(l!=k) {M.swaprow(l,k); eps*=-1;}
//       if(k>0) d=M(k-1,k-1); 
//       for(size_type i=k+1; i<M.nbrow(); i++){
//      value_type c;
//         c= M(i,k);
//      for(size_type j=M.nbcol()-1; j>k; j--) {
//        M(i,j)*=M(k,k);
//        M(i,j)-=c*M(k,j);
//        if ( (M(i,j)!=0) && (d!=0) && (d!=1) ) M(i,j)/=d;
//      }
//      M(i,k)=0;
//       }
//     }
//   }
//   // M(M.nbrow()-1,M.nbcol()-1)*=eps;
//   return M(M.nbrow()-1,M.nbcol()-1)*eps;
// }
//------------------------------------------------------------------------
/* Pseudo-triangulation of the matrix @M@ by rows and columns permutations.
  The matrix @M@ is modified on output, and the number of rows and 
  columns transpositions is returned.
*/
template <class R>
typename R::size_type pseudo_triang( R & M)
{
  typedef typename R::size_type size_type;

  size_type * colind=new size_type[M.nbcol()], * rowind=new size_type[M.nbrow()];
  size_type nbtrans=0;

  for (size_type i=0; i<M.nbrow(); i++) {
    colind[i]=0;
    rowind[i]=0;
  }

  for (size_type i=0; i<M.nbrow(); i++)
    for (size_type j=0; j<M.nbcol(); j++) 
        if (IsNull(M(i,j))) {
          rowind[i]+=1;
          colind[j]+=1;
        }

  bool finish=false;
  while (!finish) {
    finish=true;
    for (size_type i=0; i<M.nbrow()-1; i++)
      {
        if (rowind[i]>rowind[i+1]) {
          M.swaprow(i,i+1);
          swap(rowind[i],rowind[i+1]);
          nbtrans++;
          finish=false;
        }
      }
  }

  finish=false;
  while (!finish) {
    finish=true;
    for (size_type i=0; i<M.nbcol()-1; i++)
      {
        if (colind[i]<colind[i+1]) {
          M.swapcol(i,i+1);
          swap(colind[i],colind[i+1]);
          nbtrans++;
          finish=false;
        }
      }
  }
  return nbtrans;
}
//--------------------------------------------------------------------------
template <class R>
typename R::value_type Det(const R & M){
  using namespace MATRIX;
  assert(M.nbcol()==M.nbrow());

  if(M.nbcol()==2)
    return  M(0,0)*M(1,1)-M(1,0)*M(0,1);
  else if(M.nbcol()==3) {
    typename R::value_type res;
    res = M(0,0)*M(1,1)*M(2,2);
    res -=M(0,0)*M(1,2)*M(2,1);
    res -=M(1,0)*M(0,1)*M(2,2);
    res +=M(1,0)*M(0,2)*M(2,1);
    res +=M(2,0)*M(0,1)*M(1,2);
    res -=M(2,0)*M(0,2)*M(1,1);
    return res;
  } else 
    {
      R N(M.rep());
      return decomp(Bareiss(),N.rep());
    }
}

//---------------------------------------------------------------------------
template <class R>
typename R::value_type Det(const Mixed & mth, R & M){
  assert(M.nbcol()==M.nbrow());
  
  typedef typename R::value_type value_type;
  typedef typename R::size_type  size_type;
  
  value_type delta(0);
  size_type opp(0);

  if (M.nbcol()==1) return M(0,0);
  if (M.nbcol()==2)
    return (M(0,0)*M(1,1)-M(0,1)*M(1,0));

  opp=pseudo_triang(M);
  //  cout <<M<<endl;
  for(size_type k=0; k< M.nbcol();k++) {
    if (!IsNull(M(0,k))) {
      R N(M.nbrow(),M.nbcol());
      for(size_type i=0; i< M.nbrow();i++)
        for(size_type j=0; j< M.nbcol();j++)
          N(i,j) = M(i,j);
      
      submatrix(N,0,0,k,k);

      if ( M.nbcol() >= 4) {
        bool prob=true, altern=false;
        size_type i=0;
        for (size_type j=0;j<N.nbrow();j++) {
          while ((prob) && (i<=j)) {
            if (IsNull(N(j,i))) i++;
            else prob=false;
          }
          if (prob) {
            size_type l=j+1;
            while (l<N.nbrow() && IsNull(N(l,j)) ) l++;
            if (l<N.nbrow()) N.swaprow(j,l);
            else if (((j+1)<N.nbrow()) && (!IsNull(N(j,j+1)))) N.swapcol(j,j+1);
            if (!altern) altern=true;
            else altern=false;
          }
          i=0;
          prob=true;
        }

        Triang<Bareiss>(N);
        if (altern) delta-= value_type(eps(k))*M(0,k)*N(M.nbrow()-2,M.nbcol()-2);
        else delta+= value_type(eps(k))*M(0,k)*N(M.nbrow()-2,M.nbcol()-2);
      }
      else 
        if (M.nbcol()==4) {
          value_type minor1,minor2,minor3;
          minor1 =N(1,1)*N(2,2);
          minor1-=N(1,2)*N(2,1);
          minor1*=N(0,0);
          minor2 =N(1,0)*N(2,2);
          minor2-=N(1,2)*N(2,0);
          minor2*=N(0,1); 
          minor2*=(value_type)(-1);
          minor3 =N(1,0)*N(1,2);
          minor3-=N(1,1)*N(2,0);
          minor3*=N(0,2);
          delta+=value_type(eps(k))*M(0,k)*(minor1-minor2+minor3);
        }
        else 
          if (M.nbcol()==3)
            delta+= value_type(eps(k))*M(0,k)*(N(0,0)*N(1,1)-N(0,1)*N(1,0));
    }
  }
  
  if (opp%2 == 0) 
    return delta;
  else
    return -delta;
}
//--------------------------------------------------------------------------
template<class R>
typename R::value_type Det(const Berko & mth, const R & M)
{
  using MATRIX::trace;
  assert(M.nbcol()==M.nbrow());
  typedef typename R::value_type value_type;
  typedef typename R::size_type  size_type;

  R SM(M);
  value_type d; int eps=1;
  d=trace(SM);
  for (size_type i=2; i<=M.nbcol(); i++)
    {
      SM*=M;
      SM-=(M*d);
      d=trace(SM);
      d/=i; eps*=-1;
    }
  d*=eps;
  return d;
}


__END_NAMESPACE_SYNAPS

#endif // SYNAPS_LINALG_DET_H