"Fossies" - the Fresh Open Source Software Archive 
Member "laspack/mlsolv.c" (27 Mar 1995, 11413 Bytes) of package /linux/privat/old/laspack.tgz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
1 /****************************************************************************/
2 /* mlsolv.c */
3 /****************************************************************************/
4 /* */
5 /* Multi-Level SOLVers */
6 /* */
7 /* Copyright (C) 1992-1995 Tomas Skalicky. All rights reserved. */
8 /* */
9 /****************************************************************************/
10 /* */
11 /* ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS */
12 /* OF THE COPYRIGHT NOTICE (SEE FILE COPYRGHT.H) */
13 /* */
14 /****************************************************************************/
15
16 #include <math.h>
17 #include <stdio.h>
18 #include <string.h>
19
20 #include "laspack/mlsolv.h"
21 #include "laspack/errhandl.h"
22 #include "laspack/operats.h"
23 #include "laspack/rtc.h"
24 #include "laspack/copyrght.h"
25
26 Vector *MGStep(int NoLevels, QMatrix *A, Vector *x, Vector *b,
27 Matrix *R, Matrix *P, int Level, int Gamma,
28 IterProcType SmoothProc, int Nu1, int Nu2,
29 PrecondProcType PrecondProc, double Omega,
30 IterProcType SolvProc, int NuC,
31 PrecondProcType PrecondProcC, double OmegaC)
32 /* one multigrid iteration */
33 {
34 int CoarseMGIter; /* multi grid iteration counter for coarser grid */
35
36 if (Level == 0) {
37 /* solving of system of equations for the residual on the coarsest grid */
38 (*SolvProc)(&A[Level], &x[Level], &b[Level], NuC, PrecondProcC, OmegaC);
39 } else {
40 /* pre-smoothing - Nu1 iterations */
41 (*SmoothProc)(&A[Level], &x[Level], &b[Level], Nu1, PrecondProc, Omega);
42 /* restiction of the residual to the coarser grid */
43 Asgn_VV(&b[Level - 1], Mul_MV(&R[Level - 1],
44 Sub_VV(&b[Level], Mul_QV(&A[Level], &x[Level]))));
45 /* initialisation of vector of unknowns on the coarser grid */
46 V_SetAllCmp(&x[Level - 1], 0.0);
47 /* solving of system of equations for the residual on the coarser grid */
48 for (CoarseMGIter = 1; CoarseMGIter <= Gamma; CoarseMGIter++)
49 MGStep(NoLevels, A, x, b, R, P, Level - 1, Gamma,
50 SmoothProc, Nu1, Nu2, PrecondProc, Omega,
51 SolvProc, NuC, PrecondProcC, OmegaC);
52 /* interpolation of the solution from the coarser grid */
53 if (P != NULL)
54 AddAsgn_VV(&x[Level], Mul_MV(&P[Level], &x[Level - 1]));
55 else
56 AddAsgn_VV(&x[Level], Mul_MV(Transp_M(&R[Level - 1]), &x[Level - 1]));
57 /* post-smoothing - Nu2 iterations */
58 (*SmoothProc)(&A[Level], &x[Level], &b[Level], Nu2, PrecondProc, Omega);
59 }
60
61 return(&x[Level]);
62 }
63
64 Vector *MGIter(int NoLevels, QMatrix *A, Vector *x, Vector *b,
65 Matrix *R, Matrix *P, int MaxIter, int Gamma,
66 IterProcType SmoothProc, int Nu1, int Nu2,
67 PrecondProcType PrecondProc, double Omega,
68 IterProcType SolvProc, int NuC,
69 PrecondProcType PrecondProcC, double OmegaC)
70 /* multigrid method with residual termination control */
71 {
72 int Iter;
73 double bNorm;
74 size_t Dim;
75 Vector r;
76
77 Dim = Q_GetDim(&A[NoLevels - 1]);
78 V_Constr(&r, "r", Dim, Normal, True);
79
80 if (LASResult() == LASOK) {
81 bNorm = l2Norm_V(&b[NoLevels - 1]);
82
83 Iter = 0;
84 /* r = b - A x(i) at NoLevels - 1 */
85 Asgn_VV(&r, Sub_VV(&b[NoLevels - 1], Mul_QV(&A[NoLevels - 1], &x[NoLevels - 1])));
86 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, MGIterId)
87 && Iter < MaxIter) {
88 Iter++;
89 /* one multigrid step */
90 MGStep(NoLevels, A, x, b, R, P, NoLevels - 1, Gamma,
91 SmoothProc, Nu1, Nu2, PrecondProc, Omega,
92 SolvProc, NuC, PrecondProcC, OmegaC);
93 /* r = b - A x(i) at NoLevels - 1 */
94 Asgn_VV(&r, Sub_VV(&b[NoLevels - 1], Mul_QV(&A[NoLevels - 1], &x[NoLevels - 1])));
95 }
96 }
97
98 V_Destr(&r);
99
100 return(&x[NoLevels - 1]);
101 }
102
103 Vector *NestedMGIter(int NoLevels, QMatrix *A, Vector *x, Vector *b,
104 Matrix *R, Matrix *P, int Gamma,
105 IterProcType SmoothProc, int Nu1, int Nu2,
106 PrecondProcType PrecondProc, double Omega,
107 IterProcType SolvProc, int NuC,
108 PrecondProcType PrecondProcC, double OmegaC)
109 /* nested multigrid method */
110 {
111 int Level;
112
113 /* solution of system of equations on coarsest grid */
114 V_SetAllCmp(&x[0], 0.0);
115 MGStep(NoLevels, A, x, b, R, P, 0, Gamma,
116 SmoothProc, Nu1, Nu2, PrecondProc, Omega,
117 SolvProc, NuC, PrecondProcC, OmegaC);
118
119 for (Level = 1; Level < NoLevels; Level++) {
120 /* prolongation of solution to finer grid */
121 if (P != NULL)
122 Asgn_VV(&x[Level], Mul_MV(&P[Level], &x[Level - 1]));
123 else
124 Asgn_VV(&x[Level], Mul_MV(Transp_M(&R[Level - 1]), &x[Level - 1]));
125 /* solution of system of equations on finer grid with
126 multigrid method */
127 MGStep(NoLevels, A, x, b, R, P, Level, Gamma,
128 SmoothProc, Nu1, Nu2, PrecondProc, Omega,
129 SolvProc, NuC, PrecondProcC, OmegaC);
130 }
131
132 /* submission of reached accuracy to RTC */
133 RTCResult(1, l2Norm_V(Sub_VV(&b[NoLevels - 1],
134 Mul_QV(&A[NoLevels - 1], &x[NoLevels - 1]))),
135 l2Norm_V(&b[NoLevels - 1]), NestedMGIterId);
136
137 return(&x[NoLevels - 1]);
138 }
139
140 Vector *MGPCGIter(int NoLevels, QMatrix *A, Vector *z, Vector *r,
141 Matrix *R, Matrix *P, int MaxIter, int NoMGIter, int Gamma,
142 IterProcType SmoothProc, int Nu1, int Nu2,
143 PrecondProcType PrecondProc, double Omega,
144 IterProcType SolvProc, int NuC,
145 PrecondProcType PrecondProcC, double OmegaC)/* multigrid preconditioned CG method */
146 {
147 int Iter, MGIter;
148 double Alpha, Beta, Rho, RhoOld = 0.0;
149 double bNorm;
150 size_t Dim;
151 Vector x, p, q, b;
152
153 Dim = Q_GetDim(&A[NoLevels - 1]);
154 V_Constr(&x, "x", Dim, Normal, True);
155 V_Constr(&p, "p", Dim, Normal, True);
156 V_Constr(&q, "q", Dim, Normal, True);
157 V_Constr(&b, "b", Dim, Normal, True);
158
159 if (LASResult() == LASOK) {
160 /* copy solution and right hand side stored in parameters z and r */
161 Asgn_VV(&x, &z[NoLevels - 1]);
162 Asgn_VV(&b, &r[NoLevels - 1]);
163
164 bNorm = l2Norm_V(&b);
165
166 Iter = 0;
167 Asgn_VV(&r[NoLevels - 1], Sub_VV(&b, Mul_QV(&A[NoLevels - 1], &x)));
168 while (!RTCResult(Iter, l2Norm_V(&r[NoLevels - 1]), bNorm, MGPCGIterId)
169 && Iter < MaxIter) {
170 Iter++;
171 /* multigrid preconditioner */
172 V_SetAllCmp(&z[NoLevels - 1], 0.0);
173 for (MGIter = 1; MGIter <= NoMGIter; MGIter++)
174 MGStep(NoLevels, A, z, r, R, P, NoLevels - 1, Gamma,
175 SmoothProc, Nu1, Nu2, PrecondProc, Omega,
176 SolvProc, NuC, PrecondProcC, OmegaC);
177 Rho = Mul_VV(&r[NoLevels - 1], &z[NoLevels - 1]);
178 if (Iter == 1) {
179 Asgn_VV(&p, &z[NoLevels - 1]);
180 } else {
181 Beta = Rho / RhoOld;
182 Asgn_VV(&p, Add_VV(&z[NoLevels - 1], Mul_SV(Beta, &p)));
183 }
184 Asgn_VV(&q, Mul_QV(&A[NoLevels - 1], &p));
185 Alpha = Rho / Mul_VV(&p, &q);
186 AddAsgn_VV(&x, Mul_SV(Alpha, &p));
187 SubAsgn_VV(&r[NoLevels - 1], Mul_SV(Alpha, &q));
188 RhoOld = Rho;
189 }
190
191 /* put solution and right hand side vectors back */
192 Asgn_VV(&z[NoLevels - 1], &x);
193 Asgn_VV(&r[NoLevels - 1], &b);
194 }
195
196 V_Destr(&x);
197 V_Destr(&p);
198 V_Destr(&q);
199 V_Destr(&b);
200
201 return(&z[NoLevels - 1]);
202 }
203
204 Vector *BPXPrecond(int NoLevels, QMatrix *A, Vector *y, Vector *c,
205 Matrix *R, Matrix *P, int Level,
206 IterProcType SmoothProc, int Nu,
207 PrecondProcType PrecondProc, double Omega,
208 IterProcType SmoothProcC, int NuC,
209 PrecondProcType PrecondProcC, double OmegaC)
210 /* BPX preconditioner (recursively defined) */
211 {
212 if (Level == 0) {
213 /* smoothing on the coarsest grid - NuC iterations */
214 V_SetAllCmp(&y[Level], 0.0);
215 (*SmoothProcC)(&A[Level], &y[Level], &c[Level], NuC, PrecondProcC, OmegaC);
216 } else {
217 /* smoothing - Nu iterations */
218 V_SetAllCmp(&y[Level], 0.0);
219 (*SmoothProc)(&A[Level], &y[Level], &c[Level], Nu, PrecondProc, Omega);
220 /* restiction of the residual to the coarser grid */
221 Asgn_VV(&c[Level - 1], Mul_MV(&R[Level - 1], &c[Level]));
222 /* smoothing on the coarser grid */
223 BPXPrecond(NoLevels, A, y, c, R, P, Level - 1,
224 SmoothProc, Nu, PrecondProc, Omega, SmoothProcC, NuC, PrecondProcC, OmegaC);
225 /* interpolation of the solution from coarser grid */
226 if (P != NULL)
227 AddAsgn_VV(&y[Level], Mul_MV(&P[Level], &y[Level - 1]));
228 else
229 AddAsgn_VV(&y[Level], Mul_MV(Transp_M(&R[Level - 1]), &y[Level - 1]));
230 }
231
232 return(&y[Level]);
233 }
234
235 Vector *BPXPCGIter(int NoLevels, QMatrix *A, Vector *z, Vector *r,
236 Matrix *R, Matrix *P, int MaxIter,
237 IterProcType SmoothProc, int Nu,
238 PrecondProcType PrecondProc, double Omega,
239 IterProcType SmoothProcC, int NuC,
240 PrecondProcType PrecondProcC, double OmegaC)
241 /* BPX preconditioned CG method */
242 {
243 int Iter;
244 double Alpha, Beta, Rho, RhoOld = 0.0;
245 double bNorm;
246 size_t Dim;
247 Vector x, p, q, b;
248
249 Dim = Q_GetDim(&A[NoLevels - 1]);
250 V_Constr(&x, "x", Dim, Normal, True);
251 V_Constr(&p, "p", Dim, Normal, True);
252 V_Constr(&q, "q", Dim, Normal, True);
253 V_Constr(&b, "b", Dim, Normal, True);
254
255 if (LASResult() == LASOK) {
256 /* copy solution and right hand side stored in parameters z and r */
257 Asgn_VV(&x, &z[NoLevels - 1]);
258 Asgn_VV(&b, &r[NoLevels - 1]);
259
260 bNorm = l2Norm_V(&b);
261
262 Iter = 0;
263 Asgn_VV(&r[NoLevels - 1], Sub_VV(&b, Mul_QV(&A[NoLevels - 1], &x)));
264 while (!RTCResult(Iter, l2Norm_V(&r[NoLevels - 1]), bNorm, BPXPCGIterId)
265 && Iter < MaxIter) {
266 Iter++;
267 /* BPX preconditioner */
268 BPXPrecond(NoLevels, A, z, r, R, P, NoLevels - 1,
269 SmoothProc, Nu, PrecondProc, Omega, SmoothProcC, NuC, PrecondProcC, OmegaC);
270 Rho = Mul_VV(&r[NoLevels - 1], &z[NoLevels - 1]);
271 if (Iter == 1) {
272 Asgn_VV(&p, &z[NoLevels - 1]);
273 } else {
274 Beta = Rho / RhoOld;
275 Asgn_VV(&p, Add_VV(&z[NoLevels - 1], Mul_SV(Beta, &p)));
276 }
277 Asgn_VV(&q, Mul_QV(&A[NoLevels - 1], &p));
278 Alpha = Rho / Mul_VV(&p, &q);
279 AddAsgn_VV(&x, Mul_SV(Alpha, &p));
280 SubAsgn_VV(&r[NoLevels - 1], Mul_SV(Alpha, &q));
281 RhoOld = Rho;
282 }
283
284 /* put solution and right hand side vectors back */
285 Asgn_VV(&z[NoLevels - 1], &x);
286 Asgn_VV(&r[NoLevels - 1], &b);
287 }
288
289 V_Destr(&x);
290 V_Destr(&p);
291 V_Destr(&q);
292 V_Destr(&b);
293
294 return(&z[NoLevels - 1]);
295 }