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