Méthode du gradient conjugué en langage c

Méthode du gradient conjugué en langage c

#include <stdio.h>
#include <stdlib.h>
/// la saisie de la matrice A
void LireMatrice(double (*Q)[10],int *n)
{
int i,j;
printf("\n Entrer l'ordre de la matrice A \n");
scanf("%d",n);
for(i=0;i<*n;i++)
{
for(j=0;j<*n;j++)
{
printf("\n entrer Q[%d][%d] \n ",i+1,j+1);
scanf("%lf",&Q[i][j]);
}
}
}
/// l'affichage de la matrice Q
void AfficheMatrice(double (*Q)[10],int n)
{
int i,j;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
printf("\t %lf ",Q[i][j]);
}
printf("\n");

}
}

/// la saisie de b
void SaisieMatrice_B(double *b,int n)
{
int i;
printf("\n");
for(i=0;i<n;i++)
{
printf("\t Entrez les elements de b[%d] ",i+1);
scanf("%lf",&b[i]);
}
printf("\n");
}
/// l'affichage de b
void Aff_B(double *b,int n)
{
int i;
for(i=0;i<n;i++)
{
printf("%lf\n",b[i]);
}
}
/// la saisie de x0
void Saisie_X0(double (*x)[10],int n)
{
int i;
printf("\n");
for(i=0;i<n;i++)
{
printf("\t Entrez les elements de x0[%d] ",i+1);
scanf("%lf",&x[0][i]);
}
printf("\n");
}
/// l'affichage de x0
void aff_X0(double (*x)[10],int n)
{
int i;
for(i=0;i<n;i++)
{
printf("%lf\n",x[0][i]);
}
}
/// calcul avant d'entrer dans la boucle

///calcule et affichage de Qx0
void Calcu_Aff_QX0(double *Qx,double (*x)[10],double (*Q)[10],int n)
{
int j,i1;
double c;
for(j=0;j<n;j++)
{
c=0;
for(i1=0;i1<n;i1++)
{

c+=Q[j][i1]*x[0][i1];
}

Qx[j]=c;
}
for(j=0;j<n;j++)
{

printf("\n les Ax0 %lf \n ",Qx[j]);
}
}
int main()
{
double Q[10][10],x[10][10],Qx[10],d[10][10],g[10][10],Qd[10];
double s[10],ro2,ro[10],gd,KK,L,gg,b[10];
double c;
int i=0,i1=0,j=0,n;

/// la saisie de la matrice A
LireMatrice(Q,&n);
/// l'affichage de la matrice Q
AfficheMatrice(Q,n);
/// la saisie de b
SaisieMatrice_B(b,n);
/// l'affichage de b
Aff_B(b,n);
/// la saisie de x0
Saisie_X0(x,n);
/// l'affichage de x0
aff_X0(x,n);
/// calcul avant d'entrer dans la boucle
///calcule et affichage de Qx0
Calcu_Aff_QX0(Qx,x,Q,n);
///calcule et affichage de g0
for(j=0;j<n;j++)
{

g[0][j]=Qx[j]-b[j];
}

for(j=0;j<n;j++)
{
printf("\n les g0 %lf \n ",g[0][j]);
}

/// calcule et affichage de d0
for(j=0;j<n;j++)
{
d[0][j]=-g[0][j];

}
for(j=0;j<n;j++)
{

printf("\n les h0 %lf \n ",d[0][j]);
}
/// boucle
for(i=0;i<n;i++)
{
/// calcule de gi*hi
ro[i]=0;
ro2=0;
gd=0;
for(j=0;j<n;j++)
{
gd+=g[i][j]*d[i][j];
}
printf("\n le gd%d est: %lf \n ",i,gd);
/// calcule de Ahi
for(j=0;j<n;j++)
{
Qd[j]=0;
}

for(j=0;j<n;j++)
{c=0;
for(i1=0;i1<n;i1++)
{
c+=Q[j][i1]*d[i][i1];
}
Qd[j]=c;
}
for(j=0;j<n;j++)
{
printf("\n les Qd%d est: %lf \n ",i,Qd[j]);
}

for(j=0;j<n;j++)
{
ro2+=d[i][j]*Qd[j];
}
printf("\n ro2 est %d,%lf \n ",i,ro2);
ro[i]=(-1)*(gd/ro2);
printf("\n le ro%d fin est: %lf \n\n ",i,ro[i]);
double rod[10][10];
/// calcule et affichage de rohi

for(j=0;j<n;j++)
{
rod[i][j]=0;
}
for(i1=0;i1<n;i1++)
{
rod[i][i1]=ro[i]*d[i][i1];
}
for(i1=0;i1<n;i1++)
{
printf("\n roh%d est: %lf \n ",i,rod[i][i1]);
}
/// calcule et affichage de xi+1
for(j=0;j<n;j++)
{
x[i+1][j] =rod[i][j]+x[i][j];
}
/// calcule et affichage de Qxi

/// calcule de Qhi
for(j=0;j<n;j++)
{
Qx[j]=0;
}

for(j=0;j<n;j++)
{
c=0;
for(i1=0;i1<n;i1++)
{

c+=Q[j][i1]*x[i+1][i1];
}
Qx[j]=c;
}
for(j=0;j<n;j++)
{
printf("\n les Qx%d est: %lf \n ",i+1,Qx[j]);
}

for(i1=0;i1<n;i1++)
{
g[i+1][i1]=Qx[i1]-b[i1];

}
for(i1=0;i1<n;i1++)
{
printf("\n g%d est: %lf \n ",i+1,g[i+1][i1]);
}

/// calcule et affichage de GGi et ggi
KK=0;
gg=0;

for(i1=0;i1<n;i1++)
{
KK+=g[i+1][i1]*Qd[i1];
gg+=d[i][i1]*Qd[i1];
}

printf("\n gg%d est: %lf \n ",i,gg);

printf("\n GG%d est: %lf \n ",i,KK);
/// calcule et affichage de Li
L=KK/gg;
printf("\n L%d est: %lf \n ",i,L);
for(i1=0;i1<n;i1++)
{
s[i1]=0;
s[i1]=L*d[i][i1];
s[i1]+=(-1)*g[i+1][i1];
d[i+1][i1] = s[i1];
}
for(i1=0;i1<n;i1++)
{
printf("\n h%d est: %lf \n ",i+1,d[i+1][i1]);
}
}
for(j=0;j<n;j++)
{
printf("\n x%d est %lf : \n ",j+1,x[n][j]);
}
system("pause");
return 0;
}