38291ba43584a419bcdc1232814065b5e7a1d091
[occt.git] / src / math / math_Uzawa.cxx
1 //static const char* sccsid = "@(#)math_Uzawa.cxx       3.2 95/01/10"; // Do not delete this line. Used by sccs.
2 // File math_Uzawa.cxx
3 // lpa le 8/08/91
4
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.
13
14 //#ifndef DEB
15 #define No_Standard_RangeError
16 #define No_Standard_OutOfRange
17 #define No_Standard_DimensionError
18 //#endif
19
20 #include <math_Uzawa.ixx>
21 #include <math_Crout.hxx>
22 #include <Standard_DimensionError.hxx>
23
24
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(),
34                               1, Cont.RowNumber()){
35
36  Perform(Cont, Secont, StartingPoint, Cont.RowNumber(), 0, EpsLix,
37          EpsLic, NbIterations);
38 }
39
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())  {
51
52   Perform(Cont, Secont, StartingPoint, Nce, Nci, EpsLix, EpsLic, NbIterations);
53
54 }
55
56
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)  {
62
63   Standard_Integer i,j, k;
64   Standard_Real scale;
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();
69
70   Standard_DimensionError_Raise_if((Secont.Length() != Nlig) ||
71                                    ((Nce+Nci) != Nlig), " ");
72
73   // Calcul du vecteur Cont*X0 - D:  (erreur initiale)
74   //==================================================
75
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);
80       }
81     }
82
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);
92            }
93          }
94        }
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);
104       }
105       for (j =i+1; j <= Nlig; j++) {
106         scale += CTCinv(j,i)*Errinit(j);
107       }
108       Vardua(i) = scale;
109     }
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);
114       }
115     }
116     // restitution des valeurs calculees:
117     //===================================
118     Resul = StartingPoint + Erruza;
119     Done = Standard_True;
120     return;
121   } // Fin de la resolution directe.
122     //==============================
123
124   else {  
125     // Initialisation des variables duales.
126     //=====================================
127     for (i = 1; i <= Nlig; i++) {
128       if (i <= Nce) {
129         Vardua(i) = 0.0;
130       }
131       else {
132         Vardua(i) = 1.;
133       }
134     }
135     
136     // Calcul du coefficient Rho:
137     //===========================
138     Normat = 0.0;
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);
143       }
144       Normat += Normli;
145     }
146     Rho = Coef/Normat;
147     
148     // Boucle des iterations de la methode d Uzawa.
149     //=============================================
150     for (NbIter = 1; NbIter <= NbIterations; NbIter++) {
151       for (i = 1; i <= Ncol; i++) {
152         Xian = Erruza(i);
153         Erruza(i) = -Cont(1,i)*Vardua(1);
154         for(j =2; j <= Nlig; j++) {
155           Erruza(i) -= Cont(j,i)*Vardua(j);
156         }
157         if (NbIter > 1) {                      
158           if (i == 1) {
159             Xmax = Abs(Erruza(i) - Xian);
160           }
161           Xmax = Max(Xmax, Abs(Erruza(i)-Xian));
162         }
163       }
164
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);
172         }
173         if (i <= Nce) {
174           Vardua(i) += Rho * Err;
175           Err1 = Abs(Rho*Err);
176         }
177         else {
178           Xmuian = Vardua(i);
179           Vardua(i) = Max(0.0, Vardua(i)+ Rho*Err);
180           Err1 = Abs(Vardua(i) - Xmuian);
181         }
182         if (i == 1) {
183           ErrMax = Err1;
184         }
185         ErrMax = Max(ErrMax, Err1);
186       }
187
188       if (NbIter > 1) {                   
189         if (Xmax <= EpsLix) {
190           if (ErrMax <= EpsLic) {
191 //          cout <<"Convergence atteinte dans Uzawa"<<endl;
192             Done = Standard_True;
193           }
194           else {
195 //          cout <<"convergence non atteinte pour le probleme dual"<<endl;
196             Done = Standard_False;
197             return;
198           }
199           // Restitution des valeurs calculees
200           //==================================
201           Resul = StartingPoint + Erruza;
202           Done = Standard_True;
203           return;
204         }
205       }
206
207     } // fin de la boucle d iterations.
208       Done = Standard_False;
209   }
210 }
211
212
213 void math_Uzawa::Duale(math_Vector& V) const
214 {
215   V = Vardua;
216 }
217
218 void math_Uzawa::Dump(Standard_OStream& o) const {
219
220   o << "math_Uzawa";
221   if(Done) {
222     o << " Status = Done \n";
223     o << " Number of iterations = " << NbIter << endl;
224     o << " The solution vector is: " << Resul << endl;
225   }
226   else {
227     o << " Status = not Done \n";
228   }
229 }
230
231
232
233
234