INVERSION, CRAMER ,GAUSS, GAUSS PIVOT PARTIEL, GAUSS PIVOT TOTAL, JORDAN ET LU (CROOT ET DOOLITTLE), CHOLESKY, JACOBIE, GAUSS-SEIDEL.

// mini projet ; le 20/06/2008 
#include <stdio.h>
 #include <math.h>

 float norme(float x[19],int n)
 {
 float ref;int i;
 ref=0;
 for(i=0;i<n;i++) if (x[i]>ref) ref=x[i];
 return(ref);
 }

 // fonction de remplissage

 void remplissage(float a[19][19],float b[19],int n)
 {
 int i,j;
 printf("\n * La matrice A:\n\n");
 for (i=0;i<n;i++)
 {
 printf(" | ");
 for (j=0;j<n;j++) scanf("%f",&a[i][j]);
 }
 printf("\n * Le vecteur B:\n\n");
 for(i=0;i<n;i++)
 {
 printf(" | ");
 scanf("%f",&b[i]);
 }
 }

 // fonction d'affichage système

 void aff_syst(float a[19][19],float b[19],int n)
 {
 int i,j;
 printf("\n\n");
 for (i=0;i<n;i++)
 {
 printf(" [");
 for (j=0;j<n;j++)
 {
 printf(" %.4f ",a[i][j]);
 }
 printf("] [ %.4f ]",b[i]);
 printf("\n");
 }
 }

 // fonction d'affichage matrice

 void aff_mat(float a[19][19],int n)
 {
 int i,j;
 printf("\n\n");
 for (i=0;i<n;i++)
 {
 printf(" [");
 for (j=0;j<n;j++)
 {
 printf(" %.4f ",a[i][j]);
 }
 printf("]\n");
 }
 }

 // Mettre à Zero les elements qui doivent etre des zéro

 void zero(float a[19][19],float b[19],int n)
 {
 int i,j;float eps=1e-4;
 for(i=0;i<n;i++)
 {
 for (j=0;j<n;j++) if (fabs(a[i][j])<eps) a[i][j]=0;
 if (fabs(b[i])<eps) b[i]=0;
 }
 } 


 void comatrices(float a[19][19],float c[19][19],int i,int j,int n)
 {
 int l,k;
 for(l=0;l<n;l++) for(k=0;k<n;k++)
 {
 if ((l<i)&&(k<j)) c[l][k]=a[l][k];
 if ((l>i)&&(k<j)) c[l-1][k]=a[l][k];
 if ((l<i)&&(k>j)) c[l][k-1]=a[l][k];
 if ((l>i)&&(k>j)) c[l-1][k-1]=a[l][k];
 }
 }
 // calcul du determinant

 float det(float a[19][19],int n)
 {
 int k,j;float c[19][19],s;

 k=n-1;

 if(n==0) return(1);

 s=0;
 for(j=0;j<n;j++)
 {
 comatrices(a,c,k,j,n);
 s=s+pow(-1,k+j)*a[k][j]*det(c,k);
 }
 return(s);
 }

 

 // -------------------------------------------------------------
 // Résolution Par Cramer
 // -------------------------------------------------------------

 void cramer(float a[19][19],float b[19],int n)
 {
 float A[19][19],x[19],deter;int i,j,k;

 deter=det(a,n);

 if (deter==0)
 {
 printf("\n => Determinant nul, pas de solutions \n\n");
 }

 for(j=0;j<n;j++)
 {
 for(k=0;k<n;k++)
 {
 if (k==j) for(i=0;i<n;i++) A[i][k]=b[i];
 else for(i=0;i<n;i++) A[i][k]=a[i][k];
 }
 x[j]=det(A,n)/deter;
 }
 printf("\n-------------- Cramer -------------\n");
 printf("\n * Le determiant du systeme : %f \n",deter);
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x[i]);
 }

 // -------------------------------------------------------------
 // Résolution Par Inversion
 // -------------------------------------------------------------

 void inversion(float a[19][19],float b[19],int n)
 {
 float c[19][19],comat[19][19],a_1[19][19],x[19],deter,s;int i,j;

 deter=det(a,n);

 if (deter==0)
 {
 printf("\n => Matrice non inversible, pas de solutions \n\n");
 }

 for(i=0;i<n;i++) for(j=0;j<n;j++)
 {
 comatrices(a,c,i,j,n);
 comat[i][j]=(pow(-1,i+j)/deter)*(det(c,n-1));
 }

 // transposée
 for(i=0;i<n;i++) for(j=0;j<n;j++)
 a_1[i][j]=comat[j][i];

 //résolution
 for(i=0;i<n;i++)
 {s=0;
 for(j=0;j<n;j++) s=s+a_1[i][j]*b[j];
 x[i]=s;
 }

 printf("\n-------------- Inversion -------------\n");
 printf("\n * La matrice inverse A^-1 :");
 aff_mat(a_1,n);
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x[i]);
 }

 // -------------------------------------------------------------
 // Résolution Par Gauss
 // -------------------------------------------------------------

 void gauss(float a[19][19],float b[19],int n)

 {
 float x[19],p,s;int i,j,k;

 for(k=0;k<n-1;k++)
 {
 if (a[k][k]==0)
 {
 printf("\n\n * Un pivot nul ! => methode de Gauss non applicable\n\n");
 }

 //réduction
 for(i=k+1;i<n;i++)
 {
 p=a[i][k]/a[k][k];
 for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
 b[i]=b[i]-p*b[k];
 }
 }

 //Résolution
 for(i=n-1;i>=0;i--)
 {
 s=0;
 for(j=i+1;j<n;j++)s=s+a[i][j]*x[j];
 x[i]=(b[i]-s)/a[i][i];
 }
 zero(a,b,n);
 printf("\n-------------- Gauss -------------\n");
 printf("\n * La matrice reduite :");
 aff_syst(a,b,n);
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x[i]);
 }

 // -------------------------------------------------------------
 // Résolution Par Gauss pivot pariel
 // -------------------------------------------------------------

 void gauss_pivot_partiel(float a[19][19],float b[19],int n)

 {
 float x[19],p,s,ref,temp;int i,j,k,ligne;

 for(k=0;k<n-1;k++)
 {
 // max pour le pivot partiel
 ref=0;
 for(i=k;i<n;i++) if(fabs(a[i][k])>ref) {ref=fabs(a[i][k]);ligne=i;}

 // pivotations
 for(j=k;j<n;j++)
 {temp=a[k][j]; a[k][j]=a[ligne][j] ;a[ligne][j]=temp;}

 temp=b[k]; b[k]=b[ligne]; b[ligne]=temp;

 if (a[k][k]==0)
 {
 printf("\n\n* Un pivot nul ! => methode de Gauss pivot partiel non applicable\n\n");
 }


 //réduction
 for(i=k+1;i<n;i++)
 {
 p=a[i][k]/a[k][k];
 for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
 b[i]=b[i]-p*b[k];
 }
 }

 //Résolution
 for(i=n-1;i>=0;i--)
 {
 s=0;
 for(j=i+1;j<n;j++) s=s+a[i][j]*x[j];
 x[i]=(b[i]-s)/a[i][i];
 }
 zero(a,b,n);
 printf("\n-------- Gauss avec pivot partiel --------\n");
 printf("\n * La matrice reduite :");
 aff_syst(a,b,n);
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x[i]);
 printf("\n");
 }

 // -------------------------------------------------------------
 // Résolution Par Gauss pivot total
 // -------------------------------------------------------------

 void gauss_pivot_total(float a[19][19],float b[19],int n)

 {
 float x[19],p,s,ref,temp;int i,j,k,ligne,colonne,pivot_sol[19],temps;

 // vecteur de pivotation des solutions
 for(i=0;i<n;i++) pivot_sol[i]=i;

 for(k=0;k<n-1;k++)
 {
 // max pour le pivot total
 ref=0;
 for(i=k;i<n;i++) for (j=k;j<n;j++)
 if(fabs(a[i][j])>ref) {ref=fabs(a[i][j]);ligne=i;colonne=j;}

 // pivotations
 for(j=k;j<n;j++) {temp=a[k][j]; a[k][j]=a[ligne][j] ;a[ligne][j]=temp;}

 temp=b[k]; b[k]=b[ligne]; b[ligne]=temp;

 for(i=0;i<n;i++) {temp=a[i][k]; a[i][k]=a[i][colonne] ;a[i][colonne]=temp;}

 // remplissage du vecteur accordé aux pivotations
 temps=pivot_sol[k];
 pivot_sol[k]=pivot_sol[colonne];
 pivot_sol[colonne]=temps;

 if (a[k][k]==0)
 {
 printf("\n\n * Un pivot nul ! => methode de Gauss pivot total non applicable\n\n");
 }

 //réduction
 for(i=k+1;i<n;i++)
 {
 p=a[i][k]/a[k][k];
 for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
 b[i]=b[i]-p*b[k];
 }
 }

 // Résolution
 for(i=n-1;i>=0;i--)
 {
 s=0;
 for(j=i+1;j<n;j++)s=s+a[i][j]*b[j];
 b[i]=(b[i]-s)/a[i][i];
 }

 // pivotation des solutions
 for(i=0;i<n;i++) x[pivot_sol[i]]=b[i];

 zero(a,b,n);
 printf("\n--------- Gauss avec pivot total ---------\n");
 printf("\n * La matrice reduite :");
 aff_syst(a,b,n);
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x[i]);
 printf("\n");
 }

 // -------------------------------------------------------------
 // Résolution Par Jordan
 // -------------------------------------------------------------

 void jordan(float a[19][19],float b[19],int n)

 {
 float p;int i,j,k;

 for(k=0;k<n;k++)
 {
 if (a[k][k]==0)
 {
 printf("\n\n * Un pivot nul ! => methode de Jordan non applicable\n\n");
 }

 p=a[k][k];

 //normalisation
 for (j=k;j<n;j++) a[k][j]=a[k][j]/p;
 b[k]=b[k]/p;

 //réduction
 for(i=0;i<n;i++)
 {
 if (i!=k)
 {
 p=a[i][k];
 for (j=k;j<n;j++) a[i][j]=a[i][j]-p*a[k][j];
 b[i]=b[i]-p*b[k];
 }
 }
 }
 zero(a,b,n);
 printf("\n-------------- Jordan -------------\n");
 printf("\n * La matrice reduite :");
 aff_syst(a,b,n);
 printf("\n * La resolution donne :\n\n");
 for(i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,b[i]);
 printf("\n");
 }

 // -------------------------------------------------------------
 // Résolution Par décomposition LU Doolittle
 // -------------------------------------------------------------

 void LU_doolittle(float a[19][19],float b[19],int n)

 {
 float L[19][19],U[19][19],x[19],y[19],s;int i,j,k,m;

 for (i=0;i<n;i++) for (j=0;j<n;j++)
 {
 if(i==j) L[i][j]=1; else L[i][j]=0;
 U[i][j]=0;
 }

 for (m=0;m<n;m++)
 {
 for (j=m;j<n;j++)
 {
 s=0;
 for (k=0;k<m;k++) s=s+L[m][k]*U[k][j];
 U[m][j]=a[m][j]-s;
 }
 if (U[k][k]==0)
 {
 printf("\n\n * Un mineur nul ! => methode de LU non applicable\n\n");
 }

 for (i=m+1;i<n;i++)
 {
 s=0;
 for (k=0;k<m;k++) s=s+L[i][k]*U[k][m];
 L[i][m]=(a[i][m]-s)/U[m][m];
 }
 }

 // resolution
 for(i=0;i<n;i++)
 {
 s=0;
 for(j=0;j<i;j++) s=s+L[i][j]*y[j];
 y[i]=(b[i]-s)/L[i][i];
 }

 for(i=n-1;i>=0;i--)
 {
 s=0;
 for(j=i+1;j<n;j++)s=s+U[i][j]*x[j];
 x[i]=(y[i]-s)/U[i][i];
 }

 printf("\n------------ LU Doolittle -----------\n");
 printf("\n * A = L * U \n");
 printf("\n * La matrice L :");
 aff_mat(L,n);
 printf("\n * La matrice U :");
 aff_mat(U,n);
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x[i]);
 }

 // -------------------------------------------------------------
 // Résolution Par décomposition LU Crout
 // -------------------------------------------------------------

 void LU_crout(float a[19][19],float b[19],int n)

 {
 float L[19][19],U[19][19],x[19],y[19],s;int i,j,k,m;

 for (i=0;i<n;i++) for (j=0;j<n;j++)
 {
 if(i==j) U[i][j]=1; else U[i][j]=0;
 L[i][j]=0;
 }

 for (m=0;m<n;m++)
 {
 for (i=m;i<n;i++)
 {
 s=0;
 for (k=0;k<m;k++) s=s+L[i][k]*U[k][m];
 L[i][m]=a[i][m]-s;
 }

 if (L[k][k]==0)
 {
 printf("\n\n* Un mineur nul ! => methode de LU non applicable\n\n");
 }

 for (j=m+1;j<n;j++)
 {
 s=0;
 for (k=0;k<m;k++) s=s+L[m][k]*U[k][j];
 U[m][j]=(a[m][j]-s)/L[m][m];
 }
 }

 // resolution
 for(i=0;i<n;i++)
 {
 s=0;
 for(j=0;j<i;j++) s=s+L[i][j]*y[j];
 y[i]=(b[i]-s)/L[i][i];
 }

 for(i=n-1;i>=0;i--)
 {
 s=0;
 for(j=i+1;j<n;j++)s=s+U[i][j]*x[j];
 x[i]=(y[i]-s)/U[i][i];
 }

 printf("\n-------------- LU Crout -------------\n");
 printf("\n * A = L * U \n");
 printf("\n * La matrice L :");
 aff_mat(L,n);
 printf("\n * La matrice U :");
 aff_mat(U,n);
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x[i]);
 }

 // -------------------------------------------------------------
 // Résolution Par décomposition Cholesky
 // -------------------------------------------------------------

 void cholesky(float a[19][19],float b[19],int n)
 {
 float L[19][19],Lt[19][19],x[19],y[19],s,p;int i,j,k;

 // vérification de le symétrie
 for (i=0;i<n;i++) for (j=0;j<n;j++)
 if (a[i][j]!=a[j][i])
 {
 printf("\n\n * Non symetrique => methode de Cholesky non applicable\n\n");
 }

 for (i=0;i<n;i++) for (j=0;j<n;j++) L[i][j]=0;

 for (i=0;i<n;i++)
 {
 s=0;
 for (k=0;k<i;k++) s=s+pow(L[i][k],2);
 p=a[i][i]-s;

 if (p<=0)
 {
 printf("\n\n * Non definie positive => methode de Cholesky non applicable\n\n");
 }

 L[i][i]=sqrt(p);

 for(j=i+1;j<n;j++)
 {
 s=0;
 for (k=0;k<i;k++) s=s+L[i][k]*L[j][k];
 L[j][i]=(a[j][i]-s)/L[i][i];
 }
 }

 for (i=0;i<n;i++) for (j=0;j<n;j++) Lt[i][j]=L[j][i];

 // resolution
 for(i=0;i<n;i++)
 {
 s=0;
 for(j=0;j<i;j++) s=s+L[i][j]*y[j];
 y[i]=(b[i]-s)/L[i][i];
 }

 for(i=n-1;i>=0;i--)
 {
 s=0;
 for(j=i+1;j<n;j++) s=s+Lt[i][j]*x[j];
 x[i]=(y[i]-s)/Lt[i][i];
 }

 printf("\n--------------- Cholesky --------------\n");
 printf("\n * A = L * Lt \n");
 printf("\n * La matrice L :");
 aff_mat(L,n);
 printf("\n * La matrice Lt :");
 aff_mat(Lt,n);
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x[i]);
 }

 // -------------------------------------------------------------
 // Résolution Par Jacobie
 // -------------------------------------------------------------

 void jacobie(float a[19][19],float b[19],int n)
 {
 float x[19],x1[19],x2[19],s,eps=1e-4;int i,j,k,iter=0;

 //initialisation du vecteur
 printf("\n * Veuillez initialisez le vecteur solution : \n\n");
 for (i=0;i<n;i++) {printf(" X(0)[%d]= ",i+1);scanf("%f",&x1[i]);}

 do
 {
 for(i=0;i<n;i++)
 {
 s=0;
 for (j=0;j<n;j++) if (i!=j) s=s+a[i][j]*x1[j];
 x2[i]=(b[i]-s)/a[i][i];
 }
 for (k=0;k<n;k++) {x[k]=fabs(x1[k]-x2[k]);x1[k]=x2[k];}

 iter++;
 }
 while (norme(x,n)>eps) ;

 printf("\n-------------- Jacobie -------------\n");
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x2[i]);
 printf("\n * Apres %d iteration, precision 10^-4. \n",iter);
 }

 // -------------------------------------------------------------
 // Résolution Par Gauss-Seidel
 // -------------------------------------------------------------

 void gauss_seidel(float a[19][19],float b[19],int n)
 {
 float x[19],x1[19],x2[19],s,p,eps=1e-4;int i,j,k,iter=0;

 //initialisation du vecteur
 printf("\n * Veuillez initialisez le vecteur solution : \n\n");
 for (i=0;i<n;i++) {printf(" X(0)[%d]= ",i+1);scanf("%f",&x1[i]);}

 do
 {
 for(i=0;i<n;i++)
 {
 s=0;
 p=0;
 for (j=i+1;j<n;j++) s=s+a[i][j]*x1[j];
 for (j=0;j<i;j++) p=p+a[i][j]*x2[j];
 x2[i]=(b[i]-s-p)/a[i][i];
 }
 for (k=0;k<n;k++) {x[k]=fabs(x1[k]-x2[k]);x1[k]=x2[k];}

 iter++;
 }
 while (norme(x,n)>eps) ;

 printf("\n-------------- Gauss-Saidel -------------\n");
 printf("\n * La resolution donne :\n\n");
 for (i=0;i<n;i++) printf(" X_%d = %f ;\n",i+1,x2[i]);
 printf("\n * Apres %d iteration, precision 10^-4. \n",iter);
 }

 main()
 {
 float a[19][19],b[19];int n,i,j; char valider;

 // choix de la methode
 printf("\n *******************************************************\n");
 printf(" * *\n");
 printf(" * Resolution de A x = B *\n");
 printf(" * *\n");
 printf(" *******************************************************\n");
 printf(" * *\n");
 printf(" * Choix de la methode : *\n");
 printf(" * *\n");
 printf(" * A => Cramer *\n");
 printf(" * B => Inversion *\n");
 printf(" * C => Gauss *\n");
 printf(" * D => Gauss pivot partiel *\n");
 printf(" * E => Gauss pivot total *\n");
 printf(" * F => Jordan *\n");
 printf(" * G => LU Doolittle *\n");
 printf(" * H => LU Crout *\n");
 printf(" * I => Cholesky *\n");
 printf(" * J => Jacobie *\n");
 printf(" * K => Gauss-Seidel *\n");
 printf(" * *\n");
 printf(" *******************************************************\n");
 printf(" * *\n");
 printf(" * X => E X I T *\n");
 printf(" * *\n");
 printf(" *******************************************************\n\n");

 reprendre_choix:

 printf("------ => donner votre choix ------ : ");
 scanf("%s",&valider);

 

 switch(valider)
 {
 case 'A' : {
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
			
			}cramer(a,b,n); break;
 case 'B' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
		inversion(a,b,n); break;	
			} 
 case 'C' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
		gauss(a,b,n); break;	
			} 
 case 'D' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
gauss_pivot_partiel(a,b,n); break;			
			} 
 case 'E' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
			
gauss_pivot_total(a,b,n); break;			} 
 case 'F' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
			
	jordan(a,b,n); break;		} 
 case 'G' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
			
	LU_doolittle(a,b,n); break;		} 
 case 'H' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
			
	LU_crout(a,b,n); break;		} 
 case 'I' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
		cholesky(a,b,n); break;	
			} 
 case 'J' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
	jacobie(a,b,n); break;		
			} 
 case 'K' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
			
gauss_seidel(a,b,n); break;			} 
 case 'X' :{
			
 printf("\n\n * Nombre d'equation-inconues : \n\n * N = ");
 scanf("%d",&n);
 remplissage(a,b,n);
 printf("\n\n * Le systeme est : ");
 aff_syst(a,b,n);
			
printf("\nProgramme interompu par l'utilisateur.......\n\n");

 break;			} 
 default : goto reprendre_choix;
 }
 printf("\n");

 }

// kolsi ahmed 
// ahmedkolsi@gmail.com