Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
gmg_poisson.c
Go to the documentation of this file.
1 
6 #include <time.h>
7 #include <math.h>
8 
9 #include "fasp.h"
10 #include "fasp_functs.h"
11 
12 #include "gmg_util.inl"
13 
14 /*---------------------------------*/
15 /*-- Public Functions --*/
16 /*---------------------------------*/
17 
37  REAL *b,
38  const INT nx,
39  const INT maxlevel,
40  const REAL rtol,
41  const SHORT prtlvl)
42 {
43  const REAL atol = 1.0E-15;
44  const INT max_itr_num = 100;
45 
46  REAL *u0, *r0, *b0;
47  REAL norm_r, norm_r0, norm_r1, error, factor;
48  INT i, *level, count = 0;
49  REAL AMG_start, AMG_end;
50 
51 #if DEBUG_MODE > 0
52  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
53  printf("### DEBUG: nx=%d, maxlevel=%d\n", nx, maxlevel);
54 #endif
55 
56  if ( prtlvl > PRINT_NONE ) {
57  fasp_gettime(&AMG_start);
58  printf("Num of DOF's: %d\n", nx+1);
59  }
60 
61  // set level
62  level = (INT *)malloc((maxlevel+2)*sizeof(REAL));
63  level[0] = 0; level[1] = nx+1;
64  for (i = 1; i < maxlevel; i++) {
65  level[i+1] = level[i]+(level[i]-level[i-1]+1)/2;
66  }
67  level[maxlevel+1] = level[maxlevel]+1;
68 
69  // set u0, b0, r0
70  u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
71  b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
72  r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
73 
74  fasp_array_set(level[maxlevel], u0, 0.0);
75  fasp_array_set(level[maxlevel], b0, 0.0);
76  fasp_array_set(level[maxlevel], r0, 0.0);
77 
78  fasp_array_cp(nx, u, u0);
79  fasp_array_cp(nx, b, b0);
80 
81  // compute initial l2 norm of residue
82  fasp_array_set(level[1], r0, 0.0);
83  compute_r_1d(u0, b0, r0, 0, level);
84  norm_r0 = computenorm(r0, level, 0);
85  norm_r1 = norm_r0;
86  if (norm_r0 < atol) goto FINISHED;
87 
88  if ( prtlvl > PRINT_SOME ){
89  printf("-----------------------------------------------------------\n");
90  printf("It Num | ||r||/||b|| | ||r|| | Conv. Factor\n");
91  printf("-----------------------------------------------------------\n");
92  }
93 
94  // GMG solver of V-cycle
95  while (count < max_itr_num) {
96  count++;
97  multigriditeration1d(u0, b0, level, 0, maxlevel);
98  compute_r_1d(u0, b0, r0, 0, level);
99  norm_r = computenorm(r0, level, 0);
100  factor = norm_r/norm_r1;
101  error = norm_r / norm_r0;
102  norm_r1 = norm_r;
103  if ( prtlvl > PRINT_SOME ){
104  printf("%6d | %13.6e | %13.6e | %10.4f\n",count,error,norm_r,factor);
105  }
106  if (error < rtol || norm_r < atol) break;
107  }
108 
109  if ( prtlvl > PRINT_NONE ){
110  if (count >= max_itr_num) {
111  printf("### WARNING: V-cycle failed to converge.\n");
112  }
113  else {
114  printf("Num of V-cycle's: %d, Relative Residual = %e.\n", count, error);
115  }
116  }
117 
118  // Update u
119  fasp_array_cp(level[1], u0, u);
120 
121  // print out CPU time if needed
122  if ( prtlvl > PRINT_NONE ) {
123  fasp_gettime(&AMG_end);
124  print_cputime("GMG totally", AMG_end - AMG_start);
125  }
126 
127 #if DEBUG_MODE > 0
128  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
129 #endif
130 
131 FINISHED:
132  free(level);
133  free(r0);
134  free(u0);
135  free(b0);
136 
137  return count;
138 }
139 
161  REAL *b,
162  const INT nx,
163  const INT ny,
164  const INT maxlevel,
165  const REAL rtol,
166  const SHORT prtlvl)
167 {
168  const REAL atol = 1.0E-15;
169  const INT max_itr_num = 100;
170 
171  REAL *u0, *b0, *r0;
172  REAL norm_r, norm_r0, norm_r1, error, factor;
173  INT i, k, count = 0, *nxk, *nyk, *level;
174  REAL AMG_start, AMG_end;
175 
176 #if DEBUG_MODE > 0
177  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
178  printf("### DEBUG: nx=%d, ny=%d, maxlevel=%d\n", nx, ny, maxlevel);
179 #endif
180 
181  if ( prtlvl > PRINT_NONE ) {
182  fasp_gettime(&AMG_start);
183  printf("Num of DOF's: %d\n", (nx+1)*(ny+1));
184  }
185 
186  // set nxk, nyk
187  nxk = (INT *)malloc(maxlevel*sizeof(INT));
188  nyk = (INT *)malloc(maxlevel*sizeof(INT));
189  nxk[0] = nx+1; nyk[0] = ny+1;
190  for (k=1;k<maxlevel;k++) {
191  nxk[k] = (int) (nxk[k-1]+1)/2;
192  nyk[k] = (int) (nyk[k-1]+1)/2;
193  }
194 
195  // set level
196  level = (INT *)malloc((maxlevel+2)*sizeof(REAL));
197  level[0] = 0; level[1] = (nx+1)*(ny+1);
198  for (i = 1; i < maxlevel; i++) {
199  level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1);
200  }
201  level[maxlevel+1] = level[maxlevel]+1;
202 
203  // set u0, b0
204  u0 = (REAL *)malloc(level[maxlevel+1]*sizeof(REAL));
205  b0 = (REAL *)malloc(level[maxlevel+1]*sizeof(REAL));
206  r0 = (REAL *)malloc(level[maxlevel+1]*sizeof(REAL));
207 
208  fasp_array_set(level[maxlevel], u0, 0.0);
209  fasp_array_set(level[maxlevel], b0, 0.0);
210  fasp_array_set(level[maxlevel], r0, 0.0);
211 
212  fasp_array_cp(level[1], u, u0);
213  fasp_array_cp(level[1], b, b0);
214 
215  // compute initial l2 norm of residue
216  compute_r_2d(u0, b0, r0, 0, level, nxk, nyk);
217  norm_r0 = computenorm(r0, level, 0);
218  norm_r1 = norm_r0;
219  if (norm_r0 < atol) goto FINISHED;
220 
221  if ( prtlvl > PRINT_SOME ){
222  printf("-----------------------------------------------------------\n");
223  printf("It Num | ||r||/||b|| | ||r|| | Conv. Factor\n");
224  printf("-----------------------------------------------------------\n");
225  }
226 
227  // GMG solver of V-cycle
228  while ( count < max_itr_num ) {
229  count++;
230  multigriditeration2d(u0, b0, level, 0, maxlevel, nxk, nyk);
231  compute_r_2d(u0, b0, r0, 0, level, nxk, nyk);
232  norm_r = computenorm(r0, level, 0);
233  error = norm_r / norm_r0;
234  factor = norm_r/norm_r1;
235  norm_r1 = norm_r;
236  if ( prtlvl > PRINT_SOME ){
237  printf("%6d | %13.6e | %13.6e | %10.4f\n",count,error,norm_r,factor);
238  }
239  if ( error < rtol || norm_r < atol ) break;
240  }
241 
242  if ( prtlvl > PRINT_NONE ){
243  if (count >= max_itr_num) {
244  printf("### WARNING: V-cycle failed to converge.\n");
245  }
246  else {
247  printf("Num of V-cycle's: %d, Relative Residual = %e.\n", count, error);
248  }
249  }
250 
251  // update u
252  fasp_array_cp(level[1], u0, u);
253 
254  // print out CPU time if needed
255  if ( prtlvl > PRINT_NONE ) {
256  fasp_gettime(&AMG_end);
257  print_cputime("GMG totally", AMG_end - AMG_start);
258  }
259 
260 #if DEBUG_MODE > 0
261  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
262 #endif
263 
264 FINISHED:
265  free(level);
266  free(nxk);
267  free(nyk);
268  free(u0);
269  free(b0);
270  free(r0);
271 
272  return count;
273 }
274 
297  REAL *b,
298  const INT nx,
299  const INT ny,
300  const INT nz,
301  const INT maxlevel,
302  const REAL rtol,
303  const SHORT prtlvl)
304 {
305  const REAL atol = 1.0E-15;
306  const INT max_itr_num = 100;
307 
308  REAL *u0,*r0,*b0;
309  REAL norm_r,norm_r0,norm_r1, error, factor;
310  INT i, k, count = 0, *nxk, *nyk, *nzk, *level;
311  REAL AMG_start, AMG_end;
312 
313 #if DEBUG_MODE > 0
314  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
315  printf("### DEBUG: nx=%d, ny=%d, nz=%d, maxlevel=%d\n",
316  nx, ny, nz, maxlevel);
317 #endif
318 
319  if ( prtlvl > PRINT_NONE ) {
320  fasp_gettime(&AMG_start);
321  printf("Num of DOF's: %d\n", (nx+1)*(ny+1)*(nz+1));
322  }
323 
324  // set nxk, nyk, nzk
325  nxk = (INT *)malloc(maxlevel*sizeof(INT));
326  nyk = (INT *)malloc(maxlevel*sizeof(INT));
327  nzk = (INT *)malloc(maxlevel*sizeof(INT));
328  nxk[0] = nx+1; nyk[0] = ny+1; nzk[0] = nz+1;
329  for(k=1;k<maxlevel;k++){
330  nxk[k] = (int) (nxk[k-1]+1)/2;
331  nyk[k] = (int) (nyk[k-1]+1)/2;
332  nzk[k] = (int) (nyk[k-1]+1)/2;
333  }
334 
335  // set level
336  level = (INT *)malloc((maxlevel+2)*sizeof(REAL));
337  level[0] = 0; level[1] = (nx+1)*(ny+1)*(nz+1);
338  for (i = 1; i < maxlevel; i++) {
339  level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1)*(nz/pow(2.0,i)+1);
340  }
341  level[maxlevel+1] = level[maxlevel]+1;
342 
343  // set u0, b0, r0
344  u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
345  b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
346  r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
347  fasp_array_set(level[maxlevel], u0, 0.0);
348  fasp_array_set(level[maxlevel], b0, 0.0);
349  fasp_array_set(level[maxlevel], r0, 0.0);
350  fasp_array_cp(level[1], u, u0);
351  fasp_array_cp(level[1], b, b0);
352 
353  // compute initial l2 norm of residue
354  compute_r_3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
355  norm_r0 = computenorm(r0, level, 0);
356  norm_r1 = norm_r0;
357  if (norm_r0 < atol) goto FINISHED;
358 
359  if ( prtlvl > PRINT_SOME ){
360  printf("-----------------------------------------------------------\n");
361  printf("It Num | ||r||/||b|| | ||r|| | Conv. Factor\n");
362  printf("-----------------------------------------------------------\n");
363  }
364 
365  // GMG solver of V-cycle
366  while (count < max_itr_num) {
367  count++;
368  multigriditeration3d(u0, b0, level, 0, maxlevel, nxk, nyk, nzk);
369  compute_r_3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
370  norm_r = computenorm(r0, level, 0);
371  factor = norm_r/norm_r1;
372  error = norm_r / norm_r0;
373  norm_r1 = norm_r;
374  if ( prtlvl > PRINT_SOME ){
375  printf("%6d | %13.6e | %13.6e | %10.4f\n",count,error,norm_r,factor);
376  }
377  if (error < rtol || norm_r < atol) break;
378  }
379 
380  if ( prtlvl > PRINT_NONE ){
381  if (count >= max_itr_num) {
382  printf("### WARNING: V-cycle failed to converge.\n");
383  }
384  else {
385  printf("Num of V-cycle's: %d, Relative Residual = %e.\n", count, error);
386  }
387  }
388 
389  // update u
390  fasp_array_cp(level[1], u0, u);
391 
392  // print out CPU time if needed
393  if ( prtlvl > PRINT_NONE ) {
394  fasp_gettime(&AMG_end);
395  print_cputime("GMG totally", AMG_end - AMG_start);
396  }
397 
398 FINISHED:
399  free(level);
400  free(nxk);
401  free(nyk);
402  free(nzk);
403  free(r0);
404  free(u0);
405  free(b0);
406 
407 #if DEBUG_MODE > 0
408  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
409 #endif
410 
411  return count;
412 }
413 
432  REAL *b,
433  const INT nx,
434  const INT maxlevel,
435  const REAL rtol,
436  const SHORT prtlvl)
437 {
438  const REAL atol = 1.0E-15;
439  REAL *u0,*r0,*b0;
440  REAL norm_r0, norm_r;
441  int i, *level;
442  REAL AMG_start, AMG_end;
443 
444 #if DEBUG_MODE > 0
445  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
446  printf("### DEBUG: nx=%d, maxlevel=%d\n", nx, maxlevel);
447 #endif
448 
449  if ( prtlvl > PRINT_NONE ) {
450  fasp_gettime(&AMG_start);
451  printf("Num of DOF's: %d\n", (nx+1));
452  }
453 
454  // set level
455  level = (INT *)malloc((maxlevel+2)*sizeof(REAL));
456  level[0] = 0; level[1] = nx+1;
457  for (i = 1; i < maxlevel; i++) {
458  level[i+1] = level[i]+(level[i]-level[i-1]+1)/2;
459  }
460  level[maxlevel+1] = level[maxlevel]+1;
461 
462  // set u0, b0, r0
463  u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
464  b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
465  r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
466  fasp_array_set(level[maxlevel], u0, 0.0);
467  fasp_array_set(level[maxlevel], b0, 0.0);
468  fasp_array_set(level[maxlevel], r0, 0.0);
469  fasp_array_cp(nx, u, u0);
470  fasp_array_cp(nx, b, b0);
471 
472  // compute initial l2 norm of residue
473  fasp_array_set(level[1], r0, 0.0);
474  compute_r_1d(u0, b0, r0, 0, level);
475  norm_r0 = computenorm(r0, level, 0);
476  if (norm_r0 < atol) goto FINISHED;
477 
478  // Full GMG solver
479  fullmultigrid_1d(u0, b0, level, maxlevel, nx);
480 
481  // update u
482  fasp_array_cp(level[1], u0, u);
483 
484  // print out Relative Residual and CPU time if needed
485  if ( prtlvl > PRINT_NONE ) {
486  fasp_gettime(&AMG_end);
487  print_cputime("FGMG totally", AMG_end - AMG_start);
488  compute_r_1d(u0, b0, r0, 0, level);
489  norm_r = computenorm(r0, level, 0);
490  printf("Relative Residual = %e.\n", norm_r/norm_r0);
491  }
492 
493 FINISHED:
494  free(level);
495  free(r0);
496  free(u0);
497  free(b0);
498 
499 #if DEBUG_MODE > 0
500  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
501 #endif
502 
503  return;
504 }
505 
525  REAL *b,
526  const INT nx,
527  const INT ny,
528  const INT maxlevel,
529  const REAL rtol,
530  const SHORT prtlvl)
531 {
532  const REAL atol = 1.0E-15;
533  REAL *u0,*r0,*b0;
534  REAL norm_r0, norm_r;
535  INT i, k, *nxk, *nyk, *level;
536  REAL AMG_start, AMG_end;
537 
538 #if DEBUG_MODE > 0
539  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
540  printf("### DEBUG: nx=%d, ny=%d, maxlevel=%d\n", nx, ny, maxlevel);
541 #endif
542 
543  if ( prtlvl > PRINT_NONE ) {
544  fasp_gettime(&AMG_start);
545  printf("Num of DOF's: %d\n", (nx+1)*(ny+1));
546  }
547 
548  // set nxk, nyk
549  nxk = (INT *)malloc(maxlevel*sizeof(INT));
550  nyk = (INT *)malloc(maxlevel*sizeof(INT));
551 
552  nxk[0] = nx+1; nyk[0] = ny+1;
553  for(k=1;k<maxlevel;k++){
554  nxk[k] = (int) (nxk[k-1]+1)/2;
555  nyk[k] = (int) (nyk[k-1]+1)/2;
556  }
557 
558  // set level
559  level = (INT *)malloc((maxlevel+2)*sizeof(REAL));
560  level[0] = 0; level[1] = (nx+1)*(ny+1);
561  for (i = 1; i < maxlevel; i++) {
562  level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1);
563  }
564  level[maxlevel+1] = level[maxlevel] + 1;
565 
566  // set u0, b0, r0
567  u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
568  b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
569  r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
570  fasp_array_set(level[maxlevel], u0, 0.0);
571  fasp_array_set(level[maxlevel], b0, 0.0);
572  fasp_array_set(level[maxlevel], r0, 0.0);
573  fasp_array_cp(level[1], u, u0);
574  fasp_array_cp(level[1], b, b0);
575 
576  // compute initial l2 norm of residue
577  fasp_array_set(level[1], r0, 0.0);
578  compute_r_2d(u0, b0, r0, 0, level, nxk, nyk);
579  norm_r0 = computenorm(r0, level, 0);
580  if (norm_r0 < atol) goto FINISHED;
581 
582  // FMG solver
583  fullmultigrid_2d(u0, b0, level, maxlevel, nxk, nyk);
584 
585  // update u
586  fasp_array_cp(level[1], u0, u);
587 
588  // print out Relative Residual and CPU time if needed
589  if ( prtlvl > PRINT_NONE ) {
590  fasp_gettime(&AMG_end);
591  print_cputime("FGMG totally", AMG_end - AMG_start);
592  compute_r_2d(u0, b0, r0, 0, level, nxk, nyk);
593  norm_r = computenorm(r0, level, 0);
594  printf("Relative Residual = %e.\n", norm_r/norm_r0);
595  }
596 
597 FINISHED:
598  free(level);
599  free(nxk);
600  free(nyk);
601  free(r0);
602  free(u0);
603  free(b0);
604 
605 #if DEBUG_MODE > 0
606  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
607 #endif
608 
609  return;
610 }
611 
633  REAL *b,
634  const INT nx,
635  const INT ny,
636  const INT nz,
637  const INT maxlevel,
638  const REAL rtol,
639  const SHORT prtlvl)
640 {
641  const REAL atol = 1.0E-15;
642  REAL *u0,*r0,*b0;
643  REAL norm_r0, norm_r;
644  int i, k, *nxk, *nyk, *nzk, *level;
645  REAL AMG_start, AMG_end;
646 
647 #if DEBUG_MODE > 0
648  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
649  printf("### DEBUG: nx=%d, ny=%d, nz=%d, maxlevel=%d\n",
650  nx, ny, nz, maxlevel);
651 #endif
652 
653  if ( prtlvl > PRINT_NONE ) {
654  fasp_gettime(&AMG_start);
655  printf("Num of DOF's: %d\n", (nx+1)*(ny+1)*(nz+1));
656  }
657  // set nxk, nyk, nzk
658  nxk = (INT *)malloc(maxlevel*sizeof(INT));
659  nyk = (INT *)malloc(maxlevel*sizeof(INT));
660  nzk = (INT *)malloc(maxlevel*sizeof(INT));
661 
662  nxk[0] = nx+1; nyk[0] = ny+1; nzk[0] = nz+1;
663  for(k=1;k<maxlevel;k++){
664  nxk[k] = (int) (nxk[k-1]+1)/2;
665  nyk[k] = (int) (nyk[k-1]+1)/2;
666  nzk[k] = (int) (nyk[k-1]+1)/2;
667  }
668 
669  // set level
670  level = (INT *)malloc((maxlevel+2)*sizeof(REAL));
671  level[0] = 0; level[1] = (nx+1)*(ny+1)*(nz+1);
672  for (i = 1; i < maxlevel; i++) {
673  level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1)*(nz/pow(2.0,i)+1);
674  }
675  level[maxlevel+1] = level[maxlevel]+1;
676 
677  // set u0, b0, r0
678  u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
679  b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
680  r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
681  fasp_array_set(level[maxlevel], u0, 0.0);
682  fasp_array_set(level[maxlevel], b0, 0.0);
683  fasp_array_set(level[maxlevel], r0, 0.0);
684  fasp_array_cp(level[1], u, u0);
685  fasp_array_cp(level[1], b, b0);
686 
687  // compute initial l2 norm of residue
688  compute_r_3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
689  norm_r0 = computenorm(r0, level, 0);
690  if (norm_r0 < atol) goto FINISHED;
691 
692  // FMG
693  fullmultigrid_3d(u0, b0, level, maxlevel, nxk, nyk, nzk);
694 
695  // update u
696  fasp_array_cp(level[1], u0, u);
697 
698  if ( prtlvl > PRINT_NONE ) {
699  fasp_gettime(&AMG_end);
700  print_cputime("FGMG totally", AMG_end - AMG_start);
701  compute_r_3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
702  norm_r = computenorm(r0, level, 0);
703  printf("Relative Residual = %e.\n", norm_r/norm_r0);
704  }
705 
706 FINISHED:
707  free(level);
708  free(nxk);
709  free(nyk);
710  free(nzk);
711  free(r0);
712  free(u0);
713  free(b0);
714 
715 #if DEBUG_MODE > 0
716  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
717 #endif
718 
719  return;
720 }
721 
742  REAL *b,
743  const INT nx,
744  const INT maxlevel,
745  const REAL rtol,
746  const SHORT prtlvl)
747 {
748  const REAL atol = 1.0E-15;
749  const INT max_itr_num = 100;
750 
751  REAL *u0,*r0,*b0;
752  REAL norm_r0;
753  int i, *level, iter = 0;
754  REAL AMG_start, AMG_end;
755 
756 #if DEBUG_MODE > 0
757  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
758  printf("### DEBUG: nx=%d, maxlevel=%d\n", nx, maxlevel);
759 #endif
760 
761  if ( prtlvl > PRINT_NONE ) {
762  fasp_gettime(&AMG_start);
763  printf("Num of DOF's: %d\n", (nx+1));
764  }
765  // set level
766  level = (INT *)malloc((maxlevel+2)*sizeof(REAL));
767  level[0] = 0; level[1] = nx+1;
768  for (i = 1; i < maxlevel; i++) {
769  level[i+1] = level[i]+(level[i]-level[i-1]+1)/2;
770  }
771  level[maxlevel+1] = level[maxlevel]+1;
772 
773  // set u0, b0
774  u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
775  b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
776  r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
777  fasp_array_set(level[maxlevel], u0, 0.0);
778  fasp_array_set(level[maxlevel], b0, 0.0);
779  fasp_array_set(level[maxlevel], r0, 0.0);
780  fasp_array_cp(nx, u, u0);
781  fasp_array_cp(nx, b, b0);
782 
783  // compute initial l2 norm of residue
784  fasp_array_set(level[1], r0, 0.0);
785  compute_r_1d(u, b, r0, 0, level);
786  norm_r0 = computenorm(r0, level, 0);
787  if (norm_r0 < atol) goto FINISHED;
788 
789  // Preconditioned CG method
790  iter = pcg_1d(u0, b0, level, maxlevel, nx, rtol, max_itr_num, prtlvl);
791 
792  // Update u
793  fasp_array_cp(level[1], u0, u);
794 
795  // print out CPU time if needed
796  if ( prtlvl > PRINT_NONE ) {
797  fasp_gettime(&AMG_end);
798  print_cputime("GMG_PCG totally", AMG_end - AMG_start);
799  }
800 
801 FINISHED:
802  free(level);
803  free(r0);
804  free(u0);
805  free(b0);
806 
807 #if DEBUG_MODE > 0
808  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
809 #endif
810 
811  return iter;
812 }
813 
836  REAL *b,
837  const INT nx,
838  const INT ny,
839  const INT maxlevel,
840  const REAL rtol,
841  const SHORT prtlvl)
842 {
843  const REAL atol = 1.0E-15;
844  const INT max_itr_num = 100;
845 
846  REAL *u0,*r0,*b0;
847  REAL norm_r0;
848  int i, k, *nxk, *nyk, *level, iter = 0;
849  REAL AMG_start, AMG_end;
850 
851 #if DEBUG_MODE > 0
852  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
853  printf("### DEBUG: nx=%d, ny=%d, maxlevel=%d\n", nx, ny, maxlevel);
854 #endif
855 
856  if ( prtlvl > PRINT_NONE ) {
857  fasp_gettime(&AMG_start);
858  printf("Num of DOF's: %d\n", (nx+1)*(ny+1));
859  }
860  // set nxk, nyk
861  nxk = (INT *)malloc(maxlevel*sizeof(INT));
862  nyk = (INT *)malloc(maxlevel*sizeof(INT));
863 
864  nxk[0] = nx+1; nyk[0] = ny+1;
865  for(k=1;k<maxlevel;k++){
866  nxk[k] = (int) (nxk[k-1]+1)/2;
867  nyk[k] = (int) (nyk[k-1]+1)/2;
868  }
869 
870  // set level
871  level = (INT *)malloc((maxlevel+2)*sizeof(REAL));
872  level[0] = 0; level[1] = (nx+1)*(ny+1);
873  for (i = 1; i < maxlevel; i++) {
874  level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1);
875  }
876  level[maxlevel+1] = level[maxlevel]+1;
877 
878  // set u0, b0, r0
879  u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
880  b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
881  r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
882  fasp_array_set(level[maxlevel], u0, 0.0);
883  fasp_array_set(level[maxlevel], b0, 0.0);
884  fasp_array_set(level[maxlevel], r0, 0.0);
885  fasp_array_cp(level[1], u, u0);
886  fasp_array_cp(level[1], b, b0);
887 
888  // compute initial l2 norm of residue
889  fasp_array_set(level[1], r0, 0.0);
890  compute_r_2d(u0, b0, r0, 0, level, nxk, nyk);
891  norm_r0 = computenorm(r0, level, 0);
892  if (norm_r0 < atol) goto FINISHED;
893 
894  // Preconditioned CG method
895  iter = pcg_2d(u0, b0, level, maxlevel, nxk,
896  nyk, rtol, max_itr_num, prtlvl);
897 
898  // update u
899  fasp_array_cp(level[1], u0, u);
900 
901  // print out CPU time if needed
902  if ( prtlvl > PRINT_NONE ) {
903  fasp_gettime(&AMG_end);
904  print_cputime("GMG_PCG totally", AMG_end - AMG_start);
905  }
906 
907 FINISHED:
908  free(level);
909  free(nxk);
910  free(nyk);
911  free(r0);
912  free(u0);
913  free(b0);
914 
915 #if DEBUG_MODE > 0
916  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
917 #endif
918 
919  return iter;
920 }
921 
945  REAL *b,
946  const INT nx,
947  const INT ny,
948  const INT nz,
949  const INT maxlevel,
950  const REAL rtol,
951  const SHORT prtlvl)
952 {
953  const REAL atol = 1.0E-15;
954  const INT max_itr_num = 100;
955 
956  REAL *u0,*r0,*b0;
957  REAL norm_r0;
958  int i, k, *nxk, *nyk, *nzk, *level, iter = 0;
959  REAL AMG_start, AMG_end;
960 
961 #if DEBUG_MODE > 0
962  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
963  printf("### DEBUG: nx=%d, ny=%d, nz=%d, maxlevel=%d\n",
964  nx, ny, nz, maxlevel);
965 #endif
966 
967  if ( prtlvl > PRINT_NONE ) {
968  fasp_gettime(&AMG_start);
969  printf("Num of DOF's: %d\n", (nx+1)*(ny+1)*(nz+1));
970  }
971 
972  // set nxk, nyk, nzk
973  nxk = (INT *)malloc(maxlevel*sizeof(INT));
974  nyk = (INT *)malloc(maxlevel*sizeof(INT));
975  nzk = (INT *)malloc(maxlevel*sizeof(INT));
976 
977  nxk[0] = nx+1; nyk[0] = ny+1; nzk[0] = nz+1;
978  for(k=1;k<maxlevel;k++){
979  nxk[k] = (int) (nxk[k-1]+1)/2;
980  nyk[k] = (int) (nyk[k-1]+1)/2;
981  nzk[k] = (int) (nyk[k-1]+1)/2;
982  }
983 
984  // set level
985  level = (INT *)malloc((maxlevel+2)*sizeof(REAL));
986  level[0] = 0; level[1] = (nx+1)*(ny+1)*(nz+1);
987  for (i = 1; i < maxlevel; i++) {
988  level[i+1] = level[i]+(nx/pow(2.0,i)+1)*(ny/pow(2.0,i)+1)*(nz/pow(2.0,i)+1);
989  }
990  level[maxlevel+1] = level[maxlevel]+1;
991 
992  // set u0, b0
993  u0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
994  b0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
995  r0 = (REAL *)malloc(level[maxlevel]*sizeof(REAL));
996  fasp_array_set(level[maxlevel], u0, 0.0);
997  fasp_array_set(level[maxlevel], b0, 0.0);
998  fasp_array_set(level[maxlevel], r0, 0.0);
999  fasp_array_cp(level[1], u, u0);
1000  fasp_array_cp(level[1], b, b0);
1001 
1002  // compute initial l2 norm of residue
1003  compute_r_3d(u0, b0, r0, 0, level, nxk, nyk, nzk);
1004  norm_r0 = computenorm(r0, level, 0);
1005  if (norm_r0 < atol) goto FINISHED;
1006 
1007  // Preconditioned CG method
1008  iter = pcg_3d(u0, b0, level, maxlevel, nxk,
1009  nyk, nzk, rtol, max_itr_num, prtlvl);
1010 
1011  // update u
1012  fasp_array_cp(level[1], u0, u);
1013 
1014  // print out CPU time if needed
1015  if ( prtlvl > PRINT_NONE ) {
1016  fasp_gettime(&AMG_end);
1017  print_cputime("GMG_PCG totally", AMG_end - AMG_start);
1018  }
1019 
1020 FINISHED:
1021  free(level);
1022  free(nxk);
1023  free(nyk);
1024  free(nzk);
1025  free(r0);
1026  free(u0);
1027  free(b0);
1028 
1029 #if DEBUG_MODE > 0
1030  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1031 #endif
1032 
1033  return iter;
1034 }
1035 
1036 /*---------------------------------*/
1037 /*-- End of File --*/
1038 /*---------------------------------*/
#define REAL
Definition: fasp.h:67
INT fasp_poisson_gmg_2D(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: gmg_poisson.c:160
#define PRINT_SOME
Definition: fasp_const.h:81
#define INT
Definition: fasp.h:64
void fasp_array_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: array.c:48
INT fasp_poisson_gmg_3D(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: gmg_poisson.c:296
INT fasp_poisson_pcg_gmg_3D(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: gmg_poisson.c:944
Main header file for FASP.
void fasp_poisson_fgmg_3D(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 (Full Multigrid)
Definition: gmg_poisson.c:632
void fasp_gettime(REAL *time)
Get system time.
Definition: timing.c:28
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
void fasp_poisson_fgmg_1D(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 (Full Multigrid)
Definition: gmg_poisson.c:431
void print_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: message.c:165
INT count
INT fasp_poisson_gmg_1D(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: gmg_poisson.c:36
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
INT fasp_poisson_pcg_gmg_2D(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: gmg_poisson.c:835
INT fasp_poisson_pcg_gmg_1D(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: gmg_poisson.c:741
void fasp_poisson_fgmg_2D(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 (Full Multigrid)
Definition: gmg_poisson.c:524