/home/ann2.ext/paccaulf/projsem2/ContourActif/resollu.cpp

Aller à la documentation de ce fichier.

//
// resollu.h : resolution LU 
//

#include "resollu.h"
#include <math.h>

#define TINY 1.0e-20

/* ********************************************************************* */
/* Les deux fonctions suivantes ont ete empruntees a Numerical Recipes.  */
/*                                                                       */
/* Elles permettent la decomposition LU d'une matrice et la resoultion   */
/* d'un systeme d'equation matricielle du type A.x = B                   */
/* ********************************************************************* */
bool ludcmp(double a[][MAT_SIZE], int l, int indx[])
{ // Decomposition L.U de la matrice A
  // A est modifiee par la procedure
  int i,imax=0,j,k;
  double big,dum,sum,temp;
  double *vv;
  bool Inversible = true;

  vv = new double[l];     // Allocation dynamique d'un tableau


  for (i=0;i<l;i++) {
    big=0.0;
    for (j=0;j<l;j++)
      if ((temp=fabs(a[i][j]))>big)
        big=temp;
    if (big==0)
      Inversible= false;  //Erreur: inversion de matrice singuliere
    else
      vv[i]=1.0/big;
  }

  if (Inversible)
    {
      for (j=0;j<l;j++) {
        for (i=0;i<j;i++) {
          sum=a[i][j];
          for (k=0;k<i;k++)
            sum-=a[i][k]*a[k][j];
          a[i][j]=sum;
        }
        big=0.0;
        for (i=j;i<l;i++) {
          sum=a[i][j];
          for (k=0;k<j;k++)
            sum-=a[i][k]*a[k][j];
          a[i][j]=sum;
          if ((dum=vv[i]*fabs(sum))>=big) {
            big=dum;
            imax=i;
          }
        }
        if (j!=imax) {
          for (k=0;k<l;k++) {
            dum=a[imax][k];
            a[imax][k]=a[j][k];
            a[j][k]=dum;
          }

          vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (a[j][j]==0.0)
          a[j][j]= TINY;
        if (j!=l-1) {
          dum=1.0/(a[j][j]);
          for (i=j+1;i<l;i++)
            a[i][j]*=dum;
        }
      }
    }  // if (Inversible)
  delete [] vv;  // Libere l'espace memoire attribue dynamiquement
  return (Inversible);

}

void lubksb(double a[][MAT_SIZE],int l,int indx[],double b[])
{
  // Remontee de substitution apres decomposition
  int i,ii=-1,ip,j;
  double sum;

  for (i=0;i<l;i++) {
    ip=indx[i];
    sum=b[ip];
    b[ip]=b[i];
    if (ii>=0)
      for (j=ii;j<=i-1;j++)
        sum-=a[i][j]*b[j];
    else if (sum)
      ii=i;
    b[i]=sum;
  }
  for (i=l-1;i>=0;i--) {
    sum=b[i];
    for (j=i+1;j<l;j++)
      sum-=a[i][j]*b[j];
    b[i]=sum/a[i][i];
  }
}



/*---------------------------------------------------------------------*/
/* DESCRIPTION: ResoudreSysteme()                                      */
/*              Cette fonction resoud un systeme d'equations matri-    */
/*              cielles de la forme K*B=U                              */
/*              par valeur                                             */
/* PARAMETRES:       K (IN): matrice carree de krigeage                */
/*                   U (IN): vecteur des observations                  */
/*                   B (OUT): vecteur des coefficients recherches      */
/*                   Dimension (IN): dimension de la matrice et des    */
/*                                   des vecteurs                      */
/* VALEUR DE RETOUR: true si le systeme a une solution                 */
/*                   false autrement.                                  */
/* REMARQUE:         Resolution par une methode classique de diagona-  */
/*                   lisation du systeme en descente suivi d'une       */
/*                   remonte pour determiner les coefficients.         */
/*---------------------------------------------------------------------*/
bool ResoudreSysteme(double K[][MAT_SIZE], double U[], double B[], int Dimension)
{
  /* Creation de la Matrice des coefficients B*/
  int MatTampon[MAT_SIZE];
  int Indice;
  bool SolutionExiste;


  /* Inversion de K -> K originale est ecrasee */
  SolutionExiste = ludcmp(K,  Dimension, MatTampon);

  if (SolutionExiste) // S'il y a une solution au systeme matriciel
    {
      for (Indice=0;Indice<Dimension;Indice++)
        B[Indice]= U[Indice];

      lubksb(K, Dimension, MatTampon,B);
    }
  return (SolutionExiste);
}

---