1 /****************************************************************************/
2 /* mlstest.c */
3 /****************************************************************************/
4 /* */
5 /* Multi-Level Solver TEST program for 1d and 2d poisson problems */
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 <stdio.h>
17 #include <math.h>
18 #include <ctype.h>
19 #include <time.h>
20 #include <float.h>
21
22 #include <laspack/errhandl.h>
23 #include <laspack/vector.h>
24 #include <laspack/matrix.h>
25 #include <laspack/qmatrix.h>
26 #include <laspack/operats.h>
27 #include <laspack/factor.h>
28 #include <laspack/precond.h>
29 #include <laspack/eigenval.h>
30 #include <laspack/rtc.h>
31 #include <laspack/itersolv.h>
32 #include <laspack/mlsolv.h>
33 #include <laspack/version.h>
34 #include <laspack/copyrght.h>
35
36 #include "mlstest.h"
37 #include "extrsolv.h"
38
39 typedef enum {
40 Case1DDirich,
41 Case2DDirich
42 } ProblemType;
43
44 static QMatrix *L;
45 static Vector *uw, *fr;
46 static Matrix *R, *P;
47
48 static void MGParmInput(ProblemType *Problem, int *NoLevels, size_t *MaxNoInt,
49 IterIdType *MLSolverId, Boolean *NestedMG, int *NoMGIter,
50 int *Gamma, int *RestrType,
51 IterProcType *SmoothProc, int *Nu1, int *Nu2,
52 PrecondProcType *PrecondProc, double *Omea,
53 IterProcType *SolvProc, int *NuC,
54 PrecondProcType *PrecondProcC, double *OmegaC,
55 int *MaxIter, double *Eps);
56 static void QVMConstr(int NoLevels, size_t *Dim);
57 static void QVMDestr(int NoLevels);
58 static void MGEqsGen1DDirich(int NoLevels, size_t *Dim, int RestrType);
59 static void MGEqsGen2DDirich(int NoLevels, size_t *Dim1, int RestrType);
60 static void MGSolver(int NoLevels, IterIdType MLSolverId,
61 Boolean NestedMG, int NoMGIter, int Gamma,
62 IterProcType SmoothProc, int Nu1, int Nu2,
63 PrecondProcType PrecondProc, double Omega,
64 IterProcType SolvProc, int NuC,
65 PrecondProcType PrecondProcC, double OmegaC,
66 int MaxIter, double Eps);
67 static void MGResOutput1DDirich(int NoLevels);
68 static void GenL1DDirich(int Level, size_t *Dim);
69 static void Genf1DDirich(int Level, size_t *Dim);
70 static void GenRSimple1DDirich(int Level, size_t *Dim);
71 static void GenRWeight1DDirich(int Level, size_t *Dim);
72 static void GenP1DDirich(int Level, size_t *Dim);
73 static void GenL2DDirich(int Level, size_t *Dim1);
74 static void Genf2DDirich(int Level, size_t *Dim1);
75 static void GenRSimple2DDirich(int Level, size_t *Dim1);
76 static void GenRWeight2DDirich(int Level, size_t *Dim1);
77 static void GenP2DDirich(int Level, size_t *Dim1);
78 static void IterStatus(int Iter, double rNorm, double bNorm, IterIdType IterId);
79
main(void)80 int main(void)
81 {
82 int NoLevels; /* number of grid levels */
83 int RestrType; /* type of restriction operator */
84 int Gamma; /* multigrid iteration parameter */
85 int Nu1; /* number of pre-smoothing iterations */
86 int Nu2; /* number of post-smoothing iterations */
87 int NuC; /* number of iterations on coarse grid */
88 int MaxIter; /* maximum number of iterations */
89 int NoMGIter; /* number of multigrid iterations within one MGPCG step */
90 int Level;
91 double Omega; /* relaxation parameter for smoothing iterations */
92 double OmegaC; /* relaxation parameter for iterations on coarse grid */
93 double Eps; /* epsilon - break off accurary */
94 size_t MaxNoInt; /* maximal number of intervals */
95 size_t *Dim; /* dimension, point number */
96 size_t *Dim1; /* dimension, point number in one direction */
97 Boolean NestedMG; /* nested multigrid iteration */
98 IterIdType MLSolverId; /* identifier of the chosen multi-level solver */
99 IterProcType SmoothProc; /* pointer to a procedure for smoothing iterations */
100 IterProcType SolvProc; /* pointer to a procedure for solving on coarsest grid */
101 PrecondProcType PrecondProc; /* pointer to preconditioner for smoothing iterations */
102 PrecondProcType PrecondProcC; /* pointer to preconditioner on coarsest grid */
103 ProblemType Problem; /* kind of problem */
104
105 fprintf(stderr, "mlstest Version %s\n", LASPACK_VERSION);
106 fprintf(stderr, " (C) 1992-1995 Tomas Skalicky\n");
107 fprintf(stderr,"\n\n");
108 fprintf(stderr,"Multilevel solution of a Poisson problem\n");
109 fprintf(stderr,"----------------------------------------\n");
110 fprintf(stderr,"\n");
111 /* multigrid paramater input */
112 MGParmInput(&Problem, &NoLevels, &MaxNoInt,
113 &MLSolverId, &NestedMG, &NoMGIter, &Gamma, &RestrType,
114 &SmoothProc, &Nu1, &Nu2, &PrecondProc, &Omega,
115 &SolvProc, &NuC, &PrecondProcC, &OmegaC, &MaxIter, &Eps);
116
117 /* allocation of dynamic variables */
118 Dim = (size_t *)malloc(NoLevels * sizeof(size_t));
119 Dim1 = (size_t *)malloc(NoLevels * sizeof(size_t));
120 L = (QMatrix *)malloc(NoLevels * sizeof(QMatrix));
121 uw = (Vector *)malloc(NoLevels * sizeof(Vector));
122 fr = (Vector *)malloc(NoLevels * sizeof(Vector));
123 R = (Matrix *)malloc(NoLevels * sizeof(Matrix));
124 P = (Matrix *)malloc(NoLevels * sizeof(Matrix));
125 if (Dim != NULL && Dim1 != NULL && L != NULL && uw != NULL && fr != NULL &&
126 R != NULL && P != NULL) {
127 /* setting of dimensions for all grid levels */
128 Dim1[NoLevels - 1] = MaxNoInt - 1;
129 for (Level = NoLevels - 2; Level >= 0; Level--)
130 Dim1[Level] = (Dim1[Level + 1] + 1) / 2 - 1;
131 switch (Problem) {
132 case Case1DDirich:
133 for (Level = NoLevels - 1; Level >= 0; Level--)
134 Dim[Level] = Dim1[Level];
135 break;
136 case Case2DDirich:
137 for (Level = NoLevels - 1; Level >= 0; Level--)
138 Dim[Level] = Dim1[Level] * Dim1[Level];
139 break;
140 }
141 /* allocation of matrices, vectors and operators */
142 QVMConstr(NoLevels, Dim);
143 /* generation of system of equations */
144 switch (Problem) {
145 case Case1DDirich:
146 MGEqsGen1DDirich(NoLevels, Dim, RestrType);
147 break;
148 case Case2DDirich:
149 MGEqsGen2DDirich(NoLevels, Dim1, RestrType);
150 break;
151 }
152 /* solving of system of equations with multigrid method */
153 MGSolver(NoLevels, MLSolverId, NestedMG, NoMGIter, Gamma,
154 SmoothProc, Nu1, Nu2, PrecondProc, Omega,
155 SolvProc, NuC, PrecondProcC, OmegaC, MaxIter, Eps);
156 /* solution output */
157 if (LASResult() == LASOK && Problem == Case1DDirich)
158 MGResOutput1DDirich(NoLevels);
159 /* release of matrices, vectors and operators */
160 QVMDestr(NoLevels);
161 } else {
162 /* error message */
163 fprintf(stderr, "\n");
164 fprintf(stderr, "Not enought memory running mgtest.\n");
165 }
166
167 /* release of dynamic variables */
168 if (Dim != NULL)
169 free(Dim);
170 if (Dim1 != NULL)
171 free(Dim1);
172 if (L != NULL)
173 free(L);
174 if (uw != NULL)
175 free(uw);
176 if (fr != NULL)
177 free(fr);
178 if (R != NULL)
179 free(R);
180 if (P != NULL)
181 free(P);
182
183 /* LASPack error messages */
184 if (LASResult() != LASOK) {
185 fprintf(stderr, "\n");
186 fprintf(stderr, "LASPack error: ");
187 WriteLASErrDescr(stderr);
188 }
189 fprintf(stderr, "\n");
190
191 return(0);
192 }
193
MGParmInput(ProblemType * Problem,int * NoLevels,size_t * MaxNoInt,IterIdType * MLSolverId,Boolean * NestedMG,int * NoMGIter,int * Gamma,int * RestrType,IterProcType * SmoothProc,int * Nu1,int * Nu2,PrecondProcType * PrecondProc,double * Omega,IterProcType * SolvProc,int * NuC,PrecondProcType * PrecondProcC,double * OmegaC,int * MaxIter,double * Eps)194 static void MGParmInput(ProblemType *Problem, int *NoLevels, size_t *MaxNoInt,
195 IterIdType *MLSolverId, Boolean *NestedMG, int *NoMGIter,
196 int *Gamma, int *RestrType,
197 IterProcType *SmoothProc, int *Nu1, int *Nu2,
198 PrecondProcType *PrecondProc, double *Omega,
199 IterProcType *SolvProc, int *NuC,
200 PrecondProcType *PrecondProcC, double *OmegaC,
201 int *MaxIter, double *Eps)
202 /* multigrid parameter input */
203 {
204 char Key[2];
205 int Level;
206 int ProblemId; /* problem identifier */
207 int MainMethId; /* main method identifier */
208 int SmoothMethId; /* smoothing method identifier */
209 int SolvMethId; /* solving method on coarse grid identifier */
210 int PrecondId; /* preconditioning identifier */
211 unsigned long AuxSize;
212 size_t CoarseCoef; /* coarsing coefficient */
213 Boolean InputOK;
214
215 fprintf(stderr, "\n");
216 fprintf(stderr, "Problem type: 1D, Dirichlet boundary condition ... 1\n");
217 fprintf(stderr, " 2D, Dirichlet boundary condition ... 2\n");
218 fprintf(stderr, "Chosen type: ");
219 do {
220 scanf("%d", &ProblemId);
221 if (ProblemId < 1 || ProblemId > 2)
222 fprintf(stderr, "??? : ");
223 } while (ProblemId < 1 || ProblemId > 2);
224 switch (ProblemId) {
225 case 1:
226 *Problem = Case1DDirich;
227 break;
228 case 2:
229 *Problem = Case2DDirich;
230 break;
231 }
232 do {
233 fprintf(stderr, "\n");
234 fprintf(stderr, "Number of grid levels: ");
235 scanf("%d", NoLevels);
236 CoarseCoef = 1;
237 for (Level = *NoLevels - 1; Level > 0; Level--)
238 CoarseCoef = 2 * CoarseCoef;
239 fprintf(stderr, "Number of intervals: ");
240 scanf("%lu", &AuxSize);
241 *MaxNoInt = (size_t)AuxSize;
242 InputOK = (*MaxNoInt % CoarseCoef == 0 && *MaxNoInt / CoarseCoef > 2);
243 if (!InputOK) {
244 fprintf(stderr, "Discretisation by this choice is impossible!\n");
245 }
246 } while (!InputOK);
247
248 fprintf(stderr, "\n");
249 fprintf(stderr, "Solution method: multigrid ..................... 1\n");
250 fprintf(stderr, " multigrid preconditioned CG ... 2\n");
251 fprintf(stderr, " BPX preconditioned CG ......... 3\n");
252 fprintf(stderr, "Chosen method: ");
253 do {
254 scanf("%d", &MainMethId);
255 if (MainMethId < 1 || MainMethId > 3)
256 fprintf(stderr, "??? : ");
257 } while (MainMethId < 1 || MainMethId > 3);
258 switch (MainMethId) {
259 case 1:
260 *MLSolverId = MGIterId;
261 fprintf(stderr, "\n");
262 do {
263 fprintf(stderr, "Nested multigrid iterations? (y/n) ");
264 scanf("%1s", Key);
265 Key[0] = toupper(Key[0]);
266 } while (Key[0] != 'Y' && Key[0] != 'N');
267 if (Key[0] == 'Y')
268 *NestedMG = True;
269 else
270 *NestedMG = False;
271 break;
272 case 2:
273 *MLSolverId = MGPCGIterId;
274 fprintf(stderr, "\n");
275 fprintf(stderr, "Number of multigrid iterations whithin one MGPCG step: ");
276 scanf("%d", NoMGIter);
277 break;
278 case 3:
279 *MLSolverId = BPXPCGIterId;
280 break;
281 }
282 if (*MLSolverId == MGIterId || *MLSolverId == MGPCGIterId) {
283 fprintf(stderr, "\n");
284 fprintf(stderr, "Type of multigrid: V cycle ... 1\n");
285 fprintf(stderr, " W cycle ... 2\n");
286 fprintf(stderr, "Chosen type: ");
287 scanf("%d", Gamma);
288 }
289 fprintf(stderr, "\n");
290 fprintf(stderr, "Type of restriction: simple ..... 1\n");
291 fprintf(stderr, " weighted ... 2\n");
292 fprintf(stderr, "Chosen type: ");
293 do {
294 scanf("%d", RestrType);
295 if (*RestrType != 1 && *RestrType != 2)
296 fprintf(stderr, "??? : ");
297 } while (*RestrType != 1 && *RestrType != 2);
298
299 fprintf(stderr, "\n");
300 fprintf(stderr, "Iterative methods: Jacobi ................... 1\n");
301 fprintf(stderr, " SOR forward .............. 2\n");
302 fprintf(stderr, " SOR backward ............. 3\n");
303 fprintf(stderr, " SOR symmetric ............ 4\n");
304 fprintf(stderr, " Chebyshev ................ 5\n");
305 fprintf(stderr, " CG ....................... 6\n");
306 fprintf(stderr, " CGN ...................... 7\n");
307 fprintf(stderr, " GMRES(10) ................ 8\n");
308 fprintf(stderr, " BiCG ..................... 9\n");
309 fprintf(stderr, " QMR ..................... 10\n");
310 fprintf(stderr, " CGS ..................... 11\n");
311 fprintf(stderr, " Bi-CGSTAB ............... 12\n");
312 fprintf(stderr, " Test .................... 13\n");
313 fprintf(stderr, "Smoothing method: ");
314 do {
315 scanf("%d", &SmoothMethId);
316 if (SmoothMethId < 1 || SmoothMethId > 13)
317 fprintf(stderr, "??? : ");
318 } while (SmoothMethId < 1 || SmoothMethId > 13);
319 switch (SmoothMethId) {
320 case 1:
321 *SmoothProc = JacobiIter;
322 break;
323 case 2:
324 *SmoothProc = SORForwIter;
325 break;
326 case 3:
327 *SmoothProc = SORBackwIter;
328 break;
329 case 4:
330 *SmoothProc = SSORIter;
331 break;
332 case 5:
333 *SmoothProc = ChebyshevIter;
334 break;
335 case 6:
336 *SmoothProc = CGIter;
337 break;
338 case 7:
339 *SmoothProc = CGNIter;
340 break;
341 case 8:
342 *SmoothProc = GMRESIter;
343 break;
344 case 9:
345 *SmoothProc = BiCGIter;
346 break;
347 case 10:
348 *SmoothProc = QMRIter;
349 break;
350 case 11:
351 *SmoothProc = CGSIter;
352 break;
353 case 12:
354 *SmoothProc = BiCGSTABIter;
355 break;
356 case 13:
357 *SmoothProc = TestIter;
358 break;
359 }
360 if (*MLSolverId == MGIterId || *MLSolverId == MGPCGIterId) {
361 fprintf(stderr, "Number of pre-smoothing iterations: ");
362 scanf("%d", Nu1);
363 fprintf(stderr, "Number of post-smoothing iterations: ");
364 scanf("%d", Nu2);
365 } else {
366 fprintf(stderr, "Number of smoothing iterations: ");
367 scanf("%d", Nu1);
368 Nu2 = 0;
369 }
370 fprintf(stderr, "Preconditioning: none .......... 0\n");
371 fprintf(stderr, " Jacobi ........ 1\n");
372 fprintf(stderr, " SSOR ......... 2\n");
373 fprintf(stderr, " ILU/ICH ....... 3\n");
374 fprintf(stderr, "Preconditioning for smoothing iterations: ");
375 do {
376 scanf("%d", &PrecondId);
377 if (PrecondId < 0 || PrecondId > 3)
378 fprintf(stderr, "??? : ");
379 } while (PrecondId < 0 || PrecondId > 3);
380 switch (PrecondId) {
381 case 0:
382 *PrecondProc = (PrecondProcType)NULL;
383 break;
384 case 1:
385 *PrecondProc = JacobiPrecond;
386 break;
387 case 2:
388 *PrecondProc = SSORPrecond;
389 break;
390 case 3:
391 *PrecondProc = ILUPrecond;
392 break;
393 }
394 fprintf(stderr, "Relaxation parameter for smoothing: ");
395 scanf("%lf", Omega);
396
397 fprintf(stderr, "Iterative methods: Jacobi ................... 1\n");
398 fprintf(stderr, " SOR forward .............. 2\n");
399 fprintf(stderr, " SOR backward ............. 3\n");
400 fprintf(stderr, " SOR symmetric ............ 4\n");
401 fprintf(stderr, " Chebyshev ................ 5\n");
402 fprintf(stderr, " CG ....................... 6\n");
403 fprintf(stderr, " CGN ...................... 7\n");
404 fprintf(stderr, " GMRES(10) ................ 8\n");
405 fprintf(stderr, " BiCG ..................... 9\n");
406 fprintf(stderr, " QMR ..................... 10\n");
407 fprintf(stderr, " CGS ..................... 11\n");
408 fprintf(stderr, " Bi-CGSTAB ............... 12\n");
409 fprintf(stderr, " Test .................... 13\n");
410 fprintf(stderr, "Solution method on coarsest grid: ");
411 do {
412 scanf("%d", &SolvMethId);
413 if (SolvMethId < 1 || SolvMethId > 13)
414 fprintf(stderr, "??? : ");
415 } while (SolvMethId < 1 || SolvMethId > 13);
416 switch (SolvMethId) {
417 case 1:
418 *SolvProc = JacobiIter;
419 break;
420 case 2:
421 *SolvProc = SORForwIter;
422 break;
423 case 3:
424 *SolvProc = SORBackwIter;
425 break;
426 case 4:
427 *SolvProc = SSORIter;
428 break;
429 case 5:
430 *SolvProc = ChebyshevIter;
431 break;
432 case 6:
433 *SolvProc = CGIter;
434 break;
435 case 7:
436 *SolvProc = CGNIter;
437 break;
438 case 8:
439 *SolvProc = GMRESIter;
440 break;
441 case 9:
442 *SolvProc = BiCGIter;
443 break;
444 case 10:
445 *SolvProc = QMRIter;
446 break;
447 case 11:
448 *SolvProc = CGSIter;
449 break;
450 case 12:
451 *SolvProc = BiCGSTABIter;
452 break;
453 case 13:
454 *SolvProc = TestIter;
455 break;
456 }
457 fprintf(stderr, "Number of iterations on coarsest grid: ");
458 scanf("%d", NuC);
459 fprintf(stderr, "Preconditioning: none .......... 0\n");
460 fprintf(stderr, " Jacobi ........ 1\n");
461 fprintf(stderr, " SSOR ......... 2\n");
462 fprintf(stderr, " ILU/ICH ....... 3\n");
463 fprintf(stderr, "Preconditioning on coarsest grid: ");
464 do {
465 scanf("%d", &PrecondId);
466 if (PrecondId < 0 || PrecondId > 3)
467 fprintf(stderr, "??? : ");
468 } while (PrecondId < 0 || PrecondId > 3);
469 switch (PrecondId) {
470 case 0:
471 *PrecondProcC = (PrecondProcType)NULL;
472 break;
473 case 1:
474 *PrecondProcC = JacobiPrecond;
475 break;
476 case 2:
477 *PrecondProcC = SSORPrecond;
478 break;
479 case 3:
480 *PrecondProcC = ILUPrecond;
481 break;
482 }
483 fprintf(stderr, "Relaxation parameter on coarsest grid: ");
484 scanf("%lf", OmegaC);
485
486 fprintf(stderr, "\n");
487 fprintf(stderr, "Maximum number of iterations: ");
488 scanf("%d", MaxIter);
489 fprintf(stderr, "Break off accurary for the residual: ");
490 scanf("%lf", Eps);
491 }
492
QVMConstr(int NoLevels,size_t * Dim)493 static void QVMConstr(int NoLevels, size_t *Dim)
494 /* calling of constructors for matrices, vectors of unknowns,
495 vectors of right hand side, restriction and prolongation operators */
496 {
497 int Level;
498
499 char Name[10];
500
501 for (Level = NoLevels - 1; Level >= 0; Level--) {
502 sprintf(Name, "L[%d]", Level % 1000);
503 #ifdef SYM_STOR
504 Q_Constr(&L[Level], Name, Dim[Level], True, Rowws, Normal, True);
505 #else
506 Q_Constr(&L[Level], Name, Dim[Level], False, Rowws, Normal, True);
507 #endif /* SYM_STOR */
508 sprintf(Name, "uw[%d]", Level % 1000);
509 V_Constr(&uw[Level], Name, Dim[Level], Normal, True);
510 sprintf(Name, "fr[%d]", Level % 1000);
511 V_Constr(&fr[Level], Name, Dim[Level], Normal, True);
512 if (Level < NoLevels - 1) {
513 sprintf(Name, "R[%d]", Level % 1000);
514 M_Constr(&R[Level], Name, Dim[Level], Dim[Level + 1], Rowws, Normal, True);
515 }
516 if (Level > 0) {
517 sprintf(Name, "P[%d]", Level % 1000);
518 M_Constr(&P[Level], Name, Dim[Level], Dim[Level - 1], Clmws, Normal, True);
519 }
520 }
521 }
522
QVMDestr(int NoLevels)523 static void QVMDestr(int NoLevels)
524 /* calling of destructors for matrices, vectors of unknowns,
525 vectors of right hand side, restriction and prolongation operators */
526 {
527 int Level;
528 for (Level = NoLevels - 1; Level >= 0; Level--) {
529 Q_Destr(&L[Level]);
530 V_Destr(&uw[Level]);
531 V_Destr(&fr[Level]);
532 if (Level < NoLevels - 1)
533 M_Destr(&R[Level]);
534 if (Level > 0)
535 M_Destr(&P[Level]);
536 }
537 }
538
MGEqsGen1DDirich(int NoLevels,size_t * Dim,int RestrType)539 static void MGEqsGen1DDirich(int NoLevels, size_t *Dim, int RestrType)
540 /* generation of system of equations for all grid levels for 1D case,
541 Dirichlet boundary condition */
542 {
543 int Level;
544 double ClockFactor, CPUTime;
545 clock_t BegClock, EndClock;
546
547 /* start of time counting */
548 ClockFactor = 1.0 / (double)CLOCKS_PER_SEC;
549 BegClock = clock();
550 /* generation of matrices L */
551 printf("\n");
552 printf("Generation of matrices ... \n");
553 for (Level = NoLevels - 1; Level >= 0; Level--)
554 GenL1DDirich(Level, Dim);
555 /* generation of right hand side f */
556 printf("Generation of vectors ... \n");
557 for (Level = NoLevels - 1; Level >= 0; Level--)
558 Genf1DDirich(Level, Dim);
559 /* generation of restriction operators */
560 printf("Generation of restriction operators ... \n");
561 for (Level = NoLevels - 2; Level >= 0; Level--) {
562 if (RestrType == 1) {
563 /* for simple restriction operators */
564 GenRSimple1DDirich(Level, Dim);
565 }
566 if (RestrType == 2) {
567 /* for wighted restriction operators */
568 GenRWeight1DDirich(Level, Dim);
569 }
570 }
571 printf("Generation of prolongation operators ... \n");
572 for (Level = NoLevels - 1; Level >= 1; Level--)
573 GenP1DDirich(Level, Dim);
574 /* end of time counting and geting out of CPU time */
575 EndClock = clock();
576 CPUTime = (double)(EndClock - BegClock) * ClockFactor;
577 printf("\n");
578 printf(" CPU time: %7.2f s\n", CPUTime);
579 }
580
MGEqsGen2DDirich(int NoLevels,size_t * Dim1,int RestrType)581 static void MGEqsGen2DDirich(int NoLevels, size_t *Dim1, int RestrType)
582 /* generation of system of equations for all grid levels for 2D case,
583 Dirichlet boundary condition */
584 {
585 int Level;
586 double ClockFactor, CPUTime;
587 clock_t BegClock, EndClock;
588
589 /* start of time counting */
590 ClockFactor = 1.0 / (double)CLOCKS_PER_SEC;
591 BegClock = clock();
592 /* generation of matrices L */
593 printf("\n");
594 printf("Generation of matrices ... \n");
595 for (Level = NoLevels - 1; Level >= 0; Level--)
596 GenL2DDirich(Level, Dim1);
597 /* generation of right hand side f */
598 printf("Generation of vectors ... \n");
599 for (Level = NoLevels - 1; Level >= 0; Level--)
600 Genf2DDirich(Level, Dim1);
601 /* generation of restriction operators */
602 printf("Generation of restriction operators ... \n");
603 for (Level = NoLevels - 2; Level >= 0; Level--) {
604 if (RestrType == 1) {
605 /* for simple restriction operators */
606 GenRSimple2DDirich(Level, Dim1);
607 }
608 if (RestrType == 2) {
609 /* for wighted restriction operators */
610 GenRWeight2DDirich(Level, Dim1);
611 }
612 }
613 /* generation of prolongation operators P */
614 printf("Generation of prolongation operators ... \n");
615 for (Level = NoLevels - 1; Level >= 1; Level--)
616 GenP2DDirich(Level, Dim1);
617 /* end of time counting and geting out of CPU time */
618 EndClock = clock();
619 CPUTime = (double)(EndClock - BegClock) * ClockFactor;
620 printf("\n");
621 printf(" CPU time: %7.2f s\n", CPUTime);
622 }
623
MGSolver(int NoLevels,IterIdType MLSolverId,Boolean NestedMG,int NoMGIter,int Gamma,IterProcType SmoothProc,int Nu1,int Nu2,PrecondProcType PrecondProc,double Omega,IterProcType SolvProc,int NuC,PrecondProcType PrecondProcC,double OmegaC,int MaxIter,double Eps)624 static void MGSolver(int NoLevels, IterIdType MLSolverId,
625 Boolean NestedMG, int NoMGIter, int Gamma,
626 IterProcType SmoothProc, int Nu1, int Nu2,
627 PrecondProcType PrecondProc, double Omega,
628 IterProcType SolvProc, int NuC,
629 PrecondProcType PrecondProcC, double OmegaC,
630 int MaxIter, double Eps)
631 /* solution of discret problem with multigrid solver */
632 {
633 int NoIter; /* number of performed iterations */
634 int Level;
635 double AccBeg, AccEnd; /* reached accuracy (of residuum) */
636 double ContrRatePerMGIter, ContrRatePerSec;
637 double ClockFactor, CPUTime;
638 clock_t BegClock, EndClock;
639
640 /* setting of RTC parameters */
641 SetRTCAccuracy(Eps);
642 SetRTCAuxProc(IterStatus);
643
644 /* estimate extremal eigenvalues */
645 if(SmoothProc == ChebyshevIter || SolvProc == ChebyshevIter) {
646 /* start of time counting */
647 ClockFactor = 1.0 / (double)CLOCKS_PER_SEC;
648 BegClock = clock();
649
650 /* initialization random-number generator */
651 srand(1);
652
653 printf("\n");
654 printf("Estimating extremal eigenvalues ...\n");
655 if(SmoothProc == ChebyshevIter) {
656 for (Level = NoLevels - 1; Level > 0; Level--) {
657 GetMinEigenval(&L[Level], PrecondProc, Omega);
658 GetMaxEigenval(&L[Level], PrecondProc, Omega);
659 }
660 }
661 if(SolvProc == ChebyshevIter) {
662 GetMinEigenval(&L[0], PrecondProcC, OmegaC);
663 GetMaxEigenval(&L[0], PrecondProcC, OmegaC);
664 }
665
666 /* end of time counting and geting out the CPU time */
667 EndClock = clock();
668 CPUTime = (double)(EndClock - BegClock) * ClockFactor;
669 printf("\n");
670 printf(" CPU time: %7.2f s\n", CPUTime);
671 }
672
673 /* factorize matrices */
674 if(PrecondProc == ILUPrecond || PrecondProcC == ILUPrecond) {
675 /* start of time counting */
676 ClockFactor = 1.0 / (double)CLOCKS_PER_SEC;
677 BegClock = clock();
678
679 printf("\n");
680 printf("Incomplete factorization of matrices ...\n");
681 if(PrecondProc == ILUPrecond
682 && SmoothProc != JacobiIter
683 && SmoothProc != SORForwIter
684 && SmoothProc != SORBackwIter
685 && SmoothProc != SSORIter) {
686 for (Level = NoLevels - 1; Level > 0; Level--)
687 ILUFactor(&L[Level]);
688 }
689 if(PrecondProcC == ILUPrecond
690 && SolvProc != JacobiIter
691 && SolvProc != SORForwIter
692 && SolvProc != SORBackwIter
693 && SolvProc != SSORIter) {
694 ILUFactor(&L[0]);
695 }
696
697 /* end of time counting and geting out the CPU time */
698 EndClock = clock();
699 CPUTime = (double)(EndClock - BegClock) * ClockFactor;
700 printf("\n");
701 printf(" CPU time: %7.2f s\n", CPUTime);
702 }
703
704 /* solving of system of equations by nested multigrid method */
705 if (MLSolverId == MGIterId && NestedMG) {
706 /* initialisation of vector of unknowns */
707 V_SetAllCmp(&uw[0], 0.0);
708
709 /* approximation of solution by nested multigrid method */
710 printf("\n");
711 printf("Doing nested multigrid iterations ...\n");
712 NestedMGIter(NoLevels, L, uw, fr, R, P, Gamma,
713 SmoothProc, Nu1, Nu2, PrecondProc, Omega,
714 SolvProc, NuC, PrecondProcC ,OmegaC);
715
716 AccBeg = GetLastAccuracy();
717 } else {
718 /* initialisation of vector of unknowns */
719 V_SetAllCmp(&uw[NoLevels - 1], 0.0);
720
721 AccBeg = 1.0;
722 }
723
724 /* start of time counting */
725 ClockFactor = 1.0 / (double)CLOCKS_PER_SEC;
726 BegClock = clock();
727
728 /* solving of system of equations */
729 switch (MLSolverId) {
730 case MGIterId:
731 printf("\n");
732 printf("Doing multigrid iterations ... \n\n");
733 MGIter(NoLevels, L, uw, fr, R, P, MaxIter, Gamma,
734 SmoothProc, Nu1, Nu2, PrecondProc, Omega,
735 SolvProc, NuC, PrecondProcC, OmegaC);
736 break;
737 case MGPCGIterId:
738 printf("\n");
739 printf("Doing multigrid preconditioned CG iterations ... \n\n");
740 MGPCGIter(NoLevels, L, uw, fr, R, P, MaxIter, NoMGIter, Gamma,
741 SmoothProc, Nu1, Nu2, PrecondProc, Omega,
742 SolvProc, NuC, PrecondProcC, OmegaC);
743 break;
744 case BPXPCGIterId:
745 printf("\n");
746 printf("Doing BPX preconditioned CG iterations ... \n\n");
747 BPXPCGIter(NoLevels, L, uw, fr, R, P, MaxIter,
748 SmoothProc, Nu1, PrecondProc, Omega,
749 SolvProc, NuC, PrecondProcC, OmegaC);
750 break;
751 default:
752 break;
753 }
754
755 /* end of time counting and geting out of CPU time */
756 EndClock = clock();
757 CPUTime = (double)(EndClock - BegClock) * ClockFactor;
758 printf("\n");
759 printf(" CPU time: %7.2f s\n", CPUTime);
760
761 AccEnd = GetLastAccuracy();
762 NoIter = GetLastNoIter();
763
764 /* computing of middle contraction rates */
765 if (NoIter > 0)
766 ContrRatePerMGIter = pow(AccEnd / AccBeg, 1.0 / (double)NoIter);
767 else
768 ContrRatePerMGIter = 0.0;
769 if (CPUTime > DBL_EPSILON)
770 ContrRatePerSec = pow(AccEnd / AccBeg, 1.0 / CPUTime);
771 else
772 ContrRatePerSec = 0.0;
773 printf("\n");
774 printf("Middle contraction rate\n");
775 printf(" referred to one iteration: %10.3e\n", ContrRatePerMGIter);
776 printf(" referred to 1 s CPU time: %10.3e\n", ContrRatePerSec);
777 }
778
MGResOutput1DDirich(int NoLevels)779 static void MGResOutput1DDirich(int NoLevels)
780 /* solution output for 1D case, Dirichlet boundary condition */
781 {
782 char Key[2];
783 double x;
784 size_t Ind, Dim;
785
786 fprintf(stderr, "\n");
787 do {
788 fprintf(stderr, "Output solution? (y/n) ");
789 scanf("%1s", Key);
790 Key[0] = toupper(Key[0]);
791 } while (Key[0] != 'Y' && Key[0] != 'N');
792 if (Key[0] == 'Y') {
793 printf("\n");
794 printf(" x u(x) (node)\n");
795 printf("---------------------------------\n");
796 printf(" 0.00000 0.00000\n");
797 Dim = V_GetDim(&uw[NoLevels - 1]);
798 for (Ind = 1; Ind <= Dim; Ind++) {
799 x = (double)Ind / (double)(Dim + 1);
800 printf(" %8.5f %8.5f (%lu)\n", x, V_GetCmp(&uw[NoLevels - 1],Ind),
801 (unsigned long)Ind);
802 }
803 printf(" 1.00000 0.00000\n");
804 printf("---------------------------------\n");
805 }
806 }
807
GenL1DDirich(int Level,size_t * Dim)808 static void GenL1DDirich(int Level, size_t *Dim)
809 /* generation of matrix for 1D case, Dirichlet boundary
810 condition */
811 {
812 double h;
813 size_t Row;
814
815 if (LASResult() == LASOK) {
816 /* computing of the grid parameter */
817 h = 1.0 / (double)(Dim[Level] + 1);
818
819 #ifdef SYM_STOR
820
821 /* setting of matrix elements */
822 for (Row = 1; Row < Dim[Level]; Row++) {
823 Q_SetLen(&L[Level], Row, 2);
824 if (LASResult() == LASOK) {
825 Q_SetEntry(&L[Level], Row, 0, Row, 2.0 / h);
826 Q_SetEntry(&L[Level], Row, 1, Row + 1, -1.0 / h);
827 }
828 }
829 Row = Dim[Level];
830 Q_SetLen(&L[Level], Row, 1);
831 if (LASResult() == LASOK) {
832 Q_SetEntry(&L[Level], Row, 0, Row, 2.0 / h);
833 }
834
835 #else
836
837 Row = 1;
838 Q_SetLen(&L[Level], Row, 2);
839 if (LASResult() == LASOK) {
840 Q_SetEntry(&L[Level], Row, 0, Row, 2.0 / h);
841 Q_SetEntry(&L[Level], Row, 1, Row + 1, -1.0 / h);
842 }
843 for (Row = 2; Row < Dim[Level]; Row++) {
844 Q_SetLen(&L[Level], Row, 3);
845 if (LASResult() == LASOK) {
846 Q_SetEntry(&L[Level], Row, 0, Row, 2.0 / h);
847 Q_SetEntry(&L[Level], Row, 1, Row - 1, -1.0 / h);
848 Q_SetEntry(&L[Level], Row, 2, Row + 1, -1.0 / h);
849 }
850 }
851 Row = Dim[Level];
852 Q_SetLen(&L[Level], Row, 2);
853 if (LASResult() == LASOK) {
854 Q_SetEntry(&L[Level], Row, 0, Row, 2.0 / h);
855 Q_SetEntry(&L[Level], Row, 1, Row - 1, -1.0 / h);
856 }
857
858 #endif /* SYM_STOR */
859
860 }
861 }
862
Genf1DDirich(int Level,size_t * Dim)863 static void Genf1DDirich(int Level, size_t *Dim)
864 /* generation of right hand side for 1D case, Dirichlet boundary condition */
865 {
866 double h;
867 size_t Ind;
868
869 if (LASResult() == LASOK) {
870 /* computing of the grid parameter */
871 h = 1.0 / (double)(Dim[Level] + 1);
872
873 for (Ind = 1; Ind <= Dim[Level]; Ind++) {
874 V_SetCmp(&fr[Level], Ind, 1.0 * h);
875 }
876 }
877 }
878
GenRSimple1DDirich(int Level,size_t * Dim)879 static void GenRSimple1DDirich(int Level, size_t *Dim)
880 /* generation of simple restriction operator (as matrix) for 1D case,
881 Dirichlet boundary condition */
882 {
883 size_t Row;
884
885 if (LASResult() == LASOK) {
886 for (Row = 1; Row <= Dim[Level]; Row++) {
887 M_SetLen(&R[Level], Row, 1);
888 if (LASResult() == LASOK) {
889 M_SetEntry(&R[Level], Row, 0, 2 * Row, 2.0);
890 }
891 }
892 }
893 }
894
GenRWeight1DDirich(int Level,size_t * Dim)895 static void GenRWeight1DDirich(int Level, size_t *Dim)
896 /* generation of weighted restriction operator (as matrix) for 1D case,
897 Dirichlet boundary condition */
898 {
899 size_t Row;
900
901 if (LASResult() == LASOK) {
902 for (Row = 1; Row <= Dim[Level]; Row++) {
903 M_SetLen(&R[Level], Row, 3);
904 if (LASResult() == LASOK) {
905 M_SetEntry(&R[Level], Row, 0, 2 * Row - 1, 0.5);
906 M_SetEntry(&R[Level], Row, 1, 2 * Row, 1.0);
907 M_SetEntry(&R[Level], Row, 2, 2 * Row + 1, 0.5);
908 }
909 }
910 }
911 }
912
GenP1DDirich(int Level,size_t * Dim)913 static void GenP1DDirich(int Level, size_t *Dim)
914 /* generation of prolongation operator (as matrix) for 1D case, Dirichlet
915 boundary condition */
916 {
917 size_t Clm;
918
919 if (LASResult() == LASOK) {
920 for (Clm = 1; Clm <= Dim[Level - 1]; Clm++) {
921 M_SetLen(&P[Level], Clm, 3);
922 if (LASResult() == LASOK) {
923 M_SetEntry(&P[Level], Clm, 0, 2 * Clm - 1, 0.5);
924 M_SetEntry(&P[Level], Clm, 1, 2 * Clm, 1.0);
925 M_SetEntry(&P[Level], Clm, 2, 2 * Clm + 1, 0.5);
926 }
927 }
928 }
929 }
930
GenL2DDirich(int Level,size_t * Dim1)931 static void GenL2DDirich(int Level, size_t *Dim1)
932 /* generation of matrix for 2D case, Dirichlet boundary
933 condition */
934 {
935 size_t BlockRow, RowBasis, Row;
936
937 if (LASResult() == LASOK) {
938 #ifdef SYM_STOR
939
940 /* setting of matrix elements */
941 for (BlockRow = 1; BlockRow < Dim1[Level]; BlockRow++) {
942 RowBasis = Dim1[Level] * (BlockRow - 1);
943 for (Row = RowBasis + 1; Row < RowBasis + Dim1[Level]; Row++) {
944 Q_SetLen(&L[Level], Row, 3);
945 if (LASResult() == LASOK) {
946 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
947 Q_SetEntry(&L[Level], Row, 1, Row + 1, -1.0);
948 Q_SetEntry(&L[Level], Row, 2, Row + Dim1[Level], -1.0);
949 }
950 }
951 Row = RowBasis + Dim1[Level];
952 Q_SetLen(&L[Level], Row, 2);
953 if (LASResult() == LASOK) {
954 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
955 Q_SetEntry(&L[Level], Row, 1, Row + Dim1[Level], -1.0);
956 }
957 }
958 BlockRow = Dim1[Level];
959 RowBasis = Dim1[Level] * (BlockRow - 1);
960 for (Row = RowBasis + 1; Row < RowBasis + Dim1[Level]; Row++) {
961 Q_SetLen(&L[Level], Row, 2);
962 if (LASResult() == LASOK) {
963 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
964 Q_SetEntry(&L[Level], Row, 1, Row + 1, -1.0);
965 }
966 }
967 Row = RowBasis + Dim1[Level];
968 Q_SetLen(&L[Level], Row, 1);
969 if (LASResult() == LASOK) {
970 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
971 }
972
973 #else
974
975 BlockRow = 1;
976 RowBasis = Dim1[Level] * (BlockRow - 1);
977 Row = RowBasis + 1;
978 Q_SetLen(&L[Level], Row, 3);
979 if (LASResult() == LASOK) {
980 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
981 Q_SetEntry(&L[Level], Row, 1, Row + 1, -1.0);
982 Q_SetEntry(&L[Level], Row, 2, Row + Dim1[Level], -1.0);
983 }
984 for (Row = RowBasis + 2; Row < RowBasis + Dim1[Level]; Row++) {
985 Q_SetLen(&L[Level], Row, 4);
986 if (LASResult() == LASOK) {
987 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
988 Q_SetEntry(&L[Level], Row, 1, Row - 1, -1.0);
989 Q_SetEntry(&L[Level], Row, 2, Row + 1, -1.0);
990 Q_SetEntry(&L[Level], Row, 3, Row + Dim1[Level], -1.0);
991 }
992 }
993 Row = RowBasis + Dim1[Level];
994 Q_SetLen(&L[Level], Row, 3);
995 if (LASResult() == LASOK) {
996 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
997 Q_SetEntry(&L[Level], Row, 1, Row - 1, -1.0);
998 Q_SetEntry(&L[Level], Row, 2, Row + Dim1[Level], -1.0);
999 }
1000 for (BlockRow = 2; BlockRow < Dim1[Level]; BlockRow++) {
1001 RowBasis = Dim1[Level] * (BlockRow - 1);
1002 Row = RowBasis + 1;
1003 Q_SetLen(&L[Level], Row, 4);
1004 if (LASResult() == LASOK) {
1005 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
1006 Q_SetEntry(&L[Level], Row, 1, Row - Dim1[Level], -1.0);
1007 Q_SetEntry(&L[Level], Row, 2, Row + 1, -1.0);
1008 Q_SetEntry(&L[Level], Row, 3, Row + Dim1[Level], -1.0);
1009 }
1010 for (Row = RowBasis + 2; Row < RowBasis + Dim1[Level]; Row++) {
1011 Q_SetLen(&L[Level], Row, 5);
1012 if (LASResult() == LASOK) {
1013 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
1014 Q_SetEntry(&L[Level], Row, 1, Row - 1, -1.0);
1015 Q_SetEntry(&L[Level], Row, 2, Row - Dim1[Level], -1.0);
1016 Q_SetEntry(&L[Level], Row, 3, Row + 1, -1.0);
1017 Q_SetEntry(&L[Level], Row, 4, Row + Dim1[Level], -1.0);
1018 }
1019 }
1020 Row = RowBasis + Dim1[Level];
1021 Q_SetLen(&L[Level], Row, 4);
1022 if (LASResult() == LASOK) {
1023 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
1024 Q_SetEntry(&L[Level], Row, 1, Row - 1, -1.0);
1025 Q_SetEntry(&L[Level], Row, 2, Row - Dim1[Level], -1.0);
1026 Q_SetEntry(&L[Level], Row, 3, Row + Dim1[Level], -1.0);
1027 }
1028 }
1029 BlockRow = Dim1[Level];
1030 RowBasis = Dim1[Level] * (BlockRow - 1);
1031 Row = RowBasis + 1;
1032 Q_SetLen(&L[Level], Row, 3);
1033 if (LASResult() == LASOK) {
1034 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
1035 Q_SetEntry(&L[Level], Row, 1, Row - Dim1[Level], -1.0);
1036 Q_SetEntry(&L[Level], Row, 2, Row + 1, -1.0);
1037 }
1038 for (Row = RowBasis + 2; Row < RowBasis + Dim1[Level]; Row++) {
1039 Q_SetLen(&L[Level], Row, 4);
1040 if (LASResult() == LASOK) {
1041 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
1042 Q_SetEntry(&L[Level], Row, 1, Row - 1, -1.0);
1043 Q_SetEntry(&L[Level], Row, 2, Row - Dim1[Level], -1.0);
1044 Q_SetEntry(&L[Level], Row, 3, Row + 1, -1.0);
1045 }
1046 }
1047 Row = RowBasis + Dim1[Level];
1048 Q_SetLen(&L[Level], Row, 3);
1049 if (LASResult() == LASOK) {
1050 Q_SetEntry(&L[Level], Row, 0, Row, 4.0);
1051 Q_SetEntry(&L[Level], Row, 1, Row - 1, -1.0);
1052 Q_SetEntry(&L[Level], Row, 2, Row - Dim1[Level], -1.0);
1053 }
1054
1055 #endif /* SYM_STOR */
1056
1057 }
1058 }
1059
Genf2DDirich(int Level,size_t * Dim1)1060 static void Genf2DDirich(int Level, size_t *Dim1)
1061 /* generation of right hand side for 2D case, Dirichlet boundary condition */
1062 {
1063 double h, hh;
1064 size_t Block, IndBasis, Ind;
1065
1066 if (LASResult() == LASOK) {
1067 /* computing of grid parameter */
1068 h = 1.0 / (double)(Dim1[Level] + 1);
1069 hh = pow(h, 2.0);
1070
1071 for (Block = 1; Block <= Dim1[Level]; Block++) {
1072 IndBasis = Dim1[Level] * (Block - 1);
1073 for (Ind = IndBasis + 1; Ind <= IndBasis + Dim1[Level]; Ind++) {
1074 V_SetCmp(&fr[Level], Ind, 1.0 * hh);
1075 }
1076 }
1077 }
1078 }
1079
GenRSimple2DDirich(int Level,size_t * Dim1)1080 static void GenRSimple2DDirich(int Level, size_t *Dim1)
1081 /* generation of simple restriction operator (as matrix) for 2D case,
1082 Dirichlet boundary condition */
1083 {
1084 size_t BlockRow, RowBasis, RowInBlock, Row;
1085
1086 if (LASResult() == LASOK) {
1087 for (BlockRow = 1; BlockRow <= Dim1[Level]; BlockRow++) {
1088 RowBasis = Dim1[Level] * (BlockRow - 1);
1089 for (RowInBlock = 1; RowInBlock <= Dim1[Level]; RowInBlock++) {
1090 Row = RowBasis + RowInBlock;
1091 M_SetLen(&R[Level], Row, 1);
1092 if (LASResult() == LASOK) {
1093 M_SetEntry(&R[Level], Row, 0, 2 * RowInBlock
1094 + (2 * BlockRow - 1) * Dim1[Level + 1], 4.0);
1095 }
1096 }
1097 }
1098 }
1099 }
1100
GenRWeight2DDirich(int Level,size_t * Dim1)1101 static void GenRWeight2DDirich(int Level, size_t *Dim1)
1102 /* generation of weighted restriction operator (as matrix) for 2D case,
1103 Dirichlet boundary condition */
1104 {
1105 size_t BlockRow, RowBasis, RowInBlock, Row;
1106
1107 if (LASResult() == LASOK) {
1108 for (BlockRow = 1; BlockRow <= Dim1[Level]; BlockRow++) {
1109 RowBasis = Dim1[Level] * (BlockRow - 1);
1110 for (RowInBlock = 1; RowInBlock <= Dim1[Level]; RowInBlock++) {
1111 Row = RowBasis + RowInBlock;
1112 M_SetLen(&R[Level], Row, 9);
1113 if (LASResult() == LASOK) {
1114 M_SetEntry(&R[Level], Row, 0, 2 * RowInBlock
1115 + (2 * BlockRow - 2) * Dim1[Level + 1] - 1, 0.25);
1116 M_SetEntry(&R[Level], Row, 1, 2 * RowInBlock
1117 + (2 * BlockRow - 2) * Dim1[Level + 1], 0.5);
1118 M_SetEntry(&R[Level], Row, 2, 2 * RowInBlock
1119 + (2 * BlockRow - 2) * Dim1[Level + 1] + 1, 0.25);
1120 M_SetEntry(&R[Level], Row, 3, 2 * RowInBlock
1121 + (2 * BlockRow - 1) * Dim1[Level + 1] - 1, 0.5);
1122 M_SetEntry(&R[Level], Row, 4, 2 * RowInBlock
1123 + (2 * BlockRow - 1) * Dim1[Level + 1], 1.0);
1124 M_SetEntry(&R[Level], Row, 5, 2 * RowInBlock
1125 + (2 * BlockRow - 1) * Dim1[Level + 1] + 1, 0.5);
1126 M_SetEntry(&R[Level], Row, 6, 2 * RowInBlock
1127 + (2 * BlockRow) * Dim1[Level + 1] - 1, 0.25);
1128 M_SetEntry(&R[Level], Row, 7, 2 * RowInBlock
1129 + (2 * BlockRow) * Dim1[Level + 1], 0.5);
1130 M_SetEntry(&R[Level], Row, 8, 2 * RowInBlock
1131 + (2 * BlockRow) * Dim1[Level + 1] + 1, 0.25);
1132 }
1133 }
1134 }
1135 }
1136 }
1137
GenP2DDirich(int Level,size_t * Dim1)1138 static void GenP2DDirich(int Level, size_t *Dim1)
1139 /* generation of prolongation operator (as matrix) for 2D case, Dirichlet
1140 boundary condition */
1141 {
1142 size_t BlockClm, ClmBasis, ClmInBlock, Clm;
1143
1144 if (LASResult() == LASOK) {
1145 for (BlockClm = 1; BlockClm <= Dim1[Level - 1]; BlockClm++) {
1146 ClmBasis = Dim1[Level - 1] * (BlockClm - 1);
1147 for (ClmInBlock = 1; ClmInBlock <= Dim1[Level - 1]; ClmInBlock++) {
1148 Clm = ClmBasis + ClmInBlock;
1149 M_SetLen(&P[Level], Clm, 9);
1150 if (LASResult() == LASOK) {
1151 M_SetEntry(&P[Level], Clm, 0, 2 * ClmInBlock
1152 + (2 * BlockClm - 2) * Dim1[Level] - 1, 0.25);
1153 M_SetEntry(&P[Level], Clm, 1, 2 * ClmInBlock
1154 + (2 * BlockClm - 2) * Dim1[Level], 0.5);
1155 M_SetEntry(&P[Level], Clm, 2, 2 * ClmInBlock
1156 + (2 * BlockClm - 2) * Dim1[Level] + 1, 0.25);
1157 M_SetEntry(&P[Level], Clm, 3, 2 * ClmInBlock
1158 + (2 * BlockClm - 1) * Dim1[Level] - 1, 0.5);
1159 M_SetEntry(&P[Level], Clm, 4, 2 * ClmInBlock
1160 + (2 * BlockClm - 1) * Dim1[Level], 1.0);
1161 M_SetEntry(&P[Level], Clm, 5, 2 * ClmInBlock
1162 + (2 * BlockClm - 1) * Dim1[Level] + 1, 0.5);
1163 M_SetEntry(&P[Level], Clm, 6, 2 * ClmInBlock
1164 + (2 * BlockClm) * Dim1[Level] - 1, 0.25);
1165 M_SetEntry(&P[Level], Clm, 7, 2 * ClmInBlock
1166 + (2 * BlockClm) * Dim1[Level], 0.5);
1167 M_SetEntry(&P[Level], Clm, 8, 2 * ClmInBlock
1168 + (2 * BlockClm) * Dim1[Level] + 1, 0.25);
1169 }
1170 }
1171 }
1172 }
1173 }
1174
IterStatus(int Iter,double rNorm,double bNorm,IterIdType IterId)1175 static void IterStatus(int Iter, double rNorm, double bNorm, IterIdType IterId)
1176 /* put out accuracy after each multigrid iteration */
1177 {
1178 if (IterId == MGIterId || IterId == MGPCGIterId || IterId == BPXPCGIterId) {
1179 printf("%3d. iteration ... accuracy = ", Iter);
1180 if (!IsZero(bNorm))
1181 printf("%11.4e\n", rNorm / bNorm);
1182 else
1183 printf(" ---\n");
1184 }
1185 }
1186