Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
PreMGRecurAMLI.c
Go to the documentation of this file.
1
19#include <math.h>
20#include <time.h>
21
22#include "fasp.h"
23#include "fasp_functs.h"
24
25/*---------------------------------*/
26/*-- Declare Private Functions --*/
27/*---------------------------------*/
28
29#include "PreMGRecurAMLI.inl"
30#include "PreMGSmoother.inl"
31#include "PreMGUtil.inl"
32
33/*---------------------------------*/
34/*-- Public Functions --*/
35/*---------------------------------*/
36
58void fasp_solver_amli(AMG_data* mgl, AMG_param* param, INT l)
59{
60 const SHORT amg_type = param->AMG_type;
61 const SHORT prtlvl = param->print_level;
62 const SHORT smoother = param->smoother;
63 const SHORT smooth_order = param->smooth_order;
64 const SHORT coarse_solver = param->coarse_solver;
65 const SHORT degree = param->amli_degree;
66 const REAL relax = param->relaxation;
67 const REAL tol = param->tol * 1e-4;
68 const SHORT ndeg = param->polynomial_degree;
69
70 // local variables
71 REAL alpha = 1.0;
72 REAL* coef = param->amli_coef;
73
74 dvector *b0 = &mgl[l].b, *e0 = &mgl[l].x; // fine level b and x
75 dvector *b1 = &mgl[l + 1].b, *e1 = &mgl[l + 1].x; // coarse level b and x
76
77 dCSRmat* A0 = &mgl[l].A; // fine level matrix
78 dCSRmat* A1 = &mgl[l + 1].A; // coarse level matrix
79
80 const INT m0 = A0->row, m1 = A1->row;
81
82 INT* ordering = mgl[l].cfmark.val; // smoother ordering
83 ILU_data* LU_level = &mgl[l].LU; // fine level ILU decomposition
84 REAL* r = mgl[l].w.val; // work array for residual
85 REAL* r1 = mgl[l + 1].w.val + m1; // work array for residual
86
87 // Schwarz parameters
88 SWZ_param swzparam;
89 if (param->SWZ_levels > 0) {
90 swzparam.SWZ_blksolver = param->SWZ_blksolver;
91 }
92
93#if DEBUG_MODE > 0
94 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
95 printf("### DEBUG: n=%d, nnz=%d\n", mgl[0].A.row, mgl[0].A.nnz);
96#endif
97
98 if (prtlvl >= PRINT_MOST) printf("AMLI level %d, smoother %d.\n", l, smoother);
99
100 if (l < mgl[l].num_levels - 1) {
101
102 // pre smoothing
103 if (l < mgl[l].ILU_levels) {
104
105 fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
106
107 }
108
109 else if (l < mgl->SWZ_levels) {
110
111 switch (mgl[l].Schwarz.SWZ_type) {
112 case SCHWARZ_SYMMETRIC:
113 fasp_dcsr_swz_forward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
114 &mgl[l].b);
115 fasp_dcsr_swz_backward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
116 &mgl[l].b);
117 break;
118 default:
119 fasp_dcsr_swz_forward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
120 &mgl[l].b);
121 break;
122 }
123 }
124
125 else {
126#if MULTI_COLOR_ORDER
127 // printf("fasp_smoother_dcsr_gs_multicolor, %s, %d\n", __FUNCTION__,
128 // __LINE__);
129 fasp_smoother_dcsr_gs_multicolor(&mgl[l].x, &mgl[l].A, &mgl[l].b,
130 param->presmooth_iter, 1);
131#else
132 fasp_dcsr_presmoothing(smoother, A0, b0, e0, param->presmooth_iter, 0,
133 m0 - 1, 1, relax, ndeg, smooth_order, ordering);
134#endif
135 }
136
137 // form residual r = b - A x
138 fasp_darray_cp(m0, b0->val, r);
139 fasp_blas_dcsr_aAxpy(-1.0, A0, e0->val, r);
140
141 // restriction r1 = R*r0
142 switch (amg_type) {
143 case UA_AMG:
144 fasp_blas_dcsr_mxv_agg(&mgl[l].R, r, b1->val);
145 break;
146 default:
147 fasp_blas_dcsr_mxv(&mgl[l].R, r, b1->val);
148 break;
149 }
150
151 // coarse grid correction
152 {
153 INT i;
154
155 fasp_darray_cp(m1, b1->val, r1);
156
157 for (i = 1; i <= degree; i++) {
158 fasp_dvec_set(m1, e1, 0.0);
159 fasp_solver_amli(mgl, param, l + 1);
160
161 // b1 = (coef[degree-i]/coef[degree])*r1 + A1*e1;
162 // First, compute b1 = A1*e1
163 fasp_blas_dcsr_mxv(A1, e1->val, b1->val);
164 // Then, compute b1 = b1 + (coef[degree-i]/coef[degree])*r1
165 fasp_blas_darray_axpy(m1, coef[degree - i] / coef[degree], r1, b1->val);
166 }
167
168 fasp_dvec_set(m1, e1, 0.0);
169 fasp_solver_amli(mgl, param, l + 1);
170 }
171
172 // find the optimal scaling factor alpha
173 fasp_blas_darray_ax(m1, coef[degree], e1->val);
174 if (param->coarse_scaling == ON) {
175 alpha = fasp_blas_darray_dotprod(m1, e1->val, r1) /
176 fasp_blas_dcsr_vmv(A1, e1->val, e1->val);
177 alpha = MIN(alpha, 1.0);
178 }
179
180 // prolongation e0 = e0 + alpha * P * e1
181 switch (amg_type) {
182 case UA_AMG:
183 fasp_blas_dcsr_aAxpy_agg(alpha, &mgl[l].P, e1->val, e0->val);
184 break;
185 default:
186 fasp_blas_dcsr_aAxpy(alpha, &mgl[l].P, e1->val, e0->val);
187 break;
188 }
189
190 // post smoothing
191 if (l < mgl[l].ILU_levels) {
192
193 fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
194
195 }
196
197 else if (l < mgl->SWZ_levels) {
198
199 switch (mgl[l].Schwarz.SWZ_type) {
200 case SCHWARZ_SYMMETRIC:
201 fasp_dcsr_swz_backward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
202 &mgl[l].b);
203 fasp_dcsr_swz_forward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
204 &mgl[l].b);
205 break;
206 default:
207 fasp_dcsr_swz_backward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
208 &mgl[l].b);
209 break;
210 }
211 }
212
213 else {
214#if MULTI_COLOR_ORDER
215 fasp_smoother_dcsr_gs_multicolor(&mgl[l].x, &mgl[l].A, &mgl[l].b,
216 param->postsmooth_iter, -1);
217#else
218 fasp_dcsr_postsmoothing(smoother, A0, b0, e0, param->postsmooth_iter, 0,
219 m0 - 1, -1, relax, ndeg, smooth_order, ordering);
220#endif
221 }
222
223 }
224
225 else { // coarsest level solver
226
227 switch (coarse_solver) {
228
229#if WITH_PARDISO
230 case SOLVER_PARDISO:
231 {
232 /* use Intel MKL PARDISO direct solver on the coarsest level */
233 fasp_pardiso_solve(A0, b0, e0, &mgl[l].pdata, 0);
234 break;
235 }
236#endif
237
238#if WITH_SuperLU
239 case SOLVER_SUPERLU:
240 /* use SuperLU direct solver on the coarsest level */
241 fasp_solver_superlu(A0, b0, e0, 0);
242 break;
243#endif
244
245#if WITH_UMFPACK
246 case SOLVER_UMFPACK:
247 /* use UMFPACK direct solver on the coarsest level */
248 fasp_umfpack_solve(A0, b0, e0, mgl[l].Numeric, 0);
249 break;
250#endif
251
252#if WITH_MUMPS
253 case SOLVER_MUMPS:
254 /* use MUMPS direct solver on the coarsest level */
255 mgl[l].mumps.job = 2;
256 fasp_solver_mumps_steps(A0, b0, e0, &mgl[l].mumps);
257 break;
258#endif
259
260 default:
261 /* use iterative solver on the coarsest level */
262 fasp_coarse_itsolver(A0, b0, e0, tol, prtlvl);
263 }
264 }
265
266#if DEBUG_MODE > 0
267 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
268#endif
269}
270
291void fasp_solver_namli(AMG_data* mgl, AMG_param* param, INT l, INT num_levels)
292{
293 const SHORT amg_type = param->AMG_type;
294 const SHORT prtlvl = param->print_level;
295 const SHORT smoother = param->smoother;
296 const SHORT smooth_order = param->smooth_order;
297 const SHORT coarse_solver = param->coarse_solver;
298 const REAL relax = param->relaxation;
299 const REAL tol = param->tol * 1e-4;
300 const SHORT ndeg = param->polynomial_degree;
301
302 dvector *b0 = &mgl[l].b, *e0 = &mgl[l].x; // fine level b and x
303 dvector *b1 = &mgl[l + 1].b, *e1 = &mgl[l + 1].x; // coarse level b and x
304
305 dCSRmat* A0 = &mgl[l].A; // fine level matrix
306 dCSRmat* A1 = &mgl[l + 1].A; // coarse level matrix
307
308 const INT m0 = A0->row, m1 = A1->row;
309
310 INT* ordering = mgl[l].cfmark.val; // smoother ordering
311 ILU_data* LU_level = &mgl[l].LU; // fine level ILU decomposition
312 REAL* r = mgl[l].w.val; // work array for residual
313
314 dvector uH; // for coarse level correction
315 uH.row = m1;
316 uH.val = mgl[l + 1].w.val + m1;
317
318 // Schwarz parameters
319 SWZ_param swzparam;
320 if (param->SWZ_levels > 0) swzparam.SWZ_blksolver = param->SWZ_blksolver;
321
322#if DEBUG_MODE > 0
323 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
324 printf("### DEBUG: n=%d, nnz=%d\n", mgl[0].A.row, mgl[0].A.nnz);
325#endif
326
327 if (prtlvl >= PRINT_MOST)
328 printf("Nonlinear AMLI level %d, smoother %d.\n", num_levels, smoother);
329
330 if (l < num_levels - 1) {
331
332 // pre smoothing
333 if (l < mgl[l].ILU_levels) {
334
335 fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
336
337 }
338
339 else if (l < mgl->SWZ_levels) {
340
341 switch (mgl[l].Schwarz.SWZ_type) {
342 case SCHWARZ_SYMMETRIC:
343 fasp_dcsr_swz_forward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
344 &mgl[l].b);
345 fasp_dcsr_swz_backward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
346 &mgl[l].b);
347 break;
348 default:
349 fasp_dcsr_swz_forward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
350 &mgl[l].b);
351 break;
352 }
353 }
354
355 else {
356#if MULTI_COLOR_ORDER
357 // printf("fasp_smoother_dcsr_gs_multicolor, %s, %d\n", __FUNCTION__,
358 // __LINE__);
359 fasp_smoother_dcsr_gs_multicolor(&mgl[l].x, &mgl[l].A, &mgl[l].b,
360 param->presmooth_iter, 1);
361#else
362 fasp_dcsr_presmoothing(smoother, A0, b0, e0, param->presmooth_iter, 0,
363 m0 - 1, 1, relax, ndeg, smooth_order, ordering);
364#endif
365 }
366
367 // form residual r = b - A x
368 fasp_darray_cp(m0, b0->val, r);
369 fasp_blas_dcsr_aAxpy(-1.0, A0, e0->val, r);
370
371 // restriction r1 = R*r0
372 switch (amg_type) {
373 case UA_AMG:
374 fasp_blas_dcsr_mxv_agg(&mgl[l].R, r, b1->val);
375 break;
376 default:
377 fasp_blas_dcsr_mxv(&mgl[l].R, r, b1->val);
378 }
379
380 // call nonlinear AMLI-cycle recursively
381 {
382 fasp_dvec_set(m1, e1, 0.0);
383
384 // V-cycle will be enforced when needed !!!
385 if (mgl[l + 1].cycle_type <= 1) {
386
387 fasp_solver_namli(&mgl[l + 1], param, 0, num_levels - 1);
388
389 }
390
391 else { // recursively call preconditioned Krylov method on coarse grid
392
393 precond_data pcdata;
394
395 fasp_param_amg_to_prec(&pcdata, param);
396 pcdata.maxit = 1;
397 pcdata.max_levels = num_levels - 1;
398 pcdata.mgl_data = &mgl[l + 1];
399
400 precond pc;
401 pc.data = &pcdata;
402 pc.fct = fasp_precond_namli;
403
404 fasp_darray_cp(m1, e1->val, uH.val);
405
406 switch (param->nl_amli_krylov_type) {
407 case SOLVER_GCG: // Use GCG
408 Kcycle_dcsr_pgcg(A1, b1, &uH, &pc);
409 break;
410 default: // Use GCR
411 Kcycle_dcsr_pgcr(A1, b1, &uH, &pc);
412 }
413
414 fasp_darray_cp(m1, uH.val, e1->val);
415 }
416 }
417
418 // prolongation e0 = e0 + P*e1
419 switch (amg_type) {
420 case UA_AMG:
421 fasp_blas_dcsr_aAxpy_agg(1.0, &mgl[l].P, e1->val, e0->val);
422 break;
423 default:
424 fasp_blas_dcsr_aAxpy(1.0, &mgl[l].P, e1->val, e0->val);
425 }
426
427 // post smoothing
428 if (l < mgl[l].ILU_levels) {
429
430 fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
431
432 } else if (l < mgl->SWZ_levels) {
433
434 switch (mgl[l].Schwarz.SWZ_type) {
435 case SCHWARZ_SYMMETRIC:
436 fasp_dcsr_swz_backward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
437 &mgl[l].b);
438 fasp_dcsr_swz_forward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
439 &mgl[l].b);
440 break;
441 default:
442 fasp_dcsr_swz_backward(&mgl[l].Schwarz, &swzparam, &mgl[l].x,
443 &mgl[l].b);
444 }
445
446 }
447
448 else {
449#if MULTI_COLOR_ORDER
450 fasp_smoother_dcsr_gs_multicolor(&mgl[l].x, &mgl[l].A, &mgl[l].b,
451 param->postsmooth_iter, -1);
452#else
453 fasp_dcsr_postsmoothing(smoother, A0, b0, e0, param->postsmooth_iter, 0,
454 m0 - 1, -1, relax, ndeg, smooth_order, ordering);
455#endif
456 }
457
458 }
459
460 else { // coarsest level solver
461
462 switch (coarse_solver) {
463
464#if WITH_PARDISO
465 case SOLVER_PARDISO:
466 {
467 /* use Intel MKL PARDISO direct solver on the coarsest level */
468 fasp_pardiso_solve(A0, b0, e0, &mgl[l].pdata, 0);
469 break;
470 }
471#endif
472
473#if WITH_SuperLU
474 case SOLVER_SUPERLU:
475 /* use SuperLU direct solver on the coarsest level */
476 fasp_solver_superlu(A0, b0, e0, 0);
477 break;
478#endif
479
480#if WITH_UMFPACK
481 case SOLVER_UMFPACK:
482 /* use UMFPACK direct solver on the coarsest level */
483 fasp_umfpack_solve(A0, b0, e0, mgl[l].Numeric, 0);
484 break;
485#endif
486
487#if WITH_MUMPS
488 case SOLVER_MUMPS:
489 /* use MUMPS direct solver on the coarsest level */
490 mgl[l].mumps.job = 2;
491 fasp_solver_mumps_steps(A0, b0, e0, &mgl[l].mumps);
492 break;
493#endif
494
495 default:
496 /* use iterative solver on the coarsest level */
497 fasp_coarse_itsolver(A0, b0, e0, tol, prtlvl);
498 }
499 }
500
501#if DEBUG_MODE > 0
502 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
503#endif
504}
505
528void fasp_solver_namli_bsr(AMG_data_bsr* mgl, AMG_param* param, INT l, INT num_levels)
529{
530 const SHORT prtlvl = param->print_level;
531 const SHORT smoother = param->smoother;
532 const SHORT coarse_solver = param->coarse_solver;
533 const REAL relax = param->relaxation;
534 const REAL tol = param->tol;
535 INT i;
536
537 dvector *b0 = &mgl[l].b, *e0 = &mgl[l].x; // fine level b and x
538 dvector *b1 = &mgl[l + 1].b, *e1 = &mgl[l + 1].x; // coarse level b and x
539
540 dBSRmat* A0 = &mgl[l].A; // fine level matrix
541 dBSRmat* A1 = &mgl[l + 1].A; // coarse level matrix
542 const INT m0 = A0->ROW * A0->nb, m1 = A1->ROW * A1->nb;
543
544 ILU_data* LU_level = &mgl[l].LU; // fine level ILU decomposition
545 REAL* r = mgl[l].w.val; // for residual
546
547 dvector uH, bH; // for coarse level correction
548 uH.row = m1;
549 uH.val = mgl[l + 1].w.val + m1;
550 bH.row = m1;
551 bH.val = mgl[l + 1].w.val + 2 * m1;
552
553#if DEBUG_MODE > 0
554 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
555 printf("### DEBUG: n=%d, nnz=%d\n", mgl[0].A.ROW, mgl[0].A.NNZ);
556 // printf("### DEBUG: prtlvl=%d\n", prtlvl);
557 // exit(0);
558#endif
559
560 // REAL start_time, end_time; //! zhaoli
561
562 if (prtlvl >= PRINT_MOST)
563 printf("Nonlinear AMLI: level %d, smoother %d.\n", l, smoother);
564
565 if (l < num_levels - 1) {
566
567 // fasp_gettime(&start_time); //! zhaoli
568
569 // pre smoothing
570 if (l < param->ILU_levels) {
571 fasp_smoother_dbsr_ilu(A0, b0, e0, LU_level);
572 } else {
573 SHORT steps = param->presmooth_iter;
574
575 if (steps > 0) {
576 switch (smoother) {
577 case SMOOTHER_JACOBI:
578 for (i = 0; i < steps; i++)
580 // fasp_smoother_dbsr_jacobi (A0, b0, e0);
583 fasp_smoother_dbsr_jacobi1(A0, b0, e0, mgl[l].diaginv.val);
584 break;
585 case SMOOTHER_GS:
586 if (l == 0) {
587 for (i = 0; i < steps; i++)
588#if BAMG_GS0_DiagInv || 1
589 fasp_smoother_dbsr_gs1(A0, b0, e0, ASCEND, NULL,
590 mgl[l].diaginv.val);
591#else
592 fasp_smoother_dbsr_gs_ascend1(A0, b0, e0);
593#endif
594 } else {
595 for (i = 0; i < steps; i++)
596 fasp_smoother_dbsr_gs1(A0, b0, e0, ASCEND, NULL,
597 mgl[l].diaginv.val);
598 }
599
600 break;
601 case SMOOTHER_SOR:
602 for (i = 0; i < steps; i++)
603 // fasp_smoother_dbsr_sor(A0, b0, e0, ASCEND, NULL, relax);
604 fasp_smoother_dbsr_sor1(A0, b0, e0, ASCEND, NULL,
605 mgl[l].diaginv.val, relax);
606 break;
607 default:
608 printf("### ERROR: Unknown smoother type %d!\n", smoother);
609 fasp_chkerr(ERROR_SOLVER_TYPE, __FUNCTION__);
610 }
611 }
612 }
613
614 // fasp_gettime(&end_time); //! zhaoli
615 // PreSmoother_time_zl += end_time - start_time; //! zhaoli
616
617 // form residual r = b - A x
618 fasp_darray_cp(m0, b0->val, r);
619 fasp_blas_dbsr_aAxpy(-1.0, A0, e0->val, r);
620
621 fasp_blas_dbsr_mxv(&mgl[l].R, r, b1->val);
622
623 // call nonlinear AMLI-cycle recursively
624 {
625 fasp_dvec_set(m1, e1, 0.0);
626
627 // The coarsest problem is solved exactly.
628 // No need to call Krylov method on second coarsest level
629 if (l == num_levels - 2) {
630 fasp_solver_namli_bsr(&mgl[l + 1], param, 0, num_levels - 1);
631 } else { // recursively call preconditioned Krylov method on coarse grid
632 precond_data_bsr pcdata;
633
634 fasp_param_amg_to_precbsr(&pcdata, param);
635 pcdata.maxit = 1;
636 pcdata.max_levels = num_levels - 1;
637 pcdata.mgl_data = &mgl[l + 1];
638
639 precond pc;
640 pc.data = &pcdata;
641 pc.fct = fasp_precond_dbsr_namli;
642
643 fasp_darray_cp(m1, b1->val, bH.val);
644 fasp_darray_cp(m1, e1->val, uH.val);
645
646 const INT maxit = param->amli_degree + 1;
647
648 // fasp_gettime(&start_time); //! zhaoli
649
650 // fasp_solver_dbsr_pcg(A1, &bH, &uH, &pc, param->tol, param->tol *
651 // 1e-8,
652 // maxit, 1, PRINT_NONE);
653
654 // printf("tol = %e\n", param->tol);
655 // exit(0);
656 // fasp_solver_dbsr_pbcgs(A1, &bH, &uH, &pc, param->tol, param->tol *
657 // 1e-8,
658 // maxit, 1, PRINT_NONE);
659
660 fasp_solver_dbsr_pvfgmres(A1, &bH, &uH, &pc, param->tol,
661 param->tol * 1e-8, maxit, MIN(maxit, 30), 1,
662 PRINT_NONE);
663
664 // fasp_gettime(&end_time); //! zhaoli
665 // Krylov_time_zl += end_time - start_time; //! zhaoli
666
667 fasp_darray_cp(m1, bH.val, b1->val);
668 fasp_darray_cp(m1, uH.val, e1->val);
669 }
670 }
671
672 fasp_blas_dbsr_aAxpy(1.0, &mgl[l].P, e1->val, e0->val);
673
674 // fasp_gettime(&start_time); //! zhaoli
675
676 // post smoothing
677 if (l < param->ILU_levels) {
678 fasp_smoother_dbsr_ilu(A0, b0, e0, LU_level);
679 } else {
680 SHORT steps = param->postsmooth_iter;
681
682 if (steps > 0) {
683 switch (smoother) {
684 case SMOOTHER_JACOBI:
685 for (i = 0; i < steps; i++)
687 // fasp_smoother_dbsr_jacobi(A0, b0, e0);
690 fasp_smoother_dbsr_jacobi1(A0, b0, e0, mgl[l].diaginv.val);
691 break;
692 case SMOOTHER_GS:
693 // fasp_smoother_dbsr_gs(A0, b0, e0, ASCEND, NULL);
694 if (l == 0) {
695 for (i = 0; i < steps; i++)
696#if BAMG_GS0_DiagInv || 1
697 fasp_smoother_dbsr_gs1(A0, b0, e0, DESCEND, NULL,
698 mgl[l].diaginv.val);
699#else
700 fasp_smoother_dbsr_gs_descend1(A0, b0, e0);
701#endif
702 } else {
703 for (i = 0; i < steps; i++)
704 fasp_smoother_dbsr_gs1(A0, b0, e0, DESCEND, NULL,
705 mgl[l].diaginv.val);
706 }
707
708 break;
709 case SMOOTHER_SOR:
710 for (i = 0; i < steps; i++)
711 // fasp_smoother_dbsr_sor(A0, b0, e0, ASCEND, NULL, relax);
712 fasp_smoother_dbsr_sor1(A0, b0, e0, DESCEND, NULL,
713 mgl[l].diaginv.val, relax);
714 break;
715 default:
716 printf("### ERROR: Unknown smoother type %d!\n", smoother);
717 fasp_chkerr(ERROR_SOLVER_TYPE, __FUNCTION__);
718 }
719 }
720 }
721
722 // fasp_gettime(&end_time); //! zhaoli
723 // PostSmoother_time_zl += end_time - start_time; //! zhaoli
724
725 }
726
727 else { // coarsest level solver
728 // fasp_gettime(&start_time); //! zhaoli
729
730 switch (coarse_solver) {
731
732#if WITH_PARDISO
733 case SOLVER_PARDISO:
734 {
735 /* use Intel MKL PARDISO direct solver on the coarsest level */
736 fasp_pardiso_solve(&mgl[l].Ac, b0, e0, &mgl[l].pdata, 0);
737 break;
738 }
739#endif
740
741#if WITH_SuperLU
742 case SOLVER_SUPERLU:
743 /* use SuperLU direct solver on the coarsest level */
744 fasp_solver_superlu(&mgl[l].Ac, b0, e0, 0);
745 break;
746#endif
747
748#if WITH_UMFPACK
749 case SOLVER_UMFPACK:
750 /* use UMFPACK direct solver on the coarsest level */
751 fasp_umfpack_solve(&mgl[l].Ac, b0, e0, mgl[l].Numeric, 0);
752 break;
753#endif
754
755#if WITH_MUMPS
756 case SOLVER_MUMPS:
757 /* use MUMPS direct solver on the coarsest level */
758 mgl[l].mumps.job = 2;
759 fasp_solver_mumps_steps(&mgl[l].Ac, b0, e0, &mgl[l].mumps);
760 break;
761#endif
762
763 default:
764 /* use iterative solver on the coarsest level */
765 fasp_coarse_itsolver(&mgl[l].Ac, b0, e0, tol, prtlvl);
766 }
767
768 // fasp_gettime(&end_time); //! zhaoli
769 // Coarsen_time_zl += end_time - start_time; //! zhaoli
770 }
771
772#if DEBUG_MODE > 0
773 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
774#endif
775}
776
791void fasp_amg_amli_coef(const REAL lambda_max,
792 const REAL lambda_min,
793 const INT degree,
794 REAL* coef)
795{
796 const REAL mu0 = 1.0 / lambda_max, mu1 = 1.0 / lambda_min;
797 const REAL c = (sqrt(mu0) + sqrt(mu1)) * (sqrt(mu0) + sqrt(mu1));
798 const REAL a = (4 * mu0 * mu1) / (c);
799
800 const REAL kappa = lambda_max / lambda_min; // condition number
801 const REAL delta = (sqrt(kappa) - 1.0) / (sqrt(kappa) + 1.0);
802 const REAL b = delta * delta;
803
804 if (degree == 0) {
805 coef[0] = 0.5 * (mu0 + mu1);
806 }
807
808 else if (degree == 1) {
809 coef[0] = 0.5 * c;
810 coef[1] = -1.0 * mu0 * mu1;
811 }
812
813 else if (degree > 1) {
814 INT i;
815
816 // allocate memory
817 REAL* work = (REAL*)fasp_mem_calloc(2 * degree - 1, sizeof(REAL));
818 REAL *coef_k, *coef_km1;
819 coef_k = work;
820 coef_km1 = work + degree;
821
822 // get q_k
823 fasp_amg_amli_coef(lambda_max, lambda_min, degree - 1, coef_k);
824 // get q_km1
825 fasp_amg_amli_coef(lambda_max, lambda_min, degree - 2, coef_km1);
826
827 // get coef
828 coef[0] = a - b * coef_km1[0] + (1 + b) * coef_k[0];
829
830 for (i = 1; i < degree - 1; i++) {
831 coef[i] = -b * coef_km1[i] + (1 + b) * coef_k[i] - a * coef_k[i - 1];
832 }
833
834 coef[degree - 1] = (1 + b) * coef_k[degree - 1] - a * coef_k[degree - 2];
835
836 coef[degree] = -a * coef_k[degree - 1];
837
838 // clean memory
839 fasp_mem_free(work);
840 work = NULL;
841 }
842
843 else {
844 printf("### ERROR: Wrong AMLI degree %d!\n", degree);
845 fasp_chkerr(ERROR_INPUT_PAR, __FUNCTION__);
846 }
847
848 return;
849}
850
851/*---------------------------------*/
852/*-- End of File --*/
853/*---------------------------------*/
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_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_chkerr
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
fasp_param_amg_to_prec
void fasp_param_amg_to_prec(precond_data *pcdata, const AMG_param *amgparam)
Set precond_data with AMG_param.
Definition: AuxParam.c:782
fasp_param_amg_to_precbsr
void fasp_param_amg_to_precbsr(precond_data_bsr *pcdata, const AMG_param *amgparam)
Set precond_data_bsr with AMG_param.
Definition: AuxParam.c:848
fasp_dvec_set
void fasp_dvec_set(INT n, dvector *x, const REAL val)
Initialize dvector x[i]=val for i=0:n-1.
Definition: AuxVector.c:222
fasp_blas_darray_dotprod
REAL fasp_blas_darray_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies x and y.
Definition: BlaArray.c:771
fasp_blas_darray_ax
void fasp_blas_darray_ax(const INT n, const REAL a, REAL *x)
x = a*x
Definition: BlaArray.c:43
fasp_blas_darray_axpy
void fasp_blas_darray_axpy(const INT n, const REAL a, const REAL *x, REAL *y)
y = a*x + y
Definition: BlaArray.c:90
fasp_dcsr_swz_forward
void fasp_dcsr_swz_forward(SWZ_data *swzdata, SWZ_param *swzparam, dvector *x, dvector *b)
Schwarz smoother: forward sweep.
Definition: BlaSchwarzSetup.c:218
fasp_dcsr_swz_backward
void fasp_dcsr_swz_backward(SWZ_data *swzdata, SWZ_param *swzparam, dvector *x, dvector *b)
Schwarz smoother: backward sweep.
Definition: BlaSchwarzSetup.c:328
fasp_blas_dbsr_aAxpy
void fasp_blas_dbsr_aAxpy(const REAL alpha, const dBSRmat *A, const REAL *x, REAL *y)
Compute y := alpha*A*x + y.
Definition: BlaSpmvBSR.c:514
fasp_blas_dbsr_mxv
void fasp_blas_dbsr_mxv(const dBSRmat *A, const REAL *x, REAL *y)
Compute y := A*x.
Definition: BlaSpmvBSR.c:1055
fasp_blas_dcsr_mxv_agg
void fasp_blas_dcsr_mxv_agg(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x (nonzeros of A = 1)
Definition: BlaSpmvCSR.c:438
fasp_blas_dcsr_mxv
void fasp_blas_dcsr_mxv(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvCSR.c:242
fasp_blas_dcsr_aAxpy_agg
void fasp_blas_dcsr_aAxpy_agg(const REAL alpha, const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y (nonzeros of A = 1)
Definition: BlaSpmvCSR.c:727
fasp_blas_dcsr_vmv
REAL fasp_blas_dcsr_vmv(const dCSRmat *A, const REAL *x, const REAL *y)
vector-Matrix-vector multiplication alpha = y'*A*x
Definition: BlaSpmvCSR.c:839
fasp_blas_dcsr_aAxpy
void fasp_blas_dcsr_aAxpy(const REAL alpha, const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvCSR.c:494
fasp_smoother_dbsr_gs1
void fasp_smoother_dbsr_gs1(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark, REAL *diaginv)
Gauss-Seidel relaxation.
Definition: ItrSmootherBSR.c:520
fasp_smoother_dbsr_gs_descend1
void fasp_smoother_dbsr_gs_descend1(dBSRmat *A, dvector *b, dvector *u)
Gauss-Seidel relaxation in the descending order.
Definition: ItrSmootherBSR.c:751
fasp_smoother_dbsr_jacobi1
void fasp_smoother_dbsr_jacobi1(dBSRmat *A, dvector *b, dvector *u, REAL *diaginv)
Jacobi relaxation.
Definition: ItrSmootherBSR.c:263
fasp_smoother_dbsr_gs_ascend1
void fasp_smoother_dbsr_gs_ascend1(dBSRmat *A, dvector *b, dvector *u)
Gauss-Seidel relaxation in the ascending order.
Definition: ItrSmootherBSR.c:619
fasp_smoother_dbsr_ilu
void fasp_smoother_dbsr_ilu(dBSRmat *A, dvector *b, dvector *x, void *data)
ILU method as the smoother in solving Au=b with multigrid method.
Definition: ItrSmootherBSR.c:1479
fasp_smoother_dbsr_sor1
void fasp_smoother_dbsr_sor1(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark, REAL *diaginv, REAL weight)
SOR relaxation.
Definition: ItrSmootherBSR.c:1075
fasp_smoother_dcsr_ilu
void fasp_smoother_dcsr_ilu(dCSRmat *A, dvector *b, dvector *x, void *data)
ILU method as a smoother.
Definition: ItrSmootherCSR.c:1279
fasp_solver_dbsr_pvfgmres
INT fasp_solver_dbsr_pvfgmres(dBSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const REAL abstol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Solve "Ax=b" using PFGMRES(right preconditioned) iterative method in which the restart parameter can ...
Definition: KryPvfgmres.c:386
fasp_precond_dbsr_namli
void fasp_precond_dbsr_namli(REAL *r, REAL *z, void *data)
Nonlinear AMLI-cycle AMG preconditioner.
Definition: PreBSR.c:1294
fasp_precond_namli
void fasp_precond_namli(REAL *r, REAL *z, void *data)
Nonlinear AMLI AMG preconditioner.
Definition: PreCSR.c:515
fasp_solver_namli
void fasp_solver_namli(AMG_data *mgl, AMG_param *param, INT l, INT num_levels)
Solve Ax=b with recursive nonlinear AMLI-cycle.
Definition: PreMGRecurAMLI.c:291
fasp_amg_amli_coef
void fasp_amg_amli_coef(const REAL lambda_max, const REAL lambda_min, const INT degree, REAL *coef)
Compute the coefficients of the polynomial used by AMLI-cycle.
Definition: PreMGRecurAMLI.c:791
fasp_solver_namli_bsr
void fasp_solver_namli_bsr(AMG_data_bsr *mgl, AMG_param *param, INT l, INT num_levels)
Solve Ax=b with recursive nonlinear AMLI-cycle.
Definition: PreMGRecurAMLI.c:528
fasp_solver_amli
void fasp_solver_amli(AMG_data *mgl, AMG_param *param, INT l)
Solve Ax=b with recursive AMLI-cycle.
Definition: PreMGRecurAMLI.c:58
fasp_solver_mumps_steps
int fasp_solver_mumps_steps(dCSRmat *ptrA, dvector *b, dvector *u, Mumps_data *mumps)
Solve Ax=b by MUMPS in three steps.
Definition: XtrMumps.c:196
fasp_solver_superlu
int fasp_solver_superlu(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Au=b by SuperLU.
Definition: XtrSuperlu.c:47
fasp.h
Main header file for the FASP project.
MIN
#define MIN(a, b)
Definition: fasp.h:83
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
SMOOTHER_SOR
#define SMOOTHER_SOR
Definition: fasp_const.h:195
SOLVER_PARDISO
#define SOLVER_PARDISO
Definition: fasp_const.h:126
SMOOTHER_JACOBI
#define SMOOTHER_JACOBI
Definition: fasp_const.h:190
PRINT_MOST
#define PRINT_MOST
Definition: fasp_const.h:77
ASCEND
#define ASCEND
Definition: fasp_const.h:249
SOLVER_MUMPS
#define SOLVER_MUMPS
Definition: fasp_const.h:125
ERROR_SOLVER_TYPE
#define ERROR_SOLVER_TYPE
Definition: fasp_const.h:41
SOLVER_SUPERLU
#define SOLVER_SUPERLU
Definition: fasp_const.h:123
SCHWARZ_SYMMETRIC
#define SCHWARZ_SYMMETRIC
Definition: fasp_const.h:158
SMOOTHER_GS
#define SMOOTHER_GS
Definition: fasp_const.h:191
SOLVER_UMFPACK
#define SOLVER_UMFPACK
Definition: fasp_const.h:124
DESCEND
#define DESCEND
Definition: fasp_const.h:250
ERROR_INPUT_PAR
#define ERROR_INPUT_PAR
Definition: fasp_const.h:24
ON
#define ON
Definition of switch.
Definition: fasp_const.h:67
PRINT_NONE
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
SOLVER_GCG
#define SOLVER_GCG
Definition: fasp_const.h:109
UA_AMG
#define UA_AMG
Definition: fasp_const.h:165
AMG_data_bsr
Data for multigrid levels in dBSRmat format.
Definition: fasp_block.h:146
AMG_data_bsr::A
dBSRmat A
pointer to the matrix at level level_num
Definition: fasp_block.h:155
AMG_data_bsr::mumps
Mumps_data mumps
data for MUMPS
Definition: fasp_block.h:245
AMG_data_bsr::b
dvector b
pointer to the right-hand side at level level_num
Definition: fasp_block.h:164
AMG_data_bsr::x
dvector x
pointer to the iterative solution at level level_num
Definition: fasp_block.h:167
AMG_data_bsr::w
dvector w
temporary work space
Definition: fasp_block.h:242
AMG_data_bsr::LU
ILU_data LU
ILU matrix for ILU smoother.
Definition: fasp_block.h:220
AMG_data
Data for AMG methods.
Definition: fasp.h:804
AMG_data::A
dCSRmat A
pointer to the matrix at level level_num
Definition: fasp.h:817
AMG_data::mumps
Mumps_data mumps
data for MUMPS
Definition: fasp.h:866
AMG_data::b
dvector b
pointer to the right-hand side at level level_num
Definition: fasp.h:826
AMG_data::x
dvector x
pointer to the iterative solution at level level_num
Definition: fasp.h:829
AMG_data::w
dvector w
temporary work space
Definition: fasp.h:863
AMG_data::cfmark
ivector cfmark
pointer to the CF marker at level level_num
Definition: fasp.h:840
AMG_data::LU
ILU_data LU
ILU matrix for ILU smoother.
Definition: fasp.h:846
AMG_param
Parameters for AMG methods.
Definition: fasp.h:455
AMG_param::print_level
SHORT print_level
print level for AMG
Definition: fasp.h:461
AMG_param::polynomial_degree
SHORT polynomial_degree
degree of the polynomial smoother
Definition: fasp.h:497
AMG_param::nl_amli_krylov_type
SHORT nl_amli_krylov_type
type of Krylov method used by Nonlinear AMLI cycle
Definition: fasp.h:512
AMG_param::coarse_scaling
SHORT coarse_scaling
switch of scaling of the coarse grid correction
Definition: fasp.h:503
AMG_param::amli_coef
REAL * amli_coef
coefficients of the polynomial used by AMLI cycle
Definition: fasp.h:509
AMG_param::AMG_type
SHORT AMG_type
type of AMG method
Definition: fasp.h:458
AMG_param::tol
REAL tol
stopping tolerance for AMG solver
Definition: fasp.h:467
AMG_param::coarse_solver
SHORT coarse_solver
coarse solver type
Definition: fasp.h:500
AMG_param::relaxation
REAL relaxation
relaxation parameter for Jacobi and SOR smoother
Definition: fasp.h:494
AMG_param::smoother
SHORT smoother
smoother type
Definition: fasp.h:482
AMG_param::amli_degree
SHORT amli_degree
degree of the polynomial used by AMLI cycle
Definition: fasp.h:506
AMG_param::SWZ_blksolver
INT SWZ_blksolver
type of Schwarz block solver
Definition: fasp.h:590
AMG_param::postsmooth_iter
SHORT postsmooth_iter
number of postsmoothers
Definition: fasp.h:491
AMG_param::SWZ_levels
INT SWZ_levels
number of levels use Schwarz smoother
Definition: fasp.h:578
AMG_param::presmooth_iter
SHORT presmooth_iter
number of presmoothers
Definition: fasp.h:488
AMG_param::smooth_order
SHORT smooth_order
smoother order
Definition: fasp.h:485
ILU_data
Data for ILU setup.
Definition: fasp.h:651
Mumps_data::job
INT job
work for MUMPS
Definition: fasp.h:615
SWZ_param
Parameters for Schwarz method.
Definition: fasp.h:430
SWZ_param::SWZ_blksolver
INT SWZ_blksolver
type of Schwarz block solver
Definition: fasp.h:445
dBSRmat
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:34
dBSRmat::NNZ
INT NNZ
number of nonzero sub-blocks in matrix A, NNZ
Definition: fasp_block.h:43
dBSRmat::val
REAL * val
Definition: fasp_block.h:57
dBSRmat::nb
INT nb
dimension of each sub-block
Definition: fasp_block.h:46
dBSRmat::ROW
INT ROW
number of rows of sub-blocks in matrix A, M
Definition: fasp_block.h:37
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360
dvector::row
INT row
number of rows
Definition: fasp.h:357
ivector::val
INT * val
actual vector entries
Definition: fasp.h:374
precond_data_bsr
Data for preconditioners in dBSRmat format.
Definition: fasp_block.h:271
precond_data_bsr::mgl_data
AMG_data_bsr * mgl_data
AMG preconditioner data.
Definition: fasp_block.h:328
precond_data_bsr::max_levels
INT max_levels
max number of AMG levels
Definition: fasp_block.h:283
precond_data_bsr::maxit
INT maxit
max number of iterations of AMG preconditioner
Definition: fasp_block.h:280
precond_data
Data for preconditioners.
Definition: fasp.h:894
precond_data::mgl_data
AMG_data * mgl_data
AMG preconditioner data.
Definition: fasp.h:954
precond_data::max_levels
SHORT max_levels
max number of AMG levels
Definition: fasp.h:906
precond_data::maxit
INT maxit
max number of iterations of AMG preconditioner
Definition: fasp.h:903
precond
Preconditioner data and action.
Definition: fasp.h:1095
precond::data
void * data
data for preconditioner, void pointer
Definition: fasp.h:1098
precond::fct
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1101