1 //static const char* sccsid = "@(#)math_Uzawa.cxx 3.2 95/01/10"; // Do not delete this line. Used by sccs.
5 // Ce programme utilise l algorithme d Uzawa pour resoudre un systeme
6 // dans le cas de contraintes. Si ce sont des contraintes d egalite, la
7 // resolution est directe.
8 // Le programme ci-dessous utilise la methode de Crout pour trouver
9 // l inverse d une matrice symetrique (Le gain est d environ 30% par
10 // rapport a Gauss.). Les calculs sur les matrices sont faits avec chaque
11 // coordonnee car il est plus long d utiliser les methodes deja ecrites
12 // de la classe Matrix avec un passage par valeur.
15 #define No_Standard_RangeError
16 #define No_Standard_OutOfRange
17 #define No_Standard_DimensionError
20 #include <math_Uzawa.ixx>
21 #include <math_Crout.hxx>
22 #include <Standard_DimensionError.hxx>
25 math_Uzawa::math_Uzawa(const math_Matrix& Cont, const math_Vector& Secont,
26 const math_Vector& StartingPoint,
27 const Standard_Real EpsLix, const Standard_Real EpsLic,
28 const Standard_Integer NbIterations) :
29 Resul(1, Cont.ColNumber()),
30 Erruza(1, Cont.ColNumber()),
31 Errinit(1, Cont.ColNumber()),
32 Vardua(1, Cont.RowNumber()),
33 CTCinv(1, Cont.RowNumber(),
36 Perform(Cont, Secont, StartingPoint, Cont.RowNumber(), 0, EpsLix,
37 EpsLic, NbIterations);
40 math_Uzawa::math_Uzawa(const math_Matrix& Cont, const math_Vector& Secont,
41 const math_Vector& StartingPoint,
42 const Standard_Integer Nce, const Standard_Integer Nci,
43 const Standard_Real EpsLix, const Standard_Real EpsLic,
44 const Standard_Integer NbIterations) :
45 Resul(1, Cont.ColNumber()),
46 Erruza(1, Cont.ColNumber()),
47 Errinit(1, Cont.ColNumber()),
48 Vardua(1, Cont.RowNumber()),
49 CTCinv(1, Cont.RowNumber(),
50 1, Cont.RowNumber()) {
52 Perform(Cont, Secont, StartingPoint, Nce, Nci, EpsLix, EpsLic, NbIterations);
57 void math_Uzawa::Perform(const math_Matrix& Cont, const math_Vector& Secont,
58 const math_Vector& StartingPoint,
59 const Standard_Integer Nce, const Standard_Integer Nci,
60 const Standard_Real EpsLix, const Standard_Real EpsLic,
61 const Standard_Integer NbIterations) {
63 Standard_Integer i,j, k;
65 Standard_Real Normat, Normli, Xian, Xmax=0, Xmuian;
66 Standard_Real Rho, Err, Err1, ErrMax=0, Coef = 1./Sqrt(2.);
67 Standard_Integer Nlig = Cont.RowNumber();
68 Standard_Integer Ncol = Cont.ColNumber();
70 Standard_DimensionError_Raise_if((Secont.Length() != Nlig) ||
71 ((Nce+Nci) != Nlig), " ");
73 // Calcul du vecteur Cont*X0 - D: (erreur initiale)
74 //==================================================
76 for (i = 1; i<= Nlig; i++) {
77 Errinit(i) = Cont(i,1)*StartingPoint(1)-Secont(i);
78 for (j = 2; j <= Ncol; j++) {
79 Errinit(i) += Cont(i,j)*StartingPoint(j);
83 if (Nci == 0) { // cas de resolution directe
84 NbIter = 1; //==========================
85 // Calcul de Cont*T(Cont)
86 for (i = 1; i <= Nlig; i++) {
87 for (j =1; j <= i; j++) { // a utiliser avec Crout.
88 // for (j = 1; j <= Nlig; j++) { // a utiliser pour Gauss.
89 CTCinv(i,j)= Cont(i,1)*Cont(j,1);
90 for (k = 2; k <= Ncol; k++) {
91 CTCinv(i,j) += Cont(i,k)*Cont(j,k);
95 // Calcul de l inverse de CTCinv :
96 //================================
97 // CTCinv = CTCinv.Inverse(); // utilisation de Gauss.
98 math_Crout inv(CTCinv); // utilisation de Crout.
99 CTCinv = inv.Inverse();
100 for (i =1; i <= Nlig; i++) {
101 scale = CTCinv(i,1)*Errinit(1);
102 for (j = 2; j <= i; j++) {
103 scale += CTCinv(i,j)*Errinit(j);
105 for (j =i+1; j <= Nlig; j++) {
106 scale += CTCinv(j,i)*Errinit(j);
110 for (i = 1; i <= Ncol; i++) {
111 Erruza(i) = -Cont(1,i)*Vardua(1);
112 for (j = 2; j <= Nlig; j++) {
113 Erruza(i) -= Cont(j,i)*Vardua(j);
116 // restitution des valeurs calculees:
117 //===================================
118 Resul = StartingPoint + Erruza;
119 Done = Standard_True;
121 } // Fin de la resolution directe.
122 //==============================
125 // Initialisation des variables duales.
126 //=====================================
127 for (i = 1; i <= Nlig; i++) {
136 // Calcul du coefficient Rho:
137 //===========================
139 for (i = 1; i <= Nlig; i++) {
140 Normli = Cont(i,1)*Cont(i,1);
141 for (j = 2; j <= Ncol; j++) {
142 Normli += Cont(i,j)*Cont(i,j);
148 // Boucle des iterations de la methode d Uzawa.
149 //=============================================
150 for (NbIter = 1; NbIter <= NbIterations; NbIter++) {
151 for (i = 1; i <= Ncol; i++) {
153 Erruza(i) = -Cont(1,i)*Vardua(1);
154 for(j =2; j <= Nlig; j++) {
155 Erruza(i) -= Cont(j,i)*Vardua(j);
159 Xmax = Abs(Erruza(i) - Xian);
161 Xmax = Max(Xmax, Abs(Erruza(i)-Xian));
165 // Calcul de Xmu a l iteration NbIter et evaluation de l erreur sur
166 // la verification des contraintes.
167 //=================================================================
168 for (i = 1; i <= Nlig; i++) {
169 Err = Cont(i,1)*Erruza(1) + Errinit(i);
170 for (j = 2; j <= Ncol; j++) {
171 Err += Cont(i,j)*Erruza(j);
174 Vardua(i) += Rho * Err;
179 Vardua(i) = Max(0.0, Vardua(i)+ Rho*Err);
180 Err1 = Abs(Vardua(i) - Xmuian);
185 ErrMax = Max(ErrMax, Err1);
189 if (Xmax <= EpsLix) {
190 if (ErrMax <= EpsLic) {
191 // cout <<"Convergence atteinte dans Uzawa"<<endl;
192 Done = Standard_True;
195 // cout <<"convergence non atteinte pour le probleme dual"<<endl;
196 Done = Standard_False;
199 // Restitution des valeurs calculees
200 //==================================
201 Resul = StartingPoint + Erruza;
202 Done = Standard_True;
207 } // fin de la boucle d iterations.
208 Done = Standard_False;
213 void math_Uzawa::Duale(math_Vector& V) const
218 void math_Uzawa::Dump(Standard_OStream& o) const {
222 o << " Status = Done \n";
223 o << " Number of iterations = " << NbIter << endl;
224 o << " The solution vector is: " << Resul << endl;
227 o << " Status = not Done \n";