Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
mgrecur.c
Go to the documentation of this file.
1 
8 #include <time.h>
9 
10 #include "fasp.h"
11 #include "fasp_functs.h"
12 #include "mg_util.inl"
13 
14 /*---------------------------------*/
15 /*-- Public Functions --*/
16 /*---------------------------------*/
17 
34  AMG_param *param,
35  INT level)
36 {
37  const SHORT prtlvl = param->print_level;
38  const SHORT smoother = param->smoother;
39  const SHORT cycle_type = param->cycle_type;
40  const SHORT coarse_solver = param->coarse_solver;
41  const SHORT smooth_order = param->smooth_order;
42  const REAL relax = param->relaxation;
43  const REAL tol = param->tol*1e-4;
44  const SHORT ndeg = param->polynomial_degree;
45 
46  dvector *b0 = &mgl[level].b, *e0 = &mgl[level].x; // fine level b and x
47  dvector *b1 = &mgl[level+1].b, *e1 = &mgl[level+1].x; // coarse level b and x
48 
49  dCSRmat *A0 = &mgl[level].A; // fine level matrix
50  dCSRmat *A1 = &mgl[level+1].A; // coarse level matrix
51  const INT m0 = A0->row, m1 = A1->row;
52 
53  ILU_data *LU_level = &mgl[level].LU; // fine level ILU decomposition
54  REAL *r = mgl[level].w.val; // for residual
55  INT *ordering = mgl[level].cfmark.val; // for smoother ordering
56 
57 #if DEBUG_MODE > 0
58  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
59  printf("### DEBUG: n=%d, nnz=%d\n", mgl[0].A.row, mgl[0].A.nnz);
60 #endif
61 
62  if ( prtlvl >= PRINT_MOST )
63  printf("AMG level %d, smoother %d.\n", level, smoother);
64 
65  if ( level < mgl[level].num_levels-1 ) {
66 
67  // pre smoothing
68  if ( level < mgl[level].ILU_levels ) {
69  fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
70  }
71  else {
72  fasp_dcsr_presmoothing(smoother,A0,b0,e0,param->presmooth_iter,
73  0,m0-1,1,relax,ndeg,smooth_order,ordering);
74  }
75 
76  // form residual r = b - A x
77  fasp_array_cp(m0,b0->val,r);
78  fasp_blas_dcsr_aAxpy(-1.0,A0,e0->val,r);
79 
80  // restriction r1 = R*r0
81  fasp_blas_dcsr_mxv(&mgl[level].R, r, b1->val);
82 
83  { // call MG recursively: type = 1 for V cycle, type = 2 for W cycle
84  SHORT i;
85  fasp_dvec_set(m1,e1,0.0);
86  for (i=0; i<cycle_type; ++i) fasp_solver_mgrecur (mgl, param, level+1);
87  }
88 
89  // prolongation e0 = e0 + P*e1
90  fasp_blas_dcsr_aAxpy(1.0, &mgl[level].P, e1->val, e0->val);
91 
92  // post smoothing
93  if ( level < mgl[level].ILU_levels ) {
94  fasp_smoother_dcsr_ilu(A0, b0, e0, LU_level);
95  }
96  else {
97  fasp_dcsr_postsmoothing(smoother,A0,b0,e0,param->postsmooth_iter,
98  0,m0-1,-1,relax,ndeg,smooth_order,ordering);
99  }
100 
101  }
102 
103  else { // coarsest level solver
104 
105  switch (coarse_solver) {
106 
107 #if WITH_SuperLU
108  /* use SuperLU direct solver on the coarsest level */
109  case SOLVER_SUPERLU:
110  fasp_solver_superlu(A0, b0, e0, 0);
111  break;
112 #endif
113 
114 #if WITH_UMFPACK
115  /* use UMFPACK direct solver on the coarsest level */
116  case SOLVER_UMFPACK:
117  //fasp_solver_umfpack(A0, b0, e0, 0);
118  fasp_umfpack_solve(A0, b0, e0, mgl[level].Numeric, 0);
119  break;
120 #endif
121 
122 #if WITH_MUMPS
123  /* use MUMPS direct solver on the coarsest level */
124  case SOLVER_MUMPS:
125  mgl[level].mumps.job = 2;
126  fasp_solver_mumps_steps(A0, b0, e0, &mgl[level].mumps);
127  break;
128 #endif
129 
130 #if WITH_PARDISO
131  case SOLVER_PARDISO: {
132  // user Intel MKL PARDISO direct solver on the coarsest level
133  fasp_pardiso_solve(A0, b0, e0, &mgl[level].pdata, 0);
134  break;
135  }
136 #endif
137 
138  /* use iterative solver on the coarsest level */
139  default:
140  fasp_coarse_itsolver(A0, b0, e0, tol, prtlvl);
141 
142  }
143 
144  }
145 
146 #if DEBUG_MODE > 0
147  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
148 #endif
149 }
150 
151 /*---------------------------------*/
152 /*-- End of File --*/
153 /*---------------------------------*/
int fasp_solver_superlu(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Au=b by SuperLU.
void fasp_solver_mgrecur(AMG_data *mgl, AMG_param *param, INT level)
Solve Ax=b with recursive multigrid K-cycle.
Definition: mgrecur.c:33
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
SHORT smooth_order
smoother order
Definition: fasp.h:613
#define REAL
Definition: fasp.h:67
void fasp_smoother_dcsr_ilu(dCSRmat *A, dvector *b, dvector *x, void *data)
ILU method as a smoother.
SHORT presmooth_iter
number of presmoothers
Definition: fasp.h:616
SHORT polynomial_degree
degree of the polynomial smoother
Definition: fasp.h:625
#define SOLVER_SUPERLU
Definition: fasp_const.h:123
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
REAL * val
nonzero entries of A
Definition: fasp.h:166
#define PRINT_MOST
Definition: fasp_const.h:83
REAL * val
actual vector entries
Definition: fasp.h:348
Vector with n entries of REAL type.
Definition: fasp.h:342
#define INT
Definition: fasp.h:64
Data for ILU setup.
Definition: fasp.h:400
ILU_data LU
ILU matrix for ILU smoother.
Definition: fasp.h:764
dvector x
pointer to the iterative solution at level level_num
Definition: fasp.h:747
INT nnz
number of nonzero entries
Definition: fasp.h:157
SHORT coarse_solver
coarse solver type
Definition: fasp.h:628
#define SOLVER_PARDISO
Definition: fasp_const.h:126
Main header file for FASP.
#define SOLVER_MUMPS
Definition: fasp_const.h:125
INT * val
actual vector entries
Definition: fasp.h:362
REAL relaxation
relaxation parameter for SOR smoother
Definition: fasp.h:622
INT row
row number of matrix A, m
Definition: fasp.h:151
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
ivector cfmark
pointer to the CF marker at level level_num
Definition: fasp.h:758
SHORT cycle_type
type of AMG cycle
Definition: fasp.h:604
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
Mumps_data mumps
data for MUMPS
Definition: fasp.h:784
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
Data for AMG solvers.
Definition: fasp.h:722
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
REAL tol
stopping tolerance for AMG solver
Definition: fasp.h:595
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
INT job
work for MUMPS
Definition: fasp.h:467
dvector w
Temporary work space.
Definition: fasp.h:781
void fasp_blas_dcsr_mxv(dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_csr.c:225