Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
amg_solve.c
Go to the documentation of this file.
1 
10 #include <time.h>
11 
12 #include "fasp.h"
13 #include "fasp_functs.h"
14 
15 #include "itsolver_util.inl"
16 
17 /*---------------------------------*/
18 /*-- Public Functions --*/
19 /*---------------------------------*/
20 
37  AMG_param *param)
38 {
39  dCSRmat *ptrA = &mgl[0].A;
40  dvector *b = &mgl[0].b, *x = &mgl[0].x, *r = &mgl[0].w;
41 
42  const SHORT prtlvl = param->print_level;
43  const INT MaxIt = param->maxit;
44  const REAL tol = param->tol;
45  const REAL sumb = fasp_blas_dvec_norm2(b); // L2norm(b)
46 
47  // local variables
48  REAL solve_start, solve_end;
49  REAL relres1 = BIGREAL, absres0 = sumb, absres, factor;
50  INT iter = 0;
51 
52 #if DEBUG_MODE > 0
53  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
54  printf("### DEBUG: nrow = %d, ncol = %d, nnz = %d\n",
55  mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
56 #endif
57 
58  fasp_gettime(&solve_start);
59 
60  // Print iteration information if needed
61  print_itinfo(prtlvl, STOP_REL_RES, iter, 1.0, sumb, 0.0);
62 
63  // MG solver here
64  while ( (++iter <= MaxIt) & (sumb > SMALLREAL) ) {
65 
66 #if TRUE
67  // Call one multigrid cycle -- non recursive version
68  fasp_solver_mgcycle(mgl, param);
69 #else
70  // Call one multigrid cycle -- recursive version
71  fasp_solver_mgrecur(mgl, param, 0);
72 #endif
73 
74  // Form residual r = b - A*x
75  fasp_dvec_cp(b, r);
76  fasp_blas_dcsr_aAxpy(-1.0, ptrA, x->val, r->val);
77 
78  // Compute norms of r and convergence factor
79  absres = fasp_blas_dvec_norm2(r); // residual ||r||
80  relres1 = absres/sumb; // relative residual ||r||/||b||
81  factor = absres/absres0; // contraction factor
82  absres0 = absres; // prepare for next iteration
83 
84  // Print iteration information if needed
85  print_itinfo(prtlvl, STOP_REL_RES, iter, relres1, absres, factor);
86 
87  // Check convergence
88  if ( relres1 < tol ) break;
89  }
90 
91  if ( prtlvl > PRINT_NONE ) {
92  ITS_FINAL(iter, MaxIt, relres1);
93  fasp_gettime(&solve_end);
94  print_cputime("AMG solve",solve_end - solve_start);
95  }
96 
97 #if DEBUG_MODE > 0
98  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
99 #endif
100 
101  return iter;
102 }
103 
126  AMG_param *param)
127 {
128  dCSRmat *ptrA = &mgl[0].A;
129  dvector *b = &mgl[0].b, *x = &mgl[0].x, *r = &mgl[0].w;
130 
131  const INT MaxIt = param->maxit;
132  const SHORT prtlvl = param->print_level;
133  const REAL tol = param->tol;
134  const REAL sumb = fasp_blas_dvec_norm2(b); // L2norm(b)
135 
136  // local variables
137  REAL solve_start, solve_end, solve_time;
138  REAL relres1 = BIGREAL, absres0 = sumb, absres, factor;
139  INT iter = 0;
140 
141 #if DEBUG_MODE > 0
142  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
143  printf("### DEBUG: nrow = %d, ncol = %d, nnz = %d\n",
144  mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
145 #endif
146 
147  fasp_gettime(&solve_start);
148 
149  // Print iteration information if needed
150  print_itinfo(prtlvl, STOP_REL_RES, iter, 1.0, sumb, 0.0);
151 
152  // MG solver here
153  while ( (++iter <= MaxIt) & (sumb > SMALLREAL) ) {
154 
155  // Call one AMLI cycle
156  fasp_solver_amli(mgl, param, 0);
157 
158  // Form residual r = b-A*x
159  fasp_dvec_cp(b, r);
160  fasp_blas_dcsr_aAxpy(-1.0, ptrA, x->val, r->val);
161 
162  // Compute norms of r and convergence factor
163  absres = fasp_blas_dvec_norm2(r); // residual ||r||
164  relres1 = absres/sumb; // relative residual ||r||/||b||
165  factor = absres/absres0; // contraction factor
166  absres0 = absres; // prepare for next iteration
167 
168  // Print iteration information if needed
169  print_itinfo(prtlvl, STOP_REL_RES, iter, relres1, absres, factor);
170 
171  // Check convergence
172  if ( relres1 < tol ) break;
173  }
174 
175  if ( prtlvl > PRINT_NONE ) {
176  ITS_FINAL(iter, MaxIt, relres1);
177  fasp_gettime(&solve_end);
178  solve_time = solve_end - solve_start;
179  print_cputime("AMLI solve", solve_time);
180  }
181 
182 #if DEBUG_MODE > 0
183  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
184 #endif
185 
186  return iter;
187 }
188 
210  AMG_param *param)
211 {
212  dCSRmat *ptrA = &mgl[0].A;
213  dvector *b = &mgl[0].b, *x = &mgl[0].x, *r = &mgl[0].w;
214 
215  const INT MaxIt = param->maxit;
216  const SHORT prtlvl = param->print_level;
217  const REAL tol = param->tol;
218  const REAL sumb = fasp_blas_dvec_norm2(b); // L2norm(b)
219 
220  // local variables
221  REAL solve_start, solve_end;
222  REAL relres1 = BIGREAL, absres0 = BIGREAL, absres, factor;
223  INT iter = 0;
224 
225 #if DEBUG_MODE > 0
226  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
227  printf("### DEBUG: nrow = %d, ncol = %d, nnz = %d\n",
228  mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
229 #endif
230 
231  fasp_gettime(&solve_start);
232 
233  // Print iteration information if needed
234  print_itinfo(prtlvl, STOP_REL_RES, iter, 1.0, sumb, 0.0);
235 
236  while ( (++iter <= MaxIt) & (sumb > SMALLREAL) ) // MG solver here
237  {
238  // one multigrid cycle
239  fasp_solver_nl_amli(mgl, param, 0, mgl[0].num_levels);
240 
241  // r = b-A*x
242  fasp_dvec_cp(b, r);
243  fasp_blas_dcsr_aAxpy(-1.0, ptrA, x->val, r->val);
244 
245  absres = fasp_blas_dvec_norm2(r); // residual ||r||
246  relres1 = absres/sumb; // relative residual ||r||/||b||
247  factor = absres/absres0; // contraction factor
248 
249  // output iteration information if needed
250  print_itinfo(prtlvl, STOP_REL_RES, iter, relres1, absres, factor);
251 
252  if ( relres1 < tol ) break; // early exit condition
253 
254  absres0 = absres;
255  }
256 
257  if ( prtlvl > PRINT_NONE ) {
258  ITS_FINAL(iter, MaxIt, relres1);
259  fasp_gettime(&solve_end);
260  print_cputime("Nonlinear AMLI solve", solve_end - solve_start);
261  }
262 
263 #if DEBUG_MODE > 0
264  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
265 #endif
266 
267  return iter;
268 }
269 
282  AMG_param *param)
283 {
284  dCSRmat *ptrA = &mgl[0].A;
285  dvector *b = &mgl[0].b, *x = &mgl[0].x, *r = &mgl[0].w;
286 
287  const SHORT prtlvl = param->print_level;
288  const REAL sumb = fasp_blas_dvec_norm2(b); // L2norm(b)
289 
290  // local variables
291  REAL solve_start, solve_end;
292  REAL relres1 = BIGREAL, absres;
293 
294 #if DEBUG_MODE > 0
295  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
296  printf("### DEBUG: nrow = %d, ncol = %d, nnz = %d\n",
297  mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
298 #endif
299 
300  fasp_gettime(&solve_start);
301 
302  // Call one full multigrid cycle
303  fasp_solver_fmgcycle(mgl, param);
304 
305  // Form residual r = b-A*x
306  fasp_dvec_cp(b, r);
307  fasp_blas_dcsr_aAxpy(-1.0, ptrA, x->val, r->val);
308 
309  // Compute norms of r and convergence factor
310  absres = fasp_blas_dvec_norm2(r); // residual ||r||
311  relres1 = absres/sumb; // relative residual ||r||/||b||
312 
313  if ( prtlvl > PRINT_NONE ) {
314  printf("FMG finishes with relative residual %e.\n", relres1);
315  fasp_gettime(&solve_end);
316  print_cputime("FMG solve",solve_end - solve_start);
317  }
318 
319 #if DEBUG_MODE > 0
320  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
321 #endif
322 
323  return;
324 }
325 
326 /*---------------------------------*/
327 /*-- End of File --*/
328 /*---------------------------------*/
void fasp_solver_mgrecur(AMG_data *mgl, AMG_param *param, INT level)
Solve Ax=b with recursive multigrid K-cycle.
Definition: mgrecur.c:33
void fasp_solver_fmgcycle(AMG_data *mgl, AMG_param *param)
#include "forts_ns.h"
Definition: fmgcycle.c:36
Parameters for AMG solver.
Definition: fasp.h:583
INT maxit
max number of iterations of AMG
Definition: fasp.h:592
#define REAL
Definition: fasp.h:67
INT fasp_amg_solve_amli(AMG_data *mgl, AMG_param *param)
AMLI – SOLVE phase.
Definition: amg_solve.c:125
dvector b
pointer to the right-hand side at level level_num
Definition: fasp.h:744
void fasp_dvec_cp(dvector *x, dvector *y)
Copy dvector x to dvector y.
Definition: vec.c:345
REAL * val
nonzero entries of A
Definition: fasp.h:166
Vector with n entries of REAL type.
Definition: fasp.h:342
#define INT
Definition: fasp.h:64
void print_itinfo(const INT ptrlvl, const INT stop_type, const INT iter, const REAL relres, const REAL absres, const REAL factor)
Print out iteration information for iterative solvers.
Definition: message.c:36
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
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
void fasp_gettime(REAL *time)
Get system time.
Definition: timing.c:28
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
INT fasp_amg_solve_nl_amli(AMG_data *mgl, AMG_param *param)
Nonlinear AMLI – SOLVE phase.
Definition: amg_solve.c:209
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
#define STOP_REL_RES
Definition of iterative solver stopping criteria types.
Definition: fasp_const.h:131
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
INT fasp_amg_solve(AMG_data *mgl, AMG_param *param)
AMG – SOLVE phase.
Definition: amg_solve.c:36
dCSRmat A
pointer to the matrix at level level_num
Definition: fasp.h:735
void print_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: message.c:165
void fasp_solver_mgcycle(AMG_data *mgl, AMG_param *param)
#include "forts_ns.h"
Definition: mgcycle.c:41
SHORT print_level
print level for AMG
Definition: fasp.h:589
#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_famg_solve(AMG_data *mgl, AMG_param *param)
FMG – SOLVE phase.
Definition: amg_solve.c:281
REAL tol
stopping tolerance for AMG solver
Definition: fasp.h:595
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
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
#define SMALLREAL
Definition: fasp_const.h:238