Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
amlirecur.c
Go to the documentation of this file.
1 
8 #include <math.h>
9 #include <time.h>
10 
11 #include "fasp.h"
12 #include "fasp_functs.h"
13 // #include "forts_ns.h"
14 
15 #include "mg_util.inl"
16 
17 static SHORT fasp_krylov_cycle_dcsr_pgcg(dCSRmat *, dvector *, dvector *, precond *);
18 static SHORT fasp_krylov_cycle_dcsr_pgcr(dCSRmat *, dvector *, dvector *, precond *);
19 
20 /*---------------------------------*/
21 /*-- Public Functions --*/
22 /*---------------------------------*/
23 
46  AMG_param *param,
47  INT level)
48 {
49  const SHORT amg_type=param->AMG_type;
50  const SHORT prtlvl = param->print_level;
51  const SHORT smoother = param->smoother;
52  const SHORT smooth_order = param->smooth_order;
53  const SHORT coarse_solver = param->coarse_solver;
54  const SHORT degree= param->amli_degree;
55  const REAL relax = param->relaxation;
56  const REAL tol = param->tol*1e-4;
57  const SHORT ndeg = param->polynomial_degree;
58 
59  // local variables
60  REAL alpha = 1.0;
61  REAL * coef = param->amli_coef;
62 
63  dvector *b0 = &mgl[level].b, *e0 = &mgl[level].x; // fine level b and x
64  dvector *b1 = &mgl[level+1].b, *e1 = &mgl[level+1].x; // coarse level b and x
65 
66  dCSRmat *A0 = &mgl[level].A; // fine level matrix
67  dCSRmat *A1 = &mgl[level+1].A; // coarse level matrix
68 
69  const INT m0 = A0->row, m1 = A1->row;
70 
71  INT *ordering = mgl[level].cfmark.val; // smoother ordering
72  ILU_data *LU_level = &mgl[level].LU; // fine level ILU decomposition
73  REAL *r = mgl[level].w.val; // work array for residual
74  REAL *r1 = mgl[level+1].w.val+m1; // work array for residual
75 
76  // Schwarz parameters
77  Schwarz_param swzparam;
78  if ( param->Schwarz_levels > 0 ) {
79  swzparam.Schwarz_blksolver = param->Schwarz_blksolver;
80  }
81 
82 #if DEBUG_MODE > 0
83  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
84  printf("### DEBUG: n=%d, nnz=%d\n", mgl[0].A.row, mgl[0].A.nnz);
85 #endif
86 
87  if ( prtlvl >= PRINT_MOST )
88  printf("AMLI level %d, smoother %d.\n", level, smoother);
89 
90  if ( level < mgl[level].num_levels-1 ) {
91 
92  // pre smoothing
93  if ( level < mgl[level].ILU_levels ) {
94 
95  fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
96 
97  }
98 
99  else if ( level < mgl->Schwarz_levels ) {
100 
101  switch (mgl[level].Schwarz.Schwarz_type) {
102  case SCHWARZ_SYMMETRIC:
103  fasp_dcsr_Schwarz_forward_smoother(&mgl[level].Schwarz, &swzparam,
104  &mgl[level].x, &mgl[level].b);
105  fasp_dcsr_Schwarz_backward_smoother(&mgl[level].Schwarz, &swzparam,
106  &mgl[level].x, &mgl[level].b);
107  break;
108  default:
109  fasp_dcsr_Schwarz_forward_smoother(&mgl[level].Schwarz, &swzparam,
110  &mgl[level].x, &mgl[level].b);
111  break;
112  }
113  }
114 
115  else {
116  fasp_dcsr_presmoothing(smoother,A0,b0,e0,param->presmooth_iter,
117  0,m0-1,1,relax,ndeg,smooth_order,ordering);
118  }
119 
120  // form residual r = b - A x
121  fasp_array_cp(m0,b0->val,r);
122  fasp_blas_dcsr_aAxpy(-1.0,A0,e0->val,r);
123 
124  // restriction r1 = R*r0
125  switch (amg_type) {
126  case UA_AMG:
127  fasp_blas_dcsr_mxv_agg(&mgl[level].R, r, b1->val); break;
128  default:
129  fasp_blas_dcsr_mxv(&mgl[level].R, r, b1->val); break;
130  }
131 
132  // coarse grid correction
133  {
134  INT i;
135 
136  fasp_array_cp(m1,b1->val,r1);
137 
138  for ( i=1; i<=degree; i++ ) {
139  fasp_dvec_set(m1,e1,0.0);
140  fasp_solver_amli(mgl, param, level+1);
141 
142  // b1 = (coef[degree-i]/coef[degree])*r1 + A1*e1;
143  // First, compute b1 = A1*e1
144  fasp_blas_dcsr_mxv(A1, e1->val, b1->val);
145  // Then, compute b1 = b1 + (coef[degree-i]/coef[degree])*r1
146  fasp_blas_array_axpy(m1, coef[degree-i]/coef[degree], r1, b1->val);
147  }
148 
149  fasp_dvec_set(m1,e1,0.0);
150  fasp_solver_amli(mgl, param, level+1);
151  }
152 
153  // find the optimal scaling factor alpha
154  fasp_blas_array_ax(m1, coef[degree], e1->val);
155  if ( param->coarse_scaling == ON ) {
156  alpha = fasp_blas_array_dotprod(m1, e1->val, r1)
157  / fasp_blas_dcsr_vmv(A1, e1->val, e1->val);
158  alpha = MIN(alpha, 1.0);
159  }
160 
161  // prolongation e0 = e0 + alpha * P * e1
162  switch (amg_type) {
163  case UA_AMG:
164  fasp_blas_dcsr_aAxpy_agg(alpha, &mgl[level].P, e1->val, e0->val);
165  break;
166  default:
167  fasp_blas_dcsr_aAxpy(alpha, &mgl[level].P, e1->val, e0->val);
168  break;
169  }
170 
171  // post smoothing
172  if ( level < mgl[level].ILU_levels ) {
173 
174  fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
175 
176  }
177 
178  else if (level<mgl->Schwarz_levels) {
179 
180  switch (mgl[level].Schwarz.Schwarz_type) {
181  case SCHWARZ_SYMMETRIC:
182  fasp_dcsr_Schwarz_backward_smoother(&mgl[level].Schwarz, &swzparam,
183  &mgl[level].x, &mgl[level].b);
184  fasp_dcsr_Schwarz_forward_smoother(&mgl[level].Schwarz, &swzparam,
185  &mgl[level].x, &mgl[level].b);
186  break;
187  default:
188  fasp_dcsr_Schwarz_backward_smoother(&mgl[level].Schwarz, &swzparam,
189  &mgl[level].x, &mgl[level].b);
190  break;
191  }
192  }
193 
194  else {
195  fasp_dcsr_postsmoothing(smoother,A0,b0,e0,param->postsmooth_iter,
196  0,m0-1,-1,relax,ndeg,smooth_order,ordering);
197  }
198 
199  }
200 
201  else { // coarsest level solver
202 
203  switch (coarse_solver) {
204 
205 #if WITH_SuperLU
206  case SOLVER_SUPERLU:
207  /* use SuperLU direct solver on the coarsest level */
208  fasp_solver_superlu(A0, b0, e0, 0);
209  break;
210 #endif
211 
212 #if WITH_UMFPACK
213  case SOLVER_UMFPACK:
214  /* use UMFPACK direct solver on the coarsest level */
215  fasp_umfpack_solve(A0, b0, e0, mgl[level].Numeric, 0);
216  break;
217 #endif
218 
219 #if WITH_MUMPS
220  case SOLVER_MUMPS:
221  /* use MUMPS direct solver on the coarsest level */
222  mgl[level].mumps.job = 2;
223  fasp_solver_mumps_steps(A0, b0, e0, &mgl[level].mumps);
224  break;
225 #endif
226 
227 #if WITH_PARDISO
228  case SOLVER_PARDISO: {
229  // user Intel MKL PARDISO direct solver on the coarsest level
230  fasp_pardiso_solve(A0, b0, e0, &mgl[level].pdata, 0);
231  break;
232  }
233 #endif
234 
235  default:
236  /* use iterative solver on the coarsest level */
237  fasp_coarse_itsolver(A0, b0, e0, tol, prtlvl);
238 
239  }
240 
241  }
242 
243 #if DEBUG_MODE > 0
244  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
245 #endif
246 }
247 
270  AMG_param *param,
271  INT level,
272  INT num_levels)
273 {
274  const SHORT amg_type=param->AMG_type;
275  const SHORT prtlvl = param->print_level;
276  const SHORT smoother = param->smoother;
277  const SHORT smooth_order = param->smooth_order;
278  const SHORT coarse_solver = param->coarse_solver;
279  const REAL relax = param->relaxation;
280  const REAL tol = param->tol*1e-4;
281  const SHORT ndeg = param->polynomial_degree;
282 
283  dvector *b0 = &mgl[level].b, *e0 = &mgl[level].x; // fine level b and x
284  dvector *b1 = &mgl[level+1].b, *e1 = &mgl[level+1].x; // coarse level b and x
285 
286  dCSRmat *A0 = &mgl[level].A; // fine level matrix
287  dCSRmat *A1 = &mgl[level+1].A; // coarse level matrix
288 
289  const INT m0 = A0->row, m1 = A1->row;
290 
291  INT *ordering = mgl[level].cfmark.val; // smoother ordering
292  ILU_data *LU_level = &mgl[level].LU; // fine level ILU decomposition
293  REAL *r = mgl[level].w.val; // work array for residual
294 
295  dvector uH; // for coarse level correction
296  uH.row = m1; uH.val = mgl[level+1].w.val + m1;
297 
298  // Schwarz parameters
299  Schwarz_param swzparam;
300  if ( param->Schwarz_levels > 0 ) {
301  swzparam.Schwarz_blksolver = param->Schwarz_blksolver;
302  }
303 
304 #if DEBUG_MODE > 0
305  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
306  printf("### DEBUG: n=%d, nnz=%d\n", mgl[0].A.row, mgl[0].A.nnz);
307 #endif
308 
309  if ( prtlvl >= PRINT_MOST )
310  printf("Nonlinear AMLI level %d, smoother %d.\n", num_levels, smoother);
311 
312  if ( level < num_levels-1 ) {
313 
314  // pre smoothing
315  if ( level < mgl[level].ILU_levels ) {
316 
317  fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
318 
319  }
320 
321  else if ( level < mgl->Schwarz_levels ) {
322 
323  switch (mgl[level].Schwarz.Schwarz_type) {
324  case SCHWARZ_SYMMETRIC:
325  fasp_dcsr_Schwarz_forward_smoother(&mgl[level].Schwarz, &swzparam,
326  &mgl[level].x, &mgl[level].b);
327  fasp_dcsr_Schwarz_backward_smoother(&mgl[level].Schwarz, &swzparam,
328  &mgl[level].x, &mgl[level].b);
329  break;
330  default:
331  fasp_dcsr_Schwarz_forward_smoother(&mgl[level].Schwarz, &swzparam,
332  &mgl[level].x, &mgl[level].b);
333  break;
334  }
335  }
336 
337  else {
338 
339  fasp_dcsr_presmoothing(smoother,A0,b0,e0,param->presmooth_iter,
340  0,m0-1,1,relax,ndeg,smooth_order,ordering);
341 
342  }
343 
344  // form residual r = b - A x
345  fasp_array_cp(m0,b0->val,r);
346  fasp_blas_dcsr_aAxpy(-1.0,A0,e0->val,r);
347 
348  // restriction r1 = R*r0
349  switch (amg_type) {
350  case UA_AMG:
351  fasp_blas_dcsr_mxv_agg(&mgl[level].R, r, b1->val);
352  break;
353  default:
354  fasp_blas_dcsr_mxv(&mgl[level].R, r, b1->val);
355  break;
356  }
357 
358  // call nonlinear AMLI-cycle recursively
359  {
360  fasp_dvec_set(m1,e1,0.0);
361 
362  // V-cycle will be enforced when needed !!!
363  if ( mgl[level+1].cycle_type <= 1 ) {
364 
365  fasp_solver_nl_amli(&mgl[level+1], param, 0, num_levels-1);
366 
367  }
368 
369  else { // recursively call preconditioned Krylov method on coarse grid
370 
371  precond_data pcdata;
372 
373  fasp_param_amg_to_prec(&pcdata, param);
374  pcdata.maxit = 1;
375  pcdata.max_levels = num_levels-1;
376  pcdata.mgl_data = &mgl[level+1];
377 
378  precond pc;
379  pc.data = &pcdata;
381 
382  fasp_array_cp (m1, e1->val, uH.val);
383 
384  switch (param->nl_amli_krylov_type) {
385 
386  case SOLVER_GCG: // Use GCG
387  fasp_krylov_cycle_dcsr_pgcg(A1,b1,&uH,&pc);
388  break;
389 
390  default: // Use GCR
391  fasp_krylov_cycle_dcsr_pgcr(A1,b1,&uH,&pc);
392  break;
393  }
394 
395  fasp_array_cp (m1, uH.val, e1->val);
396  }
397 
398  }
399 
400  // prolongation e0 = e0 + P*e1
401  switch (amg_type) {
402  case UA_AMG:
403  fasp_blas_dcsr_aAxpy_agg(1.0, &mgl[level].P, e1->val, e0->val);
404  break;
405  default:
406  fasp_blas_dcsr_aAxpy(1.0, &mgl[level].P, e1->val, e0->val);
407  break;
408  }
409 
410  // post smoothing
411  if ( level < mgl[level].ILU_levels ) {
412 
413  fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
414 
415  }
416  else if ( level < mgl->Schwarz_levels ) {
417 
418  switch (mgl[level].Schwarz.Schwarz_type) {
419  case SCHWARZ_SYMMETRIC:
420  fasp_dcsr_Schwarz_backward_smoother(&mgl[level].Schwarz, &swzparam,
421  &mgl[level].x, &mgl[level].b);
422  fasp_dcsr_Schwarz_forward_smoother(&mgl[level].Schwarz, &swzparam,
423  &mgl[level].x, &mgl[level].b);
424  break;
425  default:
426  fasp_dcsr_Schwarz_backward_smoother(&mgl[level].Schwarz, &swzparam,
427  &mgl[level].x, &mgl[level].b);
428  break;
429  }
430 
431  }
432 
433  else {
434  fasp_dcsr_postsmoothing(smoother,A0,b0,e0,param->postsmooth_iter,
435  0,m0-1,-1,relax,ndeg,smooth_order,ordering);
436  }
437 
438  }
439 
440  else { // coarsest level solver
441 
442  switch (coarse_solver) {
443 
444 #if WITH_SuperLU
445  case SOLVER_SUPERLU:
446  /* use SuperLU direct solver on the coarsest level */
447  fasp_solver_superlu(A0, b0, e0, 0);
448  break;
449 #endif
450 
451 #if WITH_UMFPACK
452  case SOLVER_UMFPACK:
453  /* use UMFPACK direct solver on the coarsest level */
454  fasp_umfpack_solve(A0, b0, e0, mgl[level].Numeric, 0);
455  break;
456 #endif
457 
458 #if WITH_MUMPS
459  case SOLVER_MUMPS:
460  /* use MUMPS direct solver on the coarsest level */
461  mgl[level].mumps.job = 2;
462  fasp_solver_mumps_steps(A0, b0, e0, &mgl[level].mumps);
463  break;
464 #endif
465 
466 #if WITH_PARDISO
467  case SOLVER_PARDISO: {
468  // user Intel MKL PARDISO direct solver on the coarsest level
469  fasp_pardiso_solve(A0, b0, e0, &mgl[level].pdata, 0);
470  break;
471  }
472 #endif
473 
474  default:
475  /* use iterative solver on the coarsest level */
476  fasp_coarse_itsolver(A0, b0, e0, tol, prtlvl);
477 
478  }
479 
480  }
481 
482 #if DEBUG_MODE > 0
483  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
484 #endif
485 }
486 
509  AMG_param *param,
510  INT level,
511  INT num_levels)
512 {
513  const SHORT prtlvl = param->print_level;
514  const SHORT smoother = param->smoother;
515  const SHORT coarse_solver = param->coarse_solver;
516  const REAL relax = param->relaxation;
517  const REAL tol = param->tol;
518  INT i;
519 
520  dvector *b0 = &mgl[level].b, *e0 = &mgl[level].x; // fine level b and x
521  dvector *b1 = &mgl[level+1].b, *e1 = &mgl[level+1].x; // coarse level b and x
522 
523  dBSRmat *A0 = &mgl[level].A; // fine level matrix
524  dBSRmat *A1 = &mgl[level+1].A; // coarse level matrix
525  const INT m0 = A0->ROW*A0->nb, m1 = A1->ROW*A1->nb;
526 
527  ILU_data *LU_level = &mgl[level].LU; // fine level ILU decomposition
528  REAL *r = mgl[level].w.val; // for residual
529 
530  dvector uH, bH; // for coarse level correction
531  uH.row = m1; uH.val = mgl[level+1].w.val + m1;
532  bH.row = m1; bH.val = mgl[level+1].w.val + 2*m1;
533 
534 #if DEBUG_MODE > 0
535  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
536  printf("### DEBUG: n=%d, nnz=%d\n", mgl[0].A.ROW, mgl[0].A.NNZ);
537 #endif
538 
539  if (prtlvl>=PRINT_MOST)
540  printf("Nonlinear AMLI level %d, pre-smoother %d.\n", level, smoother);
541 
542  if (level < num_levels-1) {
543 
544  // pre smoothing
545  if (level<param->ILU_levels) {
546  fasp_smoother_dbsr_ilu(A0, b0, e0, LU_level);
547  }
548  else {
549  SHORT steps = param->presmooth_iter;
550 
551  if (steps > 0) {
552  switch (smoother) {
553  case SMOOTHER_JACOBI:
554  for (i=0; i<steps; i++)
555  fasp_smoother_dbsr_jacobi (A0, b0, e0);
556  break;
557  case SMOOTHER_GS:
558  for (i=0; i<steps; i++)
559  fasp_smoother_dbsr_gs (A0, b0, e0, ASCEND, NULL);
560  break;
561  case SMOOTHER_SOR:
562  for (i=0; i<steps; i++)
563  fasp_smoother_dbsr_sor (A0, b0, e0, ASCEND, NULL,relax);
564  break;
565  default:
566  printf("### ERROR: Wrong smoother type %d!\n", smoother);
567  fasp_chkerr(ERROR_INPUT_PAR, __FUNCTION__);
568  }
569  }
570  }
571 
572  // form residual r = b - A x
573  fasp_array_cp(m0,b0->val,r);
574  fasp_blas_dbsr_aAxpy(-1.0,A0,e0->val,r);
575 
576  fasp_blas_dbsr_mxv(&mgl[level].R, r, b1->val);
577 
578  // call nonlinear AMLI-cycle recursively
579  {
580  fasp_dvec_set(m1,e1,0.0);
581 
582  // The coarsest problem is solved exactly.
583  // No need to call Krylov method on second coarsest level
584  if (level == num_levels-2) {
585  fasp_solver_nl_amli_bsr(&mgl[level+1], param, 0, num_levels-1);
586  }
587  else { // recursively call preconditioned Krylov method on coarse grid
588  precond_data_bsr pcdata;
589 
590  fasp_param_amg_to_prec_bsr (&pcdata, param);
591  pcdata.maxit = 1;
592  pcdata.max_levels = num_levels-1;
593  pcdata.mgl_data = &mgl[level+1];
594 
595  precond pc;
596  pc.data = &pcdata;
598 
599  fasp_array_cp (m1, b1->val, bH.val);
600  fasp_array_cp (m1, e1->val, uH.val);
601 
602  const INT maxit = param->amli_degree+1;
603 
604  fasp_solver_dbsr_pvfgmres(A1,&bH,&uH,&pc,param->tol,
605  maxit,MIN(maxit,30),1,PRINT_NONE);
606 
607  fasp_array_cp (m1, bH.val, b1->val);
608  fasp_array_cp (m1, uH.val, e1->val);
609  }
610 
611  }
612 
613  fasp_blas_dbsr_aAxpy(1.0, &mgl[level].P, e1->val, e0->val);
614 
615  // post smoothing
616  if (level < param->ILU_levels) {
617  fasp_smoother_dbsr_ilu(A0, b0, e0, LU_level);
618  }
619  else {
620  SHORT steps = param->postsmooth_iter;
621 
622  if (steps > 0) {
623  switch (smoother) {
624  case SMOOTHER_JACOBI:
625  for (i=0; i<steps; i++)
626  fasp_smoother_dbsr_jacobi (A0, b0, e0);
627  break;
628  case SMOOTHER_GS:
629  for (i=0; i<steps; i++)
630  fasp_smoother_dbsr_gs(A0, b0, e0, ASCEND, NULL);
631  break;
632  case SMOOTHER_SOR:
633  for (i=0; i<steps; i++)
634  fasp_smoother_dbsr_sor(A0, b0, e0, ASCEND, NULL,relax);
635  break;
636  default:
637  printf("### ERROR: Wrong smoother type %d!\n", smoother);
638  fasp_chkerr(ERROR_INPUT_PAR, __FUNCTION__);
639  }
640  }
641  }
642 
643  }
644 
645  else { // coarsest level solver
646 
647  switch (coarse_solver) {
648 
649 #if WITH_SuperLU
650  case SOLVER_SUPERLU:
651  /* use SuperLU direct solver on the coarsest level */
652  fasp_solver_superlu(&mgl[level].Ac, b0, e0, 0);
653  break;
654 #endif
655 
656 #if WITH_UMFPACK
657  case SOLVER_UMFPACK:
658  /* use UMFPACK direct solver on the coarsest level */
659  fasp_umfpack_solve(&mgl[level].Ac, b0, e0, mgl[level].Numeric, 0);
660  break;
661 #endif
662 
663 #if WITH_MUMPS
664  case SOLVER_MUMPS:
665  /* use MUMPS direct solver on the coarsest level */
666  mgl[level].mumps.job = 2;
667  fasp_solver_mumps_steps(&mgl[level].Ac, b0, e0, &mgl[level].mumps);
668  break;
669 #endif
670 
671 #if WITH_PARDISO
672  case SOLVER_PARDISO: {
673  // user Intel MKL PARDISO direct solver on the coarsest level
674  fasp_pardiso_solve(&mgl[level].Ac, b0, e0, &mgl[level].pdata, 0);
675  break;
676  }
677 #endif
678 
679  default:
680  /* use iterative solver on the coarsest level */
681  fasp_coarse_itsolver(&mgl[level].Ac, b0, e0, tol, prtlvl);
682 
683  }
684 
685  }
686 
687 #if DEBUG_MODE > 0
688  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
689 #endif
690 }
691 
706 void fasp_amg_amli_coef (const REAL lambda_max,
707  const REAL lambda_min,
708  const INT degree,
709  REAL *coef)
710 {
711  const REAL mu0 = 1.0/lambda_max, mu1 = 1.0/lambda_min;
712  const REAL c = (sqrt(mu0)+sqrt(mu1))*(sqrt(mu0)+sqrt(mu1));
713  const REAL a = (4*mu0*mu1)/(c);
714 
715  const REAL kappa = lambda_max/lambda_min; // condition number
716  const REAL delta = (sqrt(kappa) - 1.0)/(sqrt(kappa)+1.0);
717  const REAL b = delta*delta;
718 
719  if (degree == 0) {
720  coef[0] = 0.5*(mu0+mu1);
721  }
722 
723  else if (degree == 1) {
724  coef[0] = 0.5*c;
725  coef[1] = -1.0*mu0*mu1;
726  }
727 
728  else if (degree > 1) {
729  INT i;
730 
731  // allocate memory
732  REAL *work = (REAL *)fasp_mem_calloc(2*degree-1, sizeof(REAL));
733  REAL *coef_k, *coef_km1;
734  coef_k = work; coef_km1 = work+degree;
735 
736  // get q_k
737  fasp_amg_amli_coef(lambda_max, lambda_min, degree-1, coef_k);
738  // get q_km1
739  fasp_amg_amli_coef(lambda_max, lambda_min, degree-2, coef_km1);
740 
741  // get coef
742  coef[0] = a - b*coef_km1[0] + (1+b)*coef_k[0];
743 
744  for (i=1; i<degree-1; i++) {
745  coef[i] = -b*coef_km1[i] + (1+b)*coef_k[i] - a*coef_k[i-1];
746  }
747 
748  coef[degree-1] = (1+b)*coef_k[degree-1] - a*coef_k[degree-2];
749 
750  coef[degree] = -a*coef_k[degree-1];
751 
752  // clean memory
753  if (work) fasp_mem_free(work);
754  }
755 
756  else {
757  printf("### ERROR: Wrong AMLI degree %d!\n", degree);
758  fasp_chkerr(ERROR_INPUT_PAR, __FUNCTION__);
759  }
760 
761  return;
762 }
763 
764 /*---------------------------------*/
765 /*-- Private Functions --*/
766 /*---------------------------------*/
767 
784 static SHORT fasp_krylov_cycle_dcsr_pgcg (dCSRmat *A,
785  dvector *b,
786  dvector *u,
787  precond *pc)
788 {
789  REAL absres, relres, normb;
790  REAL alpha1, alpha2, gamma1, gamma2, rho1, rho2, beta1, beta2, beta3, beta4;
791  REAL *work, *r, *x1, *v1, *v2;
792 
793  INT m=A->row;
794  REAL *x = u->val;
795 
796  // allocate temp memory
797  work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
798  r = work; x1 = r + m; v1 = r + 2*m; v2 = r + 3*m;
799 
800  normb = fasp_blas_array_norm2(m, b->val);
801 
802  fasp_array_cp(m, b->val, r);
803 
804  // Preconditioning
805  if (pc != NULL)
806  pc->fct(r, x, pc->data);
807  else
808  fasp_array_cp(m, r, x);
809 
810  // v1 = A*p
811  fasp_blas_dcsr_mxv(A, x, v1);
812 
813  // rho1 = (p,v1)
814  rho1 = fasp_blas_array_dotprod (m, x, v1);
815 
816  // alpha1 = (p, r)
817  alpha1 = fasp_blas_array_dotprod (m, x, r);
818 
819  beta1 = alpha1/rho1;
820 
821  // r = r - beta1 *v1
822  fasp_blas_array_axpy(m, -beta1, v1, r);
823 
824  // norm(r)
825  absres = fasp_blas_array_norm2(m, r);
826 
827  // compute relative residual
828  relres = absres/normb;
829 
830  // if relres reaches tol(0.2), pgcg will stop,
831  // otherwise, another one pgcg iteration will do.
832  if (relres < 0.2) {
833  fasp_blas_array_ax(m, beta1, x);
834  fasp_mem_free(work);
835  return FASP_SUCCESS;
836  }
837 
838  // Preconditioning
839  if (pc != NULL)
840  pc->fct(r, x1, pc->data);
841  else
842  fasp_array_cp(m, r, x1);
843 
844  //v2 = A*p
845  fasp_blas_dcsr_mxv(A, x1, v2);
846 
847  //gamma0 = (x1,v1)
848  gamma1 = fasp_blas_array_dotprod (m, x1, v1);
849 
850  //alpha2 = (x1,r)
851  alpha2 = fasp_blas_array_dotprod(m, x1, r);
852 
853  //rho2 = (x1,v2)
854  rho2 = fasp_blas_array_dotprod(m, x1, v2);
855 
856  gamma2 = gamma1;
857 
858  beta2 = rho2 - gamma1*gamma2/rho1;
859  beta3 = (alpha1 - gamma2*alpha2/beta2)/rho1;
860  beta4 = alpha2/beta2;
861 
862  fasp_blas_array_ax(m, beta3, x);
863 
864  fasp_blas_array_axpy(m, beta4, x1, x);
865 
866  // free
867  fasp_mem_free(work);
868  return FASP_SUCCESS;
869 }
870 
887 static SHORT fasp_krylov_cycle_dcsr_pgcr (dCSRmat *A,
888  dvector *b,
889  dvector *u,
890  precond *pc)
891 {
892  REAL absres = BIGREAL;
893  REAL relres = BIGREAL, normb = BIGREAL;
894  REAL alpha, alpha1, alpha2, alpha3, alpha4, beta, gamma, rho1, rho2;
895 
896  INT m=A->row;
897  REAL *x = u->val;
898 
899  // allocate temp memory
900  REAL *work, *r, *x1, *v1, *v2;
901  work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
902  r = work; x1 = r + m; v1 = r + 2*m; v2 = r + 3*m;
903 
904  normb=fasp_blas_array_norm2(m, b->val);
905  fasp_array_cp(m, b->val, r);
906 
907  // Preconditioning
908  if (pc != NULL)
909  pc->fct(r, x, pc->data);
910  else
911  fasp_array_cp(m, r, x);
912 
913  // v1 = A*x
914  fasp_blas_dcsr_mxv(A, x, v1);
915  // rho1 = (v1,v1)
916  rho1 = fasp_blas_array_dotprod (m, v1, v1);
917  // alpha1 = (r, v1)
918  alpha1 = fasp_blas_array_dotprod (m, v1, r);
919 
920  alpha = alpha1/rho1;
921 
922  // r = r - alpha *v1
923  fasp_blas_array_axpy(m, -alpha, v1, r);
924 
925  // norm(r)
926  absres = fasp_blas_array_norm2(m, r);
927 
928  // compute relative residual
929  relres = absres/normb;
930 
931  // if relres reaches tol(0.2), pgcr will stop,
932  // otherwise, another one pgcr iteration will do.
933  if (relres < 0.2) {
934  fasp_blas_array_ax(m, alpha, x);
935  fasp_mem_free(work);
936  return FASP_SUCCESS;
937  }
938 
939  // Preconditioning
940  if (pc != NULL)
941  pc->fct(r, x1, pc->data);
942  else
943  fasp_array_cp(m, r, x1);
944 
945  //v2 = A*x1
946  fasp_blas_dcsr_mxv(A, x1, v2);
947 
948  //gamma = (v1,v2)
949  gamma = fasp_blas_array_dotprod (m, v1, v2);
950  //beta = (v2,v2)
951  beta = fasp_blas_array_dotprod(m, v2, v2);
952  //alpha2 = (r,v2)
953  alpha2 = fasp_blas_array_dotprod(m, r, v2);
954 
955  rho2 = beta - gamma*gamma/rho1;
956 
957  alpha3 = alpha1/rho1 - gamma*alpha2/(rho1*rho2);
958 
959  alpha4 = alpha2/rho2;
960 
961  // x = alpha3*x + alpha4*x1
962  fasp_blas_array_ax(m, alpha3, x);
963  fasp_blas_array_axpy(m, alpha4, x1, x);
964 
965  // free
966  fasp_mem_free(work);
967  return FASP_SUCCESS;
968 }
969 
970 /*---------------------------------*/
971 /*-- End of File --*/
972 /*---------------------------------*/
int fasp_solver_superlu(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Au=b by SuperLU.
#define SCHWARZ_SYMMETRIC
Definition: fasp_const.h:157
dvector x
pointer to the iterative solution at level level_num
Definition: fasp_block.h:219
#define SMOOTHER_GS
Definition: fasp_const.h:184
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: message.c:199
void fasp_dcsr_Schwarz_forward_smoother(Schwarz_data *Schwarz, Schwarz_param *param, dvector *x, dvector *b)
Schwarz smoother: forward sweep.
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.
SHORT AMG_type
type of AMG method
Definition: fasp.h:586
Parameters for AMG solver.
Definition: fasp.h:583
SHORT postsmooth_iter
number of postsmoothers
Definition: fasp.h:619
SHORT smoother
smoother type
Definition: fasp.h:610
#define ERROR_INPUT_PAR
Definition: fasp_const.h:31
SHORT smooth_order
smoother order
Definition: fasp.h:613
#define REAL
Definition: fasp.h:67
SHORT amli_degree
degree of the polynomial used by AMLI cycle
Definition: fasp.h:634
void fasp_smoother_dcsr_ilu(dCSRmat *A, dvector *b, dvector *x, void *data)
ILU method as a smoother.
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: amlirecur.c:706
SHORT presmooth_iter
number of presmoothers
Definition: fasp.h:616
dBSRmat A
pointer to the matrix at level level_num
Definition: fasp_block.h:207
SHORT polynomial_degree
degree of the polynomial smoother
Definition: fasp.h:625
INT maxit
max number of iterations of AMG preconditioner
Definition: fasp.h:804
INT nb
dimension of each sub-block
Definition: fasp_block.h:56
#define SOLVER_SUPERLU
Definition: fasp_const.h:123
Mumps_data mumps
data for MUMPS
Definition: fasp_block.h:283
int fasp_solver_mumps_steps(dCSRmat *ptrA, dvector *b, dvector *u, Mumps_data *mumps)
Solve Ax=b by MUMPS in three steps.
dvector b
pointer to the right-hand side at level level_num
Definition: fasp.h:744
#define SOLVER_GCG
Definition: fasp_const.h:109
void fasp_dcsr_Schwarz_backward_smoother(Schwarz_data *Schwarz, Schwarz_param *param, dvector *x, dvector *b)
Schwarz smoother: backward sweep.
REAL * val
nonzero entries of A
Definition: fasp.h:166
#define PRINT_MOST
Definition: fasp_const.h:83
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:44
REAL * val
actual vector entries
Definition: fasp.h:348
ILU_data LU
ILU matrix for ILU smoother.
Definition: fasp_block.h:258
void fasp_blas_dcsr_aAxpy_agg(const REAL alpha, dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y (the entries of A are all ones)
Definition: blas_csr.c:593
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1005
Vector with n entries of REAL type.
Definition: fasp.h:342
void fasp_smoother_dbsr_jacobi(dBSRmat *A, dvector *b, dvector *u)
Jacobi relaxation.
Definition: smoother_bsr.c:33
dvector w
temporary work space
Definition: fasp_block.h:280
#define SMOOTHER_SOR
Definition: fasp_const.h:187
void fasp_precond_nl_amli(REAL *r, REAL *z, void *data)
Nonlinear AMLI AMG preconditioner.
Definition: precond_csr.c:502
INT max_levels
max number of AMG levels
Definition: fasp_block.h:323
#define INT
Definition: fasp.h:64
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
Data for ILU setup.
Definition: fasp.h:400
void fasp_solver_nl_amli_bsr(AMG_data_bsr *mgl, AMG_param *param, INT level, INT num_levels)
Solve Ax=b with recursive nonlinear AMLI-cycle.
Definition: amlirecur.c:508
Preconditioner data and action.
Definition: fasp.h:999
#define ASCEND
Definition: fasp_const.h:231
Data for multigrid levels. (BSR format)
Definition: fasp_block.h:198
INT Schwarz_levels
number of levels use Schwarz smoother
Definition: fasp.h:700
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void fasp_precond_dbsr_nl_amli(REAL *r, REAL *z, void *data)
Nonlinear AMLI-cycle AMG preconditioner.
Definition: precond_bsr.c:607
ILU_data LU
ILU matrix for ILU smoother.
Definition: fasp.h:764
REAL * val
Definition: fasp_block.h:67
void * data
data for preconditioner, void pointer
Definition: fasp.h:1002
INT NNZ
number of nonzero sub-blocks in matrix A, NNZ
Definition: fasp_block.h:53
void fasp_param_amg_to_prec(precond_data *pcdata, AMG_param *amgparam)
Set precond_data with AMG_param.
Definition: parameters.c:666
dvector x
pointer to the iterative solution at level level_num
Definition: fasp.h:747
Data passed to the preconditioners.
Definition: fasp_block.h:311
INT nnz
number of nonzero entries
Definition: fasp.h:157
dvector b
pointer to the right-hand side at level level_num
Definition: fasp_block.h:216
SHORT coarse_solver
coarse solver type
Definition: fasp.h:628
#define SOLVER_PARDISO
Definition: fasp_const.h:126
Main header file for FASP.
REAL * amli_coef
coefficients of the polynomial used by AMLI cycle
Definition: fasp.h:637
#define SOLVER_MUMPS
Definition: fasp_const.h:125
SHORT nl_amli_krylov_type
type of Krylov method used by Nonlinear AMLI cycle
Definition: fasp.h:640
#define SMOOTHER_JACOBI
Definition of standard smoother types.
Definition: fasp_const.h:183
INT * val
actual vector entries
Definition: fasp.h:362
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
REAL relaxation
relaxation parameter for SOR smoother
Definition: fasp.h:622
INT row
row number of matrix A, m
Definition: fasp.h:151
INT Schwarz_blksolver
type of Schwarz block solver
Definition: fasp.h:449
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
INT Schwarz_blksolver
type of Schwarz block solver
Definition: fasp.h:712
void fasp_blas_array_axpy(const INT n, const REAL a, REAL *x, REAL *y)
y = a*x + y
Definition: blas_array.c:87
ivector cfmark
pointer to the CF marker at level level_num
Definition: fasp.h:758
AMG_data_bsr * mgl_data
AMG preconditioner data.
Definition: fasp_block.h:368
void fasp_smoother_dbsr_sor(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark, REAL weight)
SOR relaxation.
INT row
number of rows
Definition: fasp.h:345
SHORT max_levels
max number of AMG levels
Definition: fasp.h:807
Parameters for Schwarz method.
Definition: fasp.h:434
AMG_data * mgl_data
AMG preconditioner data.
Definition: fasp.h:855
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define MIN(a, b)
Definition: fasp.h:73
Mumps_data mumps
data for MUMPS
Definition: fasp.h:784
void fasp_blas_dbsr_mxv(dBSRmat *A, REAL *x, REAL *y)
Compute y := A*x.
Definition: blas_bsr.c:895
REAL fasp_blas_array_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies (x,y)
Definition: blas_array.c:267
void fasp_param_amg_to_prec_bsr(precond_data_bsr *pcdata, AMG_param *amgparam)
Set precond_data_bsr with AMG_param.
Definition: parameters.c:733
REAL fasp_blas_array_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: blas_array.c:347
dCSRmat A
pointer to the matrix at level level_num
Definition: fasp.h:735
SHORT print_level
print level for AMG
Definition: fasp.h:589
#define ON
Definition of switch.
Definition: fasp_const.h:73
#define BIGREAL
Some global constants.
Definition: fasp_const.h:237
void fasp_smoother_dbsr_gs(dBSRmat *A, dvector *b, dvector *u, INT order, INT *mark)
Gauss-Seidel relaxation.
Definition: smoother_bsr.c:411
INT fasp_solver_dbsr_pvfgmres(dBSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT stop_type, const SHORT prtlvl)
Solve "Ax=b" using PFGMRES(right preconditioned) iterative method in which the restart parameter can ...
Definition: pvfgmres.c:382
Data for AMG solvers.
Definition: fasp.h:722
void fasp_blas_array_ax(const INT n, const REAL a, REAL *x)
x = a*x
Definition: blas_array.c:35
void fasp_blas_dcsr_aAxpy(const REAL alpha, dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_csr.c:479
void fasp_dvec_set(INT n, dvector *x, REAL val)
Initialize dvector x[i]=val for i=0:n-1.
Definition: vec.c:235
#define SOLVER_UMFPACK
Definition: fasp_const.h:124
Data passed to the preconditioners.
Definition: fasp.h:795
REAL tol
stopping tolerance for AMG solver
Definition: fasp.h:595
INT maxit
max number of iterations of AMG preconditioner
Definition: fasp_block.h:320
INT ROW
number of rows of sub-blocks in matrix A, M
Definition: fasp_block.h:47
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
void fasp_blas_dcsr_mxv_agg(dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x, where the entries of A are all ones.
Definition: blas_csr.c:423
void fasp_solver_nl_amli(AMG_data *mgl, AMG_param *param, INT level, INT num_levels)
Solve Ax=b with recursive nonlinear AMLI-cycle.
Definition: amlirecur.c:269
INT job
work for MUMPS
Definition: fasp.h:467
void fasp_solver_amli(AMG_data *mgl, AMG_param *param, INT level)
Solve Ax=b with recursive AMLI-cycle.
Definition: amlirecur.c:45
dvector w
Temporary work space.
Definition: fasp.h:781
void fasp_blas_dbsr_aAxpy(const REAL alpha, dBSRmat *A, REAL *x, REAL *y)
Compute y := alpha*A*x + y.
Definition: blas_bsr.c:337
SHORT coarse_scaling
switch of scaling of the coarse grid correction
Definition: fasp.h:631
#define UA_AMG
Definition: fasp_const.h:164
REAL fasp_blas_dcsr_vmv(dCSRmat *A, REAL *x, REAL *y)
vector-Matrix-vector multiplication alpha = y'*A*x
Definition: blas_csr.c:704
void fasp_blas_dcsr_mxv(dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_csr.c:225