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 
MGStep(int NoLevels,QMatrix * A,Vector * x,Vector * b,Matrix * R,Matrix * P,int Level,int Gamma,IterProcType SmoothProc,int Nu1,int Nu2,PrecondProcType PrecondProc,double Omega,IterProcType SolvProc,int NuC,PrecondProcType PrecondProcC,double OmegaC)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 
MGIter(int NoLevels,QMatrix * A,Vector * x,Vector * b,Matrix * R,Matrix * P,int MaxIter,int Gamma,IterProcType SmoothProc,int Nu1,int Nu2,PrecondProcType PrecondProc,double Omega,IterProcType SolvProc,int NuC,PrecondProcType PrecondProcC,double OmegaC)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 
NestedMGIter(int NoLevels,QMatrix * A,Vector * x,Vector * b,Matrix * R,Matrix * P,int Gamma,IterProcType SmoothProc,int Nu1,int Nu2,PrecondProcType PrecondProc,double Omega,IterProcType SolvProc,int NuC,PrecondProcType PrecondProcC,double OmegaC)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 
MGPCGIter(int NoLevels,QMatrix * A,Vector * z,Vector * r,Matrix * R,Matrix * P,int MaxIter,int NoMGIter,int Gamma,IterProcType SmoothProc,int Nu1,int Nu2,PrecondProcType PrecondProc,double Omega,IterProcType SolvProc,int NuC,PrecondProcType PrecondProcC,double OmegaC)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 
BPXPrecond(int NoLevels,QMatrix * A,Vector * y,Vector * c,Matrix * R,Matrix * P,int Level,IterProcType SmoothProc,int Nu,PrecondProcType PrecondProc,double Omega,IterProcType SmoothProcC,int NuC,PrecondProcType PrecondProcC,double OmegaC)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 
BPXPCGIter(int NoLevels,QMatrix * A,Vector * z,Vector * r,Matrix * R,Matrix * P,int MaxIter,IterProcType SmoothProc,int Nu,PrecondProcType PrecondProc,double Omega,IterProcType SmoothProcC,int NuC,PrecondProcType PrecondProcC,double OmegaC)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 }
296