1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
17 // Ce programme utilise l algorithme d Uzawa pour resoudre un systeme
18 // dans le cas de contraintes. Si ce sont des contraintes d egalite, la
19 // resolution est directe.
20 // Le programme ci-dessous utilise la methode de Crout pour trouver
21 // l inverse d une matrice symetrique (Le gain est d environ 30% par
22 // rapport a Gauss.). Les calculs sur les matrices sont faits avec chaque
23 // coordonnee car il est plus long d utiliser les methodes deja ecrites
24 // de la classe Matrix avec un passage par valeur.
27 #define No_Standard_RangeError
28 #define No_Standard_OutOfRange
29 #define No_Standard_DimensionError
33 #include <math_Crout.hxx>
34 #include <math_Matrix.hxx>
35 #include <math_Uzawa.hxx>
36 #include <Standard_ConstructionError.hxx>
37 #include <Standard_DimensionError.hxx>
38 #include <StdFail_NotDone.hxx>
40 math_Uzawa::math_Uzawa(const math_Matrix& Cont, const math_Vector& Secont,
41 const math_Vector& StartingPoint,
42 const Standard_Real EpsLix, const Standard_Real EpsLic,
43 const Standard_Integer NbIterations) :
44 Resul(1, Cont.ColNumber()),
45 Erruza(1, Cont.ColNumber()),
46 Errinit(1, Cont.ColNumber()),
47 Vardua(1, Cont.RowNumber()),
48 CTCinv(1, Cont.RowNumber(),
51 Perform(Cont, Secont, StartingPoint, Cont.RowNumber(), 0, EpsLix,
52 EpsLic, NbIterations);
55 math_Uzawa::math_Uzawa(const math_Matrix& Cont, const math_Vector& Secont,
56 const math_Vector& StartingPoint,
57 const Standard_Integer Nce, const Standard_Integer Nci,
58 const Standard_Real EpsLix, const Standard_Real EpsLic,
59 const Standard_Integer NbIterations) :
60 Resul(1, Cont.ColNumber()),
61 Erruza(1, Cont.ColNumber()),
62 Errinit(1, Cont.ColNumber()),
63 Vardua(1, Cont.RowNumber()),
64 CTCinv(1, Cont.RowNumber(),
65 1, Cont.RowNumber()) {
67 Perform(Cont, Secont, StartingPoint, Nce, Nci, EpsLix, EpsLic, NbIterations);
72 void math_Uzawa::Perform(const math_Matrix& Cont, const math_Vector& Secont,
73 const math_Vector& StartingPoint,
74 const Standard_Integer Nce, const Standard_Integer Nci,
75 const Standard_Real EpsLix, const Standard_Real EpsLic,
76 const Standard_Integer NbIterations) {
78 Standard_Integer i,j, k;
80 Standard_Real Normat, Normli, Xian, Xmax=0, Xmuian;
81 Standard_Real Rho, Err, Err1, ErrMax=0, Coef = 1./Sqrt(2.);
82 Standard_Integer Nlig = Cont.RowNumber();
83 Standard_Integer Ncol = Cont.ColNumber();
85 Standard_DimensionError_Raise_if((Secont.Length() != Nlig) ||
86 ((Nce+Nci) != Nlig), " ");
88 // Calcul du vecteur Cont*X0 - D: (erreur initiale)
89 //==================================================
91 for (i = 1; i<= Nlig; i++) {
92 Errinit(i) = Cont(i,1)*StartingPoint(1)-Secont(i);
93 for (j = 2; j <= Ncol; j++) {
94 Errinit(i) += Cont(i,j)*StartingPoint(j);
98 if (Nci == 0) { // cas de resolution directe
99 NbIter = 1; //==========================
100 // Calcul de Cont*T(Cont)
101 for (i = 1; i <= Nlig; i++) {
102 for (j =1; j <= i; j++) { // a utiliser avec Crout.
103 // for (j = 1; j <= Nlig; j++) { // a utiliser pour Gauss.
104 CTCinv(i,j)= Cont(i,1)*Cont(j,1);
105 for (k = 2; k <= Ncol; k++) {
106 CTCinv(i,j) += Cont(i,k)*Cont(j,k);
110 // Calcul de l inverse de CTCinv :
111 //================================
112 // CTCinv = CTCinv.Inverse(); // utilisation de Gauss.
113 math_Crout inv(CTCinv); // utilisation de Crout.
114 CTCinv = inv.Inverse();
115 for (i =1; i <= Nlig; i++) {
116 scale = CTCinv(i,1)*Errinit(1);
117 for (j = 2; j <= i; j++) {
118 scale += CTCinv(i,j)*Errinit(j);
120 for (j =i+1; j <= Nlig; j++) {
121 scale += CTCinv(j,i)*Errinit(j);
125 for (i = 1; i <= Ncol; i++) {
126 Erruza(i) = -Cont(1,i)*Vardua(1);
127 for (j = 2; j <= Nlig; j++) {
128 Erruza(i) -= Cont(j,i)*Vardua(j);
131 // restitution des valeurs calculees:
132 //===================================
133 Resul = StartingPoint + Erruza;
134 Done = Standard_True;
136 } // Fin de la resolution directe.
137 //==============================
140 // Initialisation des variables duales.
141 //=====================================
142 for (i = 1; i <= Nlig; i++) {
151 // Calcul du coefficient Rho:
152 //===========================
154 for (i = 1; i <= Nlig; i++) {
155 Normli = Cont(i,1)*Cont(i,1);
156 for (j = 2; j <= Ncol; j++) {
157 Normli += Cont(i,j)*Cont(i,j);
163 // Boucle des iterations de la methode d Uzawa.
164 //=============================================
165 for (NbIter = 1; NbIter <= NbIterations; NbIter++) {
166 for (i = 1; i <= Ncol; i++) {
168 Erruza(i) = -Cont(1,i)*Vardua(1);
169 for(j =2; j <= Nlig; j++) {
170 Erruza(i) -= Cont(j,i)*Vardua(j);
174 Xmax = Abs(Erruza(i) - Xian);
176 Xmax = Max(Xmax, Abs(Erruza(i)-Xian));
180 // Calcul de Xmu a l iteration NbIter et evaluation de l erreur sur
181 // la verification des contraintes.
182 //=================================================================
183 for (i = 1; i <= Nlig; i++) {
184 Err = Cont(i,1)*Erruza(1) + Errinit(i);
185 for (j = 2; j <= Ncol; j++) {
186 Err += Cont(i,j)*Erruza(j);
189 Vardua(i) += Rho * Err;
194 Vardua(i) = Max(0.0, Vardua(i)+ Rho*Err);
195 Err1 = Abs(Vardua(i) - Xmuian);
200 ErrMax = Max(ErrMax, Err1);
204 if (Xmax <= EpsLix) {
205 if (ErrMax <= EpsLic) {
206 // cout <<"Convergence atteinte dans Uzawa"<<endl;
207 Done = Standard_True;
210 // cout <<"convergence non atteinte pour le probleme dual"<<endl;
211 Done = Standard_False;
214 // Restitution des valeurs calculees
215 //==================================
216 Resul = StartingPoint + Erruza;
217 Done = Standard_True;
222 } // fin de la boucle d iterations.
223 Done = Standard_False;
228 void math_Uzawa::Duale(math_Vector& V) const
233 void math_Uzawa::Dump(Standard_OStream& o) const {
237 o << " Status = Done \n";
238 o << " Number of iterations = " << NbIter << endl;
239 o << " The solution vector is: " << Resul << endl;
242 o << " Status = not Done \n";