Algorithmes pour calculer le determinant d'une matrice

//-------------------------------------------------
// calcul d'un determinant
//-------------------------------------------------
 
#ifndef _UTIL_H_
      #include "util.h"
#endif // _UTIL_H_
 
#ifndef _MATH_H_
      #include "math.h"
#endif // _MATH_H_
 
#include <conio.h>
 
//-------------------------------------------------
typedef double FLT;
//-------------------------------------------------
FLT **CreateMat(int n)
{
FLT **mat; // la matrice
int y; // <y> pour les lignes
 
mat = Malloc(FLT*,n);
for(y=0;y<n;y++)
  {
  mat[y] = Malloc(FLT,n);
  }
 
return mat;
} // CreateMat()
 
 
//-------------------------------------------------
FLT **RandomMat(int n,FLT coeffMin,FLT coeffMax)
{
FLT **mat;
int x,y; // <x> pour les colonnes
            // <y> pour les lignes
 
mat = CreateMat(n);
 
for(x=0;x<n;x++)
  {
  for(y=0;y<n;y++)
    {
    mat[y][x] = (FLT)RandomDouble(coeffMin,coeffMax);
    }
  }
 
return mat;
} // RandomMat()
 
 
//-------------------------------------------------
void DeleteMat(FLT **mat,int n)
{
int y; // <y> pour les lignes
for(y=0;y<n;y++)
  {
  Free(mat[y]);
  }
Free(mat);
} // DeleteMat()
 
 
//-------------------------------------------------
// calcul d'un determinant avec le pivot de Gauss
// la complexite est en n^3
FLT FastDet(FLT **_mat,int n)
{
FLT **mat; // la matrice "tampon", que l'on pourra modifier
FLT det; // la valeur du determinant
int x,y; // <x> pour les colonnes
            // <y> pour les lignes
int i,j; // <i> pour les colonnes
            // <j> pour les lignes
 
// on copie <_mat> dans <mat> car on n'a pas le droit
// de modifier les valeurs de la matrice <_mat>
mat = CreateMat(n);
for(y=0;y<n;y++)
  {
  for(x=0;x<n;x++)
    {
    mat[y][x] = _mat[y][x];
    }
  }
 
// on calcule le determinant par la methode des pivots de Gauss
// on rend la matrice triangulaire superieur tout en conservant
// le determinant de la matrice a chaque iteration
det = 1.;
 
// on balaye les lignes
for(j=0;j<n-1;j++)
  {
  // Le principe est simple, montrons(le sur une exemple 6x6
  // La matrice est sous cette forme :
  //
  // ( a1 * * * * * ) [L1]
  // ( 0 a2 * * * * ) [L2] <-- noms des lignes
  // ( 0 0 a3 * * * ) [L3]
  // ( 0 0 0 #1 * * ) [L4]
  // ( 0 0 0 #2 * * ) [L5]
  // ( 0 0 0 #3 * * ) [L6]
  //
  // Le determinant courant <det> vaut a1.a2.a3
  // Si la matrice est non-singuliere, (i.e. que sont determinant est non-nulle,
  // autrement dit elle est inversible) alors parmis les trois nombres #1 #2 et #3
  // il y en a un qui est non-nul, on prend le plus grand de ceux la en module pour
  // une raison de stabilite numerique.
  // Si #1=#2=#3=0 alors <det> est nul car la matrice est singuliere !
  // Par exemple si le plus grand est #2 alors on echange [L4] et [L5], ce qui va modifier
  // le determinant d'un signe.
  //
  // La matrice sera donc sous cette forme :
  //
  // ( a1 * * * * * ) [L1]
  // ( 0 a2 * * * * ) [L2]
  // ( 0 0 a3 * * * ) [L3]
  // ( 0 0 0 a4 * * ) [L4] (avec a4=#2)
  // ( 0 0 0 #1 * * ) [L5]
  // ( 0 0 0 #3 * * ) [L6]
  //
  // Ensuite il suffit de retrancher #1/a4.[L4] a [L5] et #3/a4.[L4] a [L6] pour obtenir :
  //
  // ( a1 * * * * * ) [L1]
  // ( 0 a2 * * * * ) [L2]
  // ( 0 0 a3 * * * ) [L3]
  // ( 0 0 0 a4 * * ) [L4]
  // ( 0 0 0 0 * * ) [L5]
  // ( 0 0 0 0 * * ) [L6]
  //
  // Le determinant est inchange par cette derniere operation.
  // Le determinant courant vaut maintenant a1.a2.a3.a4
  // A la fin de l'operation on aura <mat> :
  //
  // ( a1 * * * * * ) [L1]
  // ( 0 a2 * * * * ) [L2]
  // ( 0 0 a3 * * * ) [L3]
  // ( 0 0 0 a4 * * ) [L4]
  // ( 0 0 0 0 a5 * ) [L5]
  // ( 0 0 0 0 0 a6 ) [L6]
  //
  // Remarquons qu'il suffit de faire l'operation <n-1> fois
  // car a la fin on aura a6 directement
  //
  // En resume voici les differentes etapes :
  // ( etape 1 ) : on recherche du plus grand coefficient (en module)
  // si il est nul, alors le determinant est 0 et l'algorithme est fini !
  // ( etape 2 ) : on echange des lignes avec eventuellement changement du signe de <det>
  // et le determinant courant est multiplie par ledit coefficient
  // ( etape 3 ) : on retranche la lignes a toutes les lignes du dessous
  // pour annuler les coefficients
  int rankMax,rank;
  FLT coeffMax;
 
  // ( etape 1 )
  rankMax = j;
  for(rank=j+1;rank<n;rank++)
    {
    if(fabs(mat[rankMax][j]) < fabs(mat[rank][j]))
      {
      rankMax = rank;
      }
    }
 
  coeffMax = mat[rankMax][j];
  if(fabs(coeffMax) <= MATH_EPSILON)
    {
    det = 0.;
    goto label_end;
    }
  // ( etape 2 )
  if(rankMax != j)
    {
    for(i=j;i<n;i++)
      {
      FLT tmp;
      tmp = mat[j][i];
      mat[j][i] = mat[rankMax][i];
      mat[rankMax][i] = tmp;
      }
    det *= -1.;
    }
 
  det *= coeffMax;
  // ( etape 3 )
  for(rank=j+1;rank<n;rank++)
    {
    FLT coeff; // c'est dans l'exemple le rapport #3/a4
    coeff = mat[rank][j]/coeffMax;
    for(i=j;i<n;i++)
      {
      mat[rank][i] -= coeff*mat[j][i];
      }
    }
   
  }
 
det *= mat[n-1][n-1];
 
// on libere <mat>
label_end:
DeleteMat(mat,n);
return det;
} // FastDet()
 
 
//-------------------------------------------------
// calcul d'un determinant avec les cofacteurs
// la complexite est en n!
FLT StupidDet(FLT **mat,int n)
{
// cas trivial d'une matrice 1x1
if(1 == n)
  {
  return mat[0][0];
  }
// cas trivial d'une matrice 2x2
//
// | a b |
// | c d | = a.d - c.d
//
else if(2 == n)
  {
  return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
  }
else
  {
  FLT **cof; // matrice qui represente un cofacteur
  FLT det; // le determinant
  int x,y; // <x> pour les colonnes
                // <y> pour les lignes
  int i,j; // <i> pour les colonnes
                // <j> pour les lignes
  FLT sgn;
 
  // on cree une matrice temporaire pour le calcul des cofacteurs
  cof = CreateMat(n-1);
  for(y=0;y<n-1;y++)
    {
    for(x=0;x<n-1;x++)
      {
      cof[y][x] = mat[y][x];
      }
    }
 
  // quelques initialisations
  det = 0.0;
  sgn = +1.0;
  // on balaye la derniere colonne (de bas en haut)
  for(j=n-1;j>=0;j--)
    {
    // Le principe est tres simple, il s'agit de developper le determinant
    // suivant la derniere colonne. Cet algorithme est evidemment recursif
    // puisqu'il faut les sous-determinants
    //
    // Exemple sur une matrice 4x4 :
    //
    // | * * * a1 |
    // | * * * a2 |
    // | * * * a3 |
    // | * * * a4 |
    //
    // A chaque coefficient de la derniere colonne, on associe une sous-dertermiant ("sd") :
    // (notation : '#')
    //
    // ds1 : | * * * a1 | ds2 : | # # # a1 |
    // | # # # a2 | | * * * a2 |
    // | # # # a3 | | # # # a3 |
    // | # # # a4 | | # # # a4 |
    //
    // ds3 : | # # # a1 | ds4 : | # # # a1 |
    // | # # # a2 | | # # # a2 |
    // | * * * a3 | | # # # a3 |
    // | # # # a4 | | * * * a4 |
    //
    // La formule du calcul du determinant par les cofacteurs est donnee par :
    //
    // det = a4.ds4 - a3.ds3 + a2.ds2 - a1.ds1
    //
    // (notons l'alternance des signes)
    //
    // On doit calculer les "ds", donc creer les matrices correspondants, et calculer
    // par recurrence leurs determinants. On remarque facilement que les matrices se
    // ressemblent à (toujours) 2 lignes pres. Donc pour eviter de toujours les recalculer,
    // on les initilise de facon "incrementale".
    //
   
    // on modifie la bonne ligne pour la matrice des cofacteurs
    if(j < n-1)
      {
      for(i=0;i<n-1;i++)
        {
        cof[j][i] = mat[j+1][i];
        }
      }
    // la recurrence en personne !
    det += sgn * mat[j][n-1] * StupidDet(cof,n-1);
    // changement de signe (l'alternance remarquee)
    sgn *= -1.0;
    }
 
  DeleteMat(cof,n-1);
  return det;
  }
} // StupidDet()