Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
SolGMGPoisson.c
Go to the documentation of this file.
1
14#include <time.h>
15#include <math.h>
16
17#include "fasp.h"
18#include "fasp_functs.h"
19
20/*---------------------------------*/
21/*-- Declare Private Functions --*/
22/*---------------------------------*/
23
24#include "PreGMG.inl"
25
26/*---------------------------------*/
27/*-- Public Functions --*/
28/*---------------------------------*/
29
48INT fasp_poisson_gmg1d (REAL *u,
49 REAL *b,
50 const INT nx,
51 const INT maxlevel,
52 const REAL rtol,
53 const SHORT prtlvl)
54{
55 const REAL atol = 1.0E-15;
56 const INT max_itr_num = 100;
57
58 REAL *u0, *r0, *b0;
59 REAL norm_r, norm_r0, norm_r1, factor, error = BIGREAL;
60 INT i, *level, count = 0;
61 REAL AMG_start = 0, AMG_end;
62
63#if DEBUG_MODE > 0
64 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
65 printf("### DEBUG: nx=%d, maxlevel=%d\n", nx, maxlevel);
66#endif
67
68 if ( prtlvl > PRINT_NONE ) {
69 fasp_gettime(&AMG_start);
70 printf("Num of DOF's: %d\n", nx+1);
71 }
72
73 // set level
74 level = (INT *)malloc((maxlevel+2)*sizeof(INT));
75 level[0] = 0; level[1] = nx+1;
76 for (i = 1; i < maxlevel; i++) {
77 level[i+1] = level[i]+(level[i]-level[i-1]+1)/2;
78 }
79 level[maxlevel+1] = level[maxlevel]+1;
80
81 // set u0, b0, r0
82 u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
83 b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
84 r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
85
86 fasp_darray_set(level[maxlevel], u0, 0.0);
87 fasp_darray_set(level[maxlevel], b0, 0.0);
88 fasp_darray_set(level[maxlevel], r0, 0.0);
89
90 fasp_darray_cp(nx, u, u0);
91 fasp_darray_cp(nx, b, b0);
92
93 // compute initial l2 norm of residue
94 fasp_darray_set(level[1], r0, 0.0);
95 residual1d(u0, b0, r0, 0, level);
96 norm_r0 = l2norm(r0, level, 0);
97 norm_r1 = norm_r0;
98 if (norm_r0 < atol) goto FINISHED;
99
100 if ( prtlvl > PRINT_SOME ){
101 printf("-----------------------------------------------------------\n");
102 printf("It Num | ||r||/||b|| | ||r|| | Conv. Factor\n");
103 printf("-----------------------------------------------------------\n");
104 }
105
106 // GMG solver of V-cycle
107 while (count < max_itr_num) {
108 count++;
109 mg1d(u0, b0, level, 0, maxlevel);
110 residual1d(u0, b0, r0, 0, level);
111 norm_r = l2norm(r0, level, 0);
112 factor = norm_r/norm_r1;
113 error = norm_r / norm_r0;
114 norm_r1 = norm_r;
115 if ( prtlvl > PRINT_SOME ){
116 printf("%6d | %13.6e | %13.6e | %10.4f\n",count,error,norm_r,factor);
117 }
118 if (error < rtol || norm_r < atol) break;
119 }
120
121 if ( prtlvl > PRINT_NONE ){
122 if (count >= max_itr_num) {
123 printf("### WARNING: V-cycle failed to converge.\n");
124 }
125 else {
126 printf("Num of V-cycle's: %d, Relative Residual = %e.\n", count, error);
127 }
128 }
129
130 // Update u
131 fasp_darray_cp(level[1], u0, u);
132
133 // print out CPU time if needed
134 if ( prtlvl > PRINT_NONE ) {
135 fasp_gettime(&AMG_end);
136 fasp_cputime("GMG totally", AMG_end - AMG_start);
137 }
138
139#if DEBUG_MODE > 0
140 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
141#endif
142
143FINISHED:
144 free(level);
145 free(r0);
146 free(u0);
147 free(b0);
148
149 return count;
150}
151
172INT fasp_poisson_gmg2d (REAL *u,
173 REAL *b,
174 const INT nx,
175 const INT ny,
176 const INT maxlevel,
177 const REAL rtol,
178 const SHORT prtlvl)
179{
180 const REAL atol = 1.0E-15;
181 const INT max_itr_num = 100;
182
183 REAL *u0, *b0, *r0;
184 REAL norm_r, norm_r0, norm_r1, factor, error = BIGREAL;
185 INT i, k, count = 0, *nxk, *nyk, *level;
186 REAL AMG_start = 0, AMG_end;
187
188#if DEBUG_MODE > 0
189 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
190 printf("### DEBUG: nx=%d, ny=%d, maxlevel=%d\n", nx, ny, maxlevel);
191#endif
192
193 if ( prtlvl > PRINT_NONE ) {
194 fasp_gettime(&AMG_start);
195 printf("Num of DOF's: %d\n", (nx+1)*(ny+1));
196 }
197
198 // set nxk, nyk
199 nxk = (INT *)malloc(maxlevel*sizeof(INT));
200 nyk = (INT *)malloc(maxlevel*sizeof(INT));
201 nxk[0] = nx+1; nyk[0] = ny+1;
202 for (k=1;k<maxlevel;k++) {
203 nxk[k] = (int) (nxk[k-1]+1)/2;
204 nyk[k] = (int) (nyk[k-1]+1)/2;
205 }
206
207 // set level
208 level = (INT *)malloc((maxlevel+2)*sizeof(INT));
209 level[0] = 0; level[1] = (nx+1)*(ny+1);
210 for (i = 1; i < maxlevel; i++) {
211 level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1);
212 }
213 level[maxlevel+1] = level[maxlevel]+1;
214
215 // set u0, b0
216 u0 = (REAL *)malloc(level[maxlevel+1]*sizeof(REAL));
217 b0 = (REAL *)malloc(level[maxlevel+1]*sizeof(REAL));
218 r0 = (REAL *)malloc(level[maxlevel+1]*sizeof(REAL));
219
220 fasp_darray_set(level[maxlevel], u0, 0.0);
221 fasp_darray_set(level[maxlevel], b0, 0.0);
222 fasp_darray_set(level[maxlevel], r0, 0.0);
223
224 fasp_darray_cp(level[1], u, u0);
225 fasp_darray_cp(level[1], b, b0);
226
227 // compute initial l2 norm of residue
228 residual2d(u0, b0, r0, 0, level, nxk, nyk);
229 norm_r0 = l2norm(r0, level, 0);
230 norm_r1 = norm_r0;
231 if (norm_r0 < atol) goto FINISHED;
232
233 if ( prtlvl > PRINT_SOME ){
234 printf("-----------------------------------------------------------\n");
235 printf("It Num | ||r||/||b|| | ||r|| | Conv. Factor\n");
236 printf("-----------------------------------------------------------\n");
237 }
238
239 // GMG solver of V-cycle
240 while ( count < max_itr_num ) {
241 count++;
242 mg2d(u0, b0, level, 0, maxlevel, nxk, nyk);
243 residual2d(u0, b0, r0, 0, level, nxk, nyk);
244 norm_r = l2norm(r0, level, 0);
245 error = norm_r / norm_r0;
246 factor = norm_r/norm_r1;
247 norm_r1 = norm_r;
248 if ( prtlvl > PRINT_SOME ){
249 printf("%6d | %13.6e | %13.6e | %10.4f\n",count,error,norm_r,factor);
250 }
251 if ( error < rtol || norm_r < atol ) break;
252 }
253
254 if ( prtlvl > PRINT_NONE ){
255 if (count >= max_itr_num) {
256 printf("### WARNING: V-cycle failed to converge.\n");
257 }
258 else {
259 printf("Num of V-cycle's: %d, Relative Residual = %e.\n", count, error);
260 }
261 }
262
263 // update u
264 fasp_darray_cp(level[1], u0, u);
265
266 // print out CPU time if needed
267 if ( prtlvl > PRINT_NONE ) {
268 fasp_gettime(&AMG_end);
269 fasp_cputime("GMG totally", AMG_end - AMG_start);
270 }
271
272#if DEBUG_MODE > 0
273 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
274#endif
275
276FINISHED:
277 free(level);
278 free(nxk);
279 free(nyk);
280 free(u0);
281 free(b0);
282 free(r0);
283
284 return count;
285}
286
308INT fasp_poisson_gmg3d (REAL *u,
309 REAL *b,
310 const INT nx,
311 const INT ny,
312 const INT nz,
313 const INT maxlevel,
314 const REAL rtol,
315 const SHORT prtlvl)
316{
317 const REAL atol = 1.0E-15;
318 const INT max_itr_num = 100;
319
320 REAL *u0, *r0, *b0;
321 REAL norm_r,norm_r0,norm_r1, factor, error = BIGREAL;
322 INT i, k, count = 0, *nxk, *nyk, *nzk, *level;
323 REAL AMG_start = 0, AMG_end;
324
325#if DEBUG_MODE > 0
326 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
327 printf("### DEBUG: nx=%d, ny=%d, nz=%d, maxlevel=%d\n",
328 nx, ny, nz, maxlevel);
329#endif
330
331 if ( prtlvl > PRINT_NONE ) {
332 fasp_gettime(&AMG_start);
333 printf("Num of DOF's: %d\n", (nx+1)*(ny+1)*(nz+1));
334 }
335
336 // set nxk, nyk, nzk
337 nxk = (INT *)malloc(maxlevel*sizeof(INT));
338 nyk = (INT *)malloc(maxlevel*sizeof(INT));
339 nzk = (INT *)malloc(maxlevel*sizeof(INT));
340 nxk[0] = nx+1; nyk[0] = ny+1; nzk[0] = nz+1;
341 for(k=1;k<maxlevel;k++){
342 nxk[k] = (int) (nxk[k-1]+1)/2;
343 nyk[k] = (int) (nyk[k-1]+1)/2;
344 nzk[k] = (int) (nyk[k-1]+1)/2;
345 }
346
347 // set level
348 level = (INT *)malloc((maxlevel+2)*sizeof(INT));
349 level[0] = 0; level[1] = (nx+1)*(ny+1)*(nz+1);
350 for (i = 1; i < maxlevel; i++) {
351 level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1)*(nz/pow(2.0,i)+1);
352 }
353 level[maxlevel+1] = level[maxlevel]+1;
354
355 // set u0, b0, r0
356 u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
357 b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
358 r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
359 fasp_darray_set(level[maxlevel], u0, 0.0);
360 fasp_darray_set(level[maxlevel], b0, 0.0);
361 fasp_darray_set(level[maxlevel], r0, 0.0);
362 fasp_darray_cp(level[1], u, u0);
363 fasp_darray_cp(level[1], b, b0);
364
365 // compute initial l2 norm of residue
366 residual3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
367 norm_r0 = l2norm(r0, level, 0);
368 norm_r1 = norm_r0;
369 if (norm_r0 < atol) goto FINISHED;
370
371 if ( prtlvl > PRINT_SOME ){
372 printf("-----------------------------------------------------------\n");
373 printf("It Num | ||r||/||b|| | ||r|| | Conv. Factor\n");
374 printf("-----------------------------------------------------------\n");
375 }
376
377 // GMG solver of V-cycle
378 while (count < max_itr_num) {
379 count++;
380 mg3d(u0, b0, level, 0, maxlevel, nxk, nyk, nzk);
381 residual3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
382 norm_r = l2norm(r0, level, 0);
383 factor = norm_r/norm_r1;
384 error = norm_r / norm_r0;
385 norm_r1 = norm_r;
386 if ( prtlvl > PRINT_SOME ){
387 printf("%6d | %13.6e | %13.6e | %10.4f\n",count,error,norm_r,factor);
388 }
389 if (error < rtol || norm_r < atol) break;
390 }
391
392 if ( prtlvl > PRINT_NONE ){
393 if (count >= max_itr_num) {
394 printf("### WARNING: V-cycle failed to converge.\n");
395 }
396 else {
397 printf("Num of V-cycle's: %d, Relative Residual = %e.\n", count, error);
398 }
399 }
400
401 // update u
402 fasp_darray_cp(level[1], u0, u);
403
404 // print out CPU time if needed
405 if ( prtlvl > PRINT_NONE ) {
406 fasp_gettime(&AMG_end);
407 fasp_cputime("GMG totally", AMG_end - AMG_start);
408 }
409
410FINISHED:
411 free(level);
412 free(nxk);
413 free(nyk);
414 free(nzk);
415 free(r0);
416 free(u0);
417 free(b0);
418
419#if DEBUG_MODE > 0
420 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
421#endif
422
423 return count;
424}
425
442void fasp_poisson_fgmg1d (REAL *u,
443 REAL *b,
444 const INT nx,
445 const INT maxlevel,
446 const REAL rtol,
447 const SHORT prtlvl)
448{
449 const REAL atol = 1.0E-15;
450 REAL *u0, *r0, *b0;
451 REAL norm_r0, norm_r;
452 INT *level;
453 REAL AMG_start = 0, AMG_end;
454 int i;
455
456#if DEBUG_MODE > 0
457 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
458 printf("### DEBUG: nx=%d, maxlevel=%d\n", nx, maxlevel);
459#endif
460
461 if ( prtlvl > PRINT_NONE ) {
462 fasp_gettime(&AMG_start);
463 printf("Num of DOF's: %d\n", (nx+1));
464 }
465
466 // set level
467 level = (INT *)malloc((maxlevel+2)*sizeof(INT));
468 level[0] = 0; level[1] = nx+1;
469 for (i = 1; i < maxlevel; i++) {
470 level[i+1] = level[i]+(level[i]-level[i-1]+1)/2;
471 }
472 level[maxlevel+1] = level[maxlevel]+1;
473
474 // set u0, b0, r0
475 u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
476 b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
477 r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
478 fasp_darray_set(level[maxlevel], u0, 0.0);
479 fasp_darray_set(level[maxlevel], b0, 0.0);
480 fasp_darray_set(level[maxlevel], r0, 0.0);
481 fasp_darray_cp(nx, u, u0);
482 fasp_darray_cp(nx, b, b0);
483
484 // compute initial l2 norm of residue
485 fasp_darray_set(level[1], r0, 0.0);
486 residual1d(u0, b0, r0, 0, level);
487 norm_r0 = l2norm(r0, level, 0);
488 if (norm_r0 < atol) goto FINISHED;
489
490 // Full GMG solver
491 fmg1d(u0, b0, level, maxlevel, nx);
492
493 // update u
494 fasp_darray_cp(level[1], u0, u);
495
496 // print out Relative Residual and CPU time if needed
497 if ( prtlvl > PRINT_NONE ) {
498 fasp_gettime(&AMG_end);
499 fasp_cputime("FGMG totally", AMG_end - AMG_start);
500 residual1d(u0, b0, r0, 0, level);
501 norm_r = l2norm(r0, level, 0);
502 printf("Relative Residual = %e.\n", norm_r/norm_r0);
503 }
504
505FINISHED:
506 free(level);
507 free(r0);
508 free(u0);
509 free(b0);
510
511#if DEBUG_MODE > 0
512 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
513#endif
514
515 return;
516}
517
536void fasp_poisson_fgmg2d (REAL *u,
537 REAL *b,
538 const INT nx,
539 const INT ny,
540 const INT maxlevel,
541 const REAL rtol,
542 const SHORT prtlvl)
543{
544 const REAL atol = 1.0E-15;
545 REAL *u0, *r0, *b0;
546 REAL norm_r0, norm_r;
547 INT *nxk, *nyk, *level;
548 int i, k;
549 REAL AMG_start = 0, AMG_end;
550
551#if DEBUG_MODE > 0
552 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
553 printf("### DEBUG: nx=%d, ny=%d, maxlevel=%d\n", nx, ny, maxlevel);
554#endif
555
556 if ( prtlvl > PRINT_NONE ) {
557 fasp_gettime(&AMG_start);
558 printf("Num of DOF's: %d\n", (nx+1)*(ny+1));
559 }
560
561 // set nxk, nyk
562 nxk = (INT *)malloc(maxlevel*sizeof(INT));
563 nyk = (INT *)malloc(maxlevel*sizeof(INT));
564
565 nxk[0] = nx+1; nyk[0] = ny+1;
566 for(k=1;k<maxlevel;k++) {
567 nxk[k] = (int) (nxk[k-1]+1)/2;
568 nyk[k] = (int) (nyk[k-1]+1)/2;
569 }
570
571 // set level
572 level = (INT *)malloc((maxlevel+2)*sizeof(INT));
573 level[0] = 0; level[1] = (nx+1)*(ny+1);
574 for (i = 1; i < maxlevel; i++) {
575 level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1);
576 }
577 level[maxlevel+1] = level[maxlevel] + 1;
578
579 // set u0, b0, r0
580 u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
581 b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
582 r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
583 fasp_darray_set(level[maxlevel], u0, 0.0);
584 fasp_darray_set(level[maxlevel], b0, 0.0);
585 fasp_darray_set(level[maxlevel], r0, 0.0);
586 fasp_darray_cp(level[1], u, u0);
587 fasp_darray_cp(level[1], b, b0);
588
589 // compute initial l2 norm of residue
590 fasp_darray_set(level[1], r0, 0.0);
591 residual2d(u0, b0, r0, 0, level, nxk, nyk);
592 norm_r0 = l2norm(r0, level, 0);
593 if (norm_r0 < atol) goto FINISHED;
594
595 // FMG solver
596 fmg2d(u0, b0, level, maxlevel, nxk, nyk);
597
598 // update u
599 fasp_darray_cp(level[1], u0, u);
600
601 // print out Relative Residual and CPU time if needed
602 if ( prtlvl > PRINT_NONE ) {
603 fasp_gettime(&AMG_end);
604 fasp_cputime("FGMG totally", AMG_end - AMG_start);
605 residual2d(u0, b0, r0, 0, level, nxk, nyk);
606 norm_r = l2norm(r0, level, 0);
607 printf("Relative Residual = %e.\n", norm_r/norm_r0);
608 }
609
610FINISHED:
611 free(level);
612 free(nxk);
613 free(nyk);
614 free(r0);
615 free(u0);
616 free(b0);
617
618#if DEBUG_MODE > 0
619 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
620#endif
621
622 return;
623}
624
644void fasp_poisson_fgmg3d (REAL *u,
645 REAL *b,
646 const INT nx,
647 const INT ny,
648 const INT nz,
649 const INT maxlevel,
650 const REAL rtol,
651 const SHORT prtlvl)
652{
653 const REAL atol = 1.0E-15;
654 REAL *u0, *r0, *b0;
655 REAL norm_r0, norm_r;
656 INT *nxk, *nyk, *nzk, *level;
657 int i, k;
658 REAL AMG_start = 0, AMG_end;
659
660#if DEBUG_MODE > 0
661 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
662 printf("### DEBUG: nx=%d, ny=%d, nz=%d, maxlevel=%d\n",
663 nx, ny, nz, maxlevel);
664#endif
665
666 if ( prtlvl > PRINT_NONE ) {
667 fasp_gettime(&AMG_start);
668 printf("Num of DOF's: %d\n", (nx+1)*(ny+1)*(nz+1));
669 }
670 // set nxk, nyk, nzk
671 nxk = (INT *)malloc(maxlevel*sizeof(INT));
672 nyk = (INT *)malloc(maxlevel*sizeof(INT));
673 nzk = (INT *)malloc(maxlevel*sizeof(INT));
674
675 nxk[0] = nx+1; nyk[0] = ny+1; nzk[0] = nz+1;
676 for(k=1;k<maxlevel;k++){
677 nxk[k] = (int) (nxk[k-1]+1)/2;
678 nyk[k] = (int) (nyk[k-1]+1)/2;
679 nzk[k] = (int) (nyk[k-1]+1)/2;
680 }
681
682 // set level
683 level = (INT *)malloc((maxlevel+2)*sizeof(INT));
684 level[0] = 0; level[1] = (nx+1)*(ny+1)*(nz+1);
685 for (i = 1; i < maxlevel; i++) {
686 level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1)*(nz/pow(2.0,i)+1);
687 }
688 level[maxlevel+1] = level[maxlevel]+1;
689
690 // set u0, b0, r0
691 u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
692 b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
693 r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
694 fasp_darray_set(level[maxlevel], u0, 0.0);
695 fasp_darray_set(level[maxlevel], b0, 0.0);
696 fasp_darray_set(level[maxlevel], r0, 0.0);
697 fasp_darray_cp(level[1], u, u0);
698 fasp_darray_cp(level[1], b, b0);
699
700 // compute initial l2 norm of residue
701 residual3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
702 norm_r0 = l2norm(r0, level, 0);
703 if (norm_r0 < atol) goto FINISHED;
704
705 // FMG
706 fmg3d(u0, b0, level, maxlevel, nxk, nyk, nzk);
707
708 // update u
709 fasp_darray_cp(level[1], u0, u);
710
711 if ( prtlvl > PRINT_NONE ) {
712 fasp_gettime(&AMG_end);
713 fasp_cputime("FGMG totally", AMG_end - AMG_start);
714 residual3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
715 norm_r = l2norm(r0, level, 0);
716 printf("Relative Residual = %e.\n", norm_r/norm_r0);
717 }
718
719FINISHED:
720 free(level);
721 free(nxk);
722 free(nyk);
723 free(nzk);
724 free(r0);
725 free(u0);
726 free(b0);
727
728#if DEBUG_MODE > 0
729 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
730#endif
731
732 return;
733}
734
754INT fasp_poisson_gmgcg1d (REAL *u,
755 REAL *b,
756 const INT nx,
757 const INT maxlevel,
758 const REAL rtol,
759 const SHORT prtlvl)
760{
761 const REAL atol = 1.0E-15;
762 const INT max_itr_num = 100;
763
764 REAL *u0, *r0, *b0;
765 REAL norm_r0;
766 INT *level;
767 int i, iter = 0;
768 REAL AMG_start = 0, AMG_end;
769
770#if DEBUG_MODE > 0
771 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
772 printf("### DEBUG: nx=%d, maxlevel=%d\n", nx, maxlevel);
773#endif
774
775 if ( prtlvl > PRINT_NONE ) {
776 fasp_gettime(&AMG_start);
777 printf("Num of DOF's: %d\n", (nx+1));
778 }
779 // set level
780 level = (INT *)malloc((maxlevel+2)*sizeof(INT));
781 level[0] = 0; level[1] = nx+1;
782 for (i = 1; i < maxlevel; i++) {
783 level[i+1] = level[i]+(level[i]-level[i-1]+1)/2;
784 }
785 level[maxlevel+1] = level[maxlevel]+1;
786
787 // set u0, b0
788 u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
789 b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
790 r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
791 fasp_darray_set(level[maxlevel], u0, 0.0);
792 fasp_darray_set(level[maxlevel], b0, 0.0);
793 fasp_darray_set(level[maxlevel], r0, 0.0);
794 fasp_darray_cp(nx, u, u0);
795 fasp_darray_cp(nx, b, b0);
796
797 // compute initial l2 norm of residue
798 fasp_darray_set(level[1], r0, 0.0);
799 residual1d(u, b, r0, 0, level);
800 norm_r0 = l2norm(r0, level, 0);
801 if (norm_r0 < atol) goto FINISHED;
802
803 // Preconditioned CG method
804 iter = pcg1d(u0, b0, level, maxlevel, nx, rtol, max_itr_num, prtlvl);
805
806 // Update u
807 fasp_darray_cp(level[1], u0, u);
808
809 // print out CPU time if needed
810 if ( prtlvl > PRINT_NONE ) {
811 fasp_gettime(&AMG_end);
812 fasp_cputime("GMG_PCG totally", AMG_end - AMG_start);
813 }
814
815FINISHED:
816 free(level);
817 free(r0);
818 free(u0);
819 free(b0);
820
821#if DEBUG_MODE > 0
822 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
823#endif
824
825 return iter;
826}
827
849INT fasp_poisson_gmgcg2d (REAL *u,
850 REAL *b,
851 const INT nx,
852 const INT ny,
853 const INT maxlevel,
854 const REAL rtol,
855 const SHORT prtlvl)
856{
857 const REAL atol = 1.0E-15;
858 const INT max_itr_num = 100;
859
860 REAL *u0,*r0,*b0;
861 REAL norm_r0;
862 INT *nxk, *nyk, *level;
863 int i, k, iter = 0;
864 REAL AMG_start = 0, AMG_end;
865
866#if DEBUG_MODE > 0
867 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
868 printf("### DEBUG: nx=%d, ny=%d, maxlevel=%d\n", nx, ny, maxlevel);
869#endif
870
871 if ( prtlvl > PRINT_NONE ) {
872 fasp_gettime(&AMG_start);
873 printf("Num of DOF's: %d\n", (nx+1)*(ny+1));
874 }
875 // set nxk, nyk
876 nxk = (INT *)malloc(maxlevel*sizeof(INT));
877 nyk = (INT *)malloc(maxlevel*sizeof(INT));
878
879 nxk[0] = nx+1; nyk[0] = ny+1;
880 for (k=1;k<maxlevel;k++) {
881 nxk[k] = (int) (nxk[k-1]+1)/2;
882 nyk[k] = (int) (nyk[k-1]+1)/2;
883 }
884
885 // set level
886 level = (INT *)malloc((maxlevel+2)*sizeof(INT));
887 level[0] = 0; level[1] = (nx+1)*(ny+1);
888 for (i = 1; i < maxlevel; i++) {
889 level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1);
890 }
891 level[maxlevel+1] = level[maxlevel]+1;
892
893 // set u0, b0, r0
894 u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
895 b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
896 r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
897 fasp_darray_set(level[maxlevel], u0, 0.0);
898 fasp_darray_set(level[maxlevel], b0, 0.0);
899 fasp_darray_set(level[maxlevel], r0, 0.0);
900 fasp_darray_cp(level[1], u, u0);
901 fasp_darray_cp(level[1], b, b0);
902
903 // compute initial l2 norm of residue
904 fasp_darray_set(level[1], r0, 0.0);
905 residual2d(u0, b0, r0, 0, level, nxk, nyk);
906 norm_r0 = l2norm(r0, level, 0);
907 if (norm_r0 < atol) goto FINISHED;
908
909 // Preconditioned CG method
910 iter = pcg2d(u0, b0, level, maxlevel, nxk,
911 nyk, rtol, max_itr_num, prtlvl);
912
913 // update u
914 fasp_darray_cp(level[1], u0, u);
915
916 // print out CPU time if needed
917 if ( prtlvl > PRINT_NONE ) {
918 fasp_gettime(&AMG_end);
919 fasp_cputime("GMG_PCG totally", AMG_end - AMG_start);
920 }
921
922FINISHED:
923 free(level);
924 free(nxk);
925 free(nyk);
926 free(r0);
927 free(u0);
928 free(b0);
929
930#if DEBUG_MODE > 0
931 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
932#endif
933
934 return iter;
935}
936
959INT fasp_poisson_gmgcg3d (REAL *u,
960 REAL *b,
961 const INT nx,
962 const INT ny,
963 const INT nz,
964 const INT maxlevel,
965 const REAL rtol,
966 const SHORT prtlvl)
967{
968 const REAL atol = 1.0E-15;
969 const INT max_itr_num = 100;
970
971 REAL *u0,*r0,*b0;
972 REAL norm_r0;
973 INT *nxk, *nyk, *nzk, *level;
974 int i, k, iter = 0;
975 REAL AMG_start = 0, AMG_end;
976
977#if DEBUG_MODE > 0
978 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
979 printf("### DEBUG: nx=%d, ny=%d, nz=%d, maxlevel=%d\n",
980 nx, ny, nz, maxlevel);
981#endif
982
983 if ( prtlvl > PRINT_NONE ) {
984 fasp_gettime(&AMG_start);
985 printf("Num of DOF's: %d\n", (nx+1)*(ny+1)*(nz+1));
986 }
987
988 // set nxk, nyk, nzk
989 nxk = (INT *)malloc(maxlevel*sizeof(INT));
990 nyk = (INT *)malloc(maxlevel*sizeof(INT));
991 nzk = (INT *)malloc(maxlevel*sizeof(INT));
992
993 nxk[0] = nx+1; nyk[0] = ny+1; nzk[0] = nz+1;
994 for (k = 1; k < maxlevel; k++ ) {
995 nxk[k] = (int) (nxk[k-1]+1)/2;
996 nyk[k] = (int) (nyk[k-1]+1)/2;
997 nzk[k] = (int) (nyk[k-1]+1)/2;
998 }
999
1000 // set level
1001 level = (INT *)malloc((maxlevel+2)*sizeof(INT));
1002 level[0] = 0; level[1] = (nx+1)*(ny+1)*(nz+1);
1003 for (i = 1; i < maxlevel; i++) {
1004 level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1)*(nz/pow(2.0,i)+1);
1005 }
1006 level[maxlevel+1] = level[maxlevel]+1;
1007
1008 // set u0, b0
1009 u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
1010 b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
1011 r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
1012 fasp_darray_set(level[maxlevel], u0, 0.0);
1013 fasp_darray_set(level[maxlevel], b0, 0.0);
1014 fasp_darray_set(level[maxlevel], r0, 0.0);
1015 fasp_darray_cp(level[1], u, u0);
1016 fasp_darray_cp(level[1], b, b0);
1017
1018 // compute initial l2 norm of residue
1019 residual3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
1020 norm_r0 = l2norm(r0, level, 0);
1021 if (norm_r0 < atol) goto FINISHED;
1022
1023 // Preconditioned CG method
1024 iter = pcg3d(u0, b0, level, maxlevel, nxk,
1025 nyk, nzk, rtol, max_itr_num, prtlvl);
1026
1027 // update u
1028 fasp_darray_cp(level[1], u0, u);
1029
1030 // print out CPU time if needed
1031 if ( prtlvl > PRINT_NONE ) {
1032 fasp_gettime(&AMG_end);
1033 fasp_cputime("GMG_PCG totally", AMG_end - AMG_start);
1034 }
1035
1036FINISHED:
1037 free(level);
1038 free(nxk);
1039 free(nyk);
1040 free(nzk);
1041 free(r0);
1042 free(u0);
1043 free(b0);
1044
1045#if DEBUG_MODE > 0
1046 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1047#endif
1048
1049 return iter;
1050}
1051
1052/*---------------------------------*/
1053/*-- End of File --*/
1054/*---------------------------------*/
fasp_darray_set
void fasp_darray_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: AuxArray.c:41
fasp_darray_cp
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:210
fasp_cputime
void fasp_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: AuxMessage.c:179
fasp_gettime
void fasp_gettime(REAL *time)
Get system time.
Definition: AuxTiming.c:37
fasp_poisson_gmg2d
INT fasp_poisson_gmg2d(REAL *u, REAL *b, const INT nx, const INT ny, const INT maxlevel, const REAL rtol, const SHORT prtlvl)
Solve Ax=b of Poisson 2D equation by Geometric Multigrid Method.
Definition: SolGMGPoisson.c:172
fasp_poisson_gmgcg3d
INT fasp_poisson_gmgcg3d(REAL *u, REAL *b, const INT nx, const INT ny, const INT nz, const INT maxlevel, const REAL rtol, const SHORT prtlvl)
Solve Ax=b of Poisson 3D equation by Geometric Multigrid Method (GMG preconditioned Conjugate Gradien...
Definition: SolGMGPoisson.c:959
fasp_poisson_fgmg1d
void fasp_poisson_fgmg1d(REAL *u, REAL *b, const INT nx, const INT maxlevel, const REAL rtol, const SHORT prtlvl)
Solve Ax=b of Poisson 1D equation by Geometric Multigrid Method (FMG)
Definition: SolGMGPoisson.c:442
fasp_poisson_fgmg2d
void fasp_poisson_fgmg2d(REAL *u, REAL *b, const INT nx, const INT ny, const INT maxlevel, const REAL rtol, const SHORT prtlvl)
Solve Ax=b of Poisson 2D equation by Geometric Multigrid Method (FMG)
Definition: SolGMGPoisson.c:536
fasp_poisson_gmg1d
INT fasp_poisson_gmg1d(REAL *u, REAL *b, const INT nx, const INT maxlevel, const REAL rtol, const SHORT prtlvl)
Solve Ax=b of Poisson 1D equation by Geometric Multigrid Method.
Definition: SolGMGPoisson.c:48
fasp_poisson_fgmg3d
void fasp_poisson_fgmg3d(REAL *u, REAL *b, const INT nx, const INT ny, const INT nz, const INT maxlevel, const REAL rtol, const SHORT prtlvl)
Solve Ax=b of Poisson 3D equation by Geometric Multigrid Method (FMG)
Definition: SolGMGPoisson.c:644
fasp_poisson_gmgcg1d
INT fasp_poisson_gmgcg1d(REAL *u, REAL *b, const INT nx, const INT maxlevel, const REAL rtol, const SHORT prtlvl)
Solve Ax=b of Poisson 1D equation by Geometric Multigrid Method (GMG preconditioned Conjugate Gradien...
Definition: SolGMGPoisson.c:754
fasp_poisson_gmg3d
INT fasp_poisson_gmg3d(REAL *u, REAL *b, const INT nx, const INT ny, const INT nz, const INT maxlevel, const REAL rtol, const SHORT prtlvl)
Solve Ax=b of Poisson 3D equation by Geometric Multigrid Method.
Definition: SolGMGPoisson.c:308
fasp_poisson_gmgcg2d
INT fasp_poisson_gmgcg2d(REAL *u, REAL *b, const INT nx, const INT ny, const INT maxlevel, const REAL rtol, const SHORT prtlvl)
Solve Ax=b of Poisson 2D equation by Geometric Multigrid Method (GMG preconditioned Conjugate Gradien...
Definition: SolGMGPoisson.c:849
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
INT
#define INT
Definition: fasp.h:72
BIGREAL
#define BIGREAL
Some global constants.
Definition: fasp_const.h:255
PRINT_SOME
#define PRINT_SOME
Definition: fasp_const.h:75
PRINT_NONE
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73