Souci sur le codage de l'algorithme du gradient conjugué
Bonjour à tous, je vais présenter un peu ce que je compte faire.
Il me faut calculer à partir de la matrice (bande), trouver la solution x par la méthode du gradient conjugué sous l'équation Mx=B (M la matrice, B le second membre un vecteur et x aussi).
Par l'algorithme que je trouve sur le net, j'aimerais savoir si parmi vous qui ont codé cet algo me dire si mon code est bon.
(Pour ceux ou celles qui veulent l'intégralité des fichiers et codes je les ai laissé ici http://www.faceweb.fr/Laplace2D.rar)
Code:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
| Vecteur MatriceBande::GradientConjugue(const Vecteur &B, double epsilon)
{
MatriceBande &M = (*this);
Vecteur r0(N), r1(N), x0(N), p0(N);
double a, b;
int n = 0, Nmax=100000;
r0 = B - M*x0;
p0 = r0;
b = 0;
while (r0.NormeCarree() >= epsilon*epsilon)
{
n++;
if (n>Nmax)
{
cout << " Le gradient conjugue n'a pas converge en moins de "<< Nmax <<" iterations." << endl;
cout << " (||r||^2 = " << r0.NormeCarree()<< " )" << endl;
break;
}
a = r0.NormeCarree() / (M*p0*p0);
x0 = x0 + a*p0;
r1 = r0;
r0 = r0 - a*(M*p0);
p0 = r0 + b*p0;
b = r0.NormeCarree() / r1.NormeCarree();
}
return x0;
} |