Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
fmgcycle.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <time.h>
8 
9 #include "fasp.h"
10 #include "fasp_functs.h"
12 
13 #include "mg_util.inl"
14 
15 /*---------------------------------*/
16 /*-- Public Functions --*/
17 /*---------------------------------*/
18 
37  AMG_param *param)
38 {
39  const SHORT maxit = 3; // Max allowed V-cycles in each level
40  const SHORT amg_type = param->AMG_type;
41  const SHORT prtlvl = param->print_level;
42  const SHORT nl = mgl[0].num_levels;
43  const SHORT smoother = param->smoother;
44  const SHORT smooth_order = param->smooth_order;
45  const SHORT coarse_solver = param->coarse_solver;
46 
47  const REAL relax = param->relaxation;
48  const SHORT ndeg = param->polynomial_degree;
49  const REAL tol = param->tol*1e-4;
50 
51  // local variables
52  INT l, i, lvl, num_cycle;
53  REAL alpha = 1.0, relerr;
54 
55  // Schwarz parameters
56  Schwarz_param swzparam;
57  if ( param->Schwarz_levels > 0 ) {
58  swzparam.Schwarz_blksolver = param->Schwarz_blksolver;
59  }
60 
61 #if DEBUG_MODE > 0
62  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
63  printf("### DEBUG: n=%d, nnz=%d\n", mgl[0].A.row, mgl[0].A.nnz);
64 #endif
65 
66  if ( prtlvl >= PRINT_MOST )
67  printf("FMG_level = %d, ILU_level = %d\n", nl, param->ILU_levels);
68 
69  // restriction r1 = R*r0
70  switch (amg_type) {
71 
72  case UA_AMG:
73  for (l=0;l<nl-1;l++)
74  fasp_blas_dcsr_mxv_agg(&mgl[l].R, mgl[l].b.val, mgl[l+1].b.val);
75  break;
76 
77  default:
78  for (l=0;l<nl-1;l++)
79  fasp_blas_dcsr_mxv(&mgl[l].R, mgl[l].b.val, mgl[l+1].b.val);
80  break;
81 
82  }
83 
84  fasp_dvec_set(mgl[l].A.row, &mgl[l].x, 0.0); // initial guess
85 
86  // If only one level, just direct solver
87  if ( nl==1 ) {
88 
89  switch (coarse_solver) {
90 
91 #if WITH_SuperLU
92  case SOLVER_SUPERLU:
93  /* use SuperLU direct solver on the coarsest level */
94  fasp_solver_superlu(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, 0);
95  break;
96 #endif
97 
98 #if WITH_UMFPACK
99  case SOLVER_UMFPACK:
100  /* use UMFPACK direct solver on the coarsest level */
101  fasp_umfpack_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, mgl[nl-1].Numeric, 0);
102  break;
103 #endif
104 
105 #if WITH_MUMPS
106  case SOLVER_MUMPS:
107  /* use MUMPS direct solver on the coarsest level */
108  mgl[nl-1].mumps.job = 2;
109  fasp_solver_mumps_steps(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].mumps);
110  break;
111 #endif
112 
113 #if WITH_PARDISO
114  case SOLVER_PARDISO: {
115  // user Intel MKL PARDISO direct solver on the coarsest level
116  fasp_pardiso_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].pdata, 0);
117  break;
118  }
119 #endif
120 
121  default:
122  /* use iterative solver on the coarest level */
123  fasp_coarse_itsolver(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, tol, prtlvl);
124 
125  }
126 
127  return;
128 
129  }
130 
131  for ( i=1; i<nl; i++ ) {
132 
133  // Coarse Space Solver:
134  switch (coarse_solver) {
135 
136 #if WITH_SuperLU
137  case SOLVER_SUPERLU:
138  /* use SuperLU direct solver on the coarsest level */
139  fasp_solver_superlu(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, 0);
140  break;
141 #endif
142 
143 #if WITH_UMFPACK
144  case SOLVER_UMFPACK:
145  /* use UMFPACK direct solver on the coarsest level */
146  fasp_umfpack_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, mgl[nl-1].Numeric, 0);
147  break;
148 #endif
149 
150 #if WITH_MUMPS
151  case SOLVER_MUMPS:
152  /* use MUMPS direct solver on the coarsest level */
153  mgl[nl-1].mumps.job = 2;
154  fasp_solver_mumps_steps(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].mumps);
155  break;
156 #endif
157 
158 #if WITH_PARDISO
159  case SOLVER_PARDISO: {
160  // user Intel MKL PARDISO direct solver on the coarsest level
161  fasp_pardiso_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].pdata, 0);
162  break;
163  }
164 #endif
165 
166  default:
167  /* use iterative solver on the coarest level */
168  fasp_coarse_itsolver(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, tol, prtlvl);
169 
170  }
171 
172  // Slash part: /-cycle
173  {
174  --l; // go back to finer level
175 
176  // find the optimal scaling factor alpha
177  if ( param->coarse_scaling == ON ) {
178  alpha = fasp_blas_array_dotprod(mgl[l+1].A.row, mgl[l+1].x.val, mgl[l+1].b.val)
179  / fasp_blas_dcsr_vmv(&mgl[l+1].A, mgl[l+1].x.val, mgl[l+1].x.val);
180  alpha = MIN(alpha, 1.0); // Add this for safty! --Chensong on 10/04/2014
181  }
182 
183  // prolongation u = u + alpha*P*e1
184  switch (amg_type) {
185  case UA_AMG:
186  fasp_blas_dcsr_aAxpy_agg(alpha, &mgl[l].P, mgl[l+1].x.val, mgl[l].x.val); break;
187  default:
188  fasp_blas_dcsr_aAxpy(alpha, &mgl[l].P, mgl[l+1].x.val, mgl[l].x.val); break;
189  }
190  }
191 
192  // initialzie rel error
193  num_cycle = 0; relerr = BIGREAL;
194 
195  while ( relerr > param->tol && num_cycle < maxit) {
196 
197  ++num_cycle;
198 
199  // form residual r = b - A x
200  fasp_array_cp(mgl[l].A.row, mgl[l].b.val, mgl[l].w.val);
201  fasp_blas_dcsr_aAxpy(-1.0,&mgl[l].A, mgl[l].x.val, mgl[l].w.val);
202  relerr = fasp_blas_dvec_norm2(&mgl[l].w) / fasp_blas_dvec_norm2(&mgl[l].b);
203 
204  // Forward Sweep
205  for ( lvl=0; lvl<i; lvl++ ) {
206 
207  // pre smoothing
208  if (l<param->ILU_levels) {
209  fasp_smoother_dcsr_ilu(&mgl[l].A, &mgl[l].b, &mgl[l].x, &mgl[l].LU);
210  }
211  else if (l<mgl->Schwarz_levels) {
212  switch (mgl[l].Schwarz.Schwarz_type) {
213  case SCHWARZ_SYMMETRIC:
214  fasp_dcsr_Schwarz_forward_smoother(&mgl[l].Schwarz, &swzparam,
215  &mgl[l].x, &mgl[l].b);
216  fasp_dcsr_Schwarz_backward_smoother(&mgl[l].Schwarz, &swzparam,
217  &mgl[l].x, &mgl[l].b);
218  break;
219  default:
220  fasp_dcsr_Schwarz_forward_smoother(&mgl[l].Schwarz, &swzparam,
221  &mgl[l].x, &mgl[l].b);
222  break;
223  }
224  }
225 
226  else {
227  fasp_dcsr_presmoothing(smoother,&mgl[l].A,&mgl[l].b,&mgl[l].x,param->presmooth_iter,
228  0,mgl[l].A.row-1,1,relax,ndeg,smooth_order,mgl[l].cfmark.val);
229  }
230 
231  // form residual r = b - A x
232  fasp_array_cp(mgl[l].A.row, mgl[l].b.val, mgl[l].w.val);
233  fasp_blas_dcsr_aAxpy(-1.0,&mgl[l].A, mgl[l].x.val, mgl[l].w.val);
234 
235  // restriction r1 = R*r0
236  switch (amg_type) {
237  case UA_AMG:
238  fasp_blas_dcsr_mxv_agg(&mgl[l].R, mgl[l].w.val, mgl[l+1].b.val);
239  break;
240  default:
241  fasp_blas_dcsr_mxv(&mgl[l].R, mgl[l].w.val, mgl[l+1].b.val);
242  break;
243  }
244 
245  ++l;
246 
247  // prepare for the next level
248  fasp_dvec_set(mgl[l].A.row, &mgl[l].x, 0.0);
249 
250  } // end for lvl
251 
252  // CoarseSpaceSolver:
253  switch (coarse_solver) {
254 
255 #if WITH_SuperLU
256  case SOLVER_SUPERLU:
257  /* use SuperLU direct solver on the coarsest level */
258  fasp_solver_superlu(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, 0);
259  break;
260 #endif
261 
262 #if WITH_UMFPACK
263  case SOLVER_UMFPACK:
264  /* use UMFPACK direct solver on the coarsest level */
265  fasp_umfpack_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, mgl[nl-1].Numeric, 0);
266  break;
267 #endif
268 
269 #if WITH_MUMPS
270  case SOLVER_MUMPS:
271  /* use MUMPS direct solver on the coarsest level */
272  mgl[nl-1].mumps.job = 2;
273  fasp_solver_mumps_steps(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].mumps);
274  break;
275 #endif
276 
277 #if WITH_PARDISO
278  case SOLVER_PARDISO: {
279  // user Intel MKL PARDISO direct solver on the coarsest level
280  fasp_pardiso_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].pdata, 0);
281  break;
282  }
283 #endif
284 
285  default:
286  /* use iterative solver on the coarest level */
287  fasp_coarse_itsolver(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, tol, prtlvl);
288 
289  }
290 
291  // Backward Sweep
292  for ( lvl=0; lvl<i; lvl++ ) {
293 
294  --l;
295 
296  // find the optimal scaling factor alpha
297  if ( param->coarse_scaling == ON ) {
298  alpha = fasp_blas_array_dotprod(mgl[l+1].A.row, mgl[l+1].x.val, mgl[l+1].b.val)
299  / fasp_blas_dcsr_vmv(&mgl[l+1].A, mgl[l+1].x.val, mgl[l+1].x.val);
300  alpha = MIN(alpha, 1.0); // Add this for safty! --Chensong on 10/04/2014
301  }
302 
303  // prolongation u = u + alpha*P*e1
304  switch (amg_type)
305  {
306  case UA_AMG:
307  fasp_blas_dcsr_aAxpy_agg(alpha, &mgl[l].P, mgl[l+1].x.val, mgl[l].x.val);
308  break;
309  default:
310  fasp_blas_dcsr_aAxpy(alpha, &mgl[l].P, mgl[l+1].x.val, mgl[l].x.val);
311  break;
312  }
313 
314  // post-smoothing
315  if (l<param->ILU_levels) {
316  fasp_smoother_dcsr_ilu(&mgl[l].A, &mgl[l].b, &mgl[l].x, &mgl[l].LU);
317  }
318  else if (l<mgl->Schwarz_levels) {
319  switch (mgl[l].Schwarz.Schwarz_type) {
320  case SCHWARZ_SYMMETRIC:
321  fasp_dcsr_Schwarz_backward_smoother(&mgl[l].Schwarz, &swzparam,
322  &mgl[l].x, &mgl[l].b);
323  fasp_dcsr_Schwarz_forward_smoother(&mgl[l].Schwarz, &swzparam,
324  &mgl[l].x, &mgl[l].b);
325  break;
326  default:
327  fasp_dcsr_Schwarz_backward_smoother(&mgl[l].Schwarz, &swzparam,
328  &mgl[l].x, &mgl[l].b);
329  break;
330  }
331  }
332 
333  else {
334  fasp_dcsr_postsmoothing(smoother,&mgl[l].A,&mgl[l].b,&mgl[l].x,param->postsmooth_iter,
335  0,mgl[l].A.row-1,-1,relax,ndeg,smooth_order,mgl[l].cfmark.val);
336  }
337 
338  } // end while
339 
340  } //end while
341 
342  } // end for
343 
344 #if DEBUG_MODE > 0
345  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
346 #endif
347 
348  return;
349 }
350 
351 /*---------------------------------*/
352 /*-- End of File --*/
353 /*---------------------------------*/
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
void fasp_solver_fmgcycle(AMG_data *mgl, AMG_param *param)
#include "forts_ns.h"
Definition: fmgcycle.c:36
void fasp_dcsr_Schwarz_forward_smoother(Schwarz_data *Schwarz, Schwarz_param *param, dvector *x, dvector *b)
Schwarz smoother: forward sweep.
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
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 num_levels
number of levels in use <= max_levels
Definition: fasp.h:730
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
void fasp_dcsr_Schwarz_backward_smoother(Schwarz_data *Schwarz, Schwarz_param *param, dvector *x, dvector *b)
Schwarz smoother: backward sweep.
#define PRINT_MOST
Definition: fasp_const.h:83
REAL * val
actual vector entries
Definition: fasp.h:348
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
#define INT
Definition: fasp.h:64
INT Schwarz_levels
number of levels use Schwarz smoother
Definition: fasp.h:700
REAL fasp_blas_dvec_norm2(dvector *x)
L2 norm of dvector x.
Definition: blas_vec.c:265
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
SHORT ILU_levels
number of levels use ILU smoother
Definition: fasp.h:682
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
INT Schwarz_blksolver
type of Schwarz block solver
Definition: fasp.h:449
INT Schwarz_blksolver
type of Schwarz block solver
Definition: fasp.h:712
ivector cfmark
pointer to the CF marker at level level_num
Definition: fasp.h:758
Parameters for Schwarz method.
Definition: fasp.h:434
#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
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
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
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
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
INT job
work for MUMPS
Definition: fasp.h:467
dvector w
Temporary work space.
Definition: fasp.h:781
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