Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
PreMGSolve.c
Go to the documentation of this file.
1
19#include <time.h>
20
21#include "fasp.h"
22#include "fasp_functs.h"
23
24/*---------------------------------*/
25/*-- Declare Private Functions --*/
26/*---------------------------------*/
27
28#include "KryUtil.inl"
29
30/*---------------------------------*/
31/*-- Public Functions --*/
32/*---------------------------------*/
33
49INT fasp_amg_solve (AMG_data *mgl,
50 AMG_param *param)
51{
52 dCSRmat *ptrA = &mgl[0].A;
53 dvector *b = &mgl[0].b, *x = &mgl[0].x, *r = &mgl[0].w;
54
55 const SHORT prtlvl = param->print_level;
56 const INT MaxIt = param->maxit;
57 const REAL tol = param->tol;
58 const REAL sumb = fasp_blas_dvec_norm2(b); // L2norm(b)
59
60 // local variables
61 REAL solve_start, solve_end;
62 REAL relres1 = 1.0, absres0 = sumb, absres, factor;
63 INT iter = 0;
64
65#if DEBUG_MODE > 0
66 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
67 printf("### DEBUG: nrow = %d, ncol = %d, nnz = %d\n",
68 mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
69#endif
70
71 fasp_gettime(&solve_start);
72
73 // Print iteration information if needed
74 fasp_itinfo(prtlvl, STOP_REL_RES, iter, relres1, sumb, 0.0);
75
76 // If b = 0, set x = 0 to be a trivial solution
77 if ( sumb <= SMALLREAL ) fasp_dvec_set(x->row, x, 0.0);
78
79 // MG solver here
80 while ( (iter++ < MaxIt) & (sumb > SMALLREAL) ) {
81
82#if TRUE
83 // Call one multigrid cycle -- non recursive version
84 fasp_solver_mgcycle(mgl, param);
85#else
86 // Call one multigrid cycle -- recursive version
87 fasp_solver_mgrecur(mgl, param, 0);
88#endif
89
90 // Form residual r = b - A*x
91 fasp_dvec_cp(b, r);
92 fasp_blas_dcsr_aAxpy(-1.0, ptrA, x->val, r->val);
93
94 // Compute norms of r and convergence factor
95 absres = fasp_blas_dvec_norm2(r); // residual ||r||
96 relres1 = absres/MAX(SMALLREAL,sumb); // relative residual ||r||/||b||
97 factor = absres/absres0; // contraction factor
98 absres0 = absres; // prepare for next iteration
99
100 // Print iteration information if needed
101 fasp_itinfo(prtlvl, STOP_REL_RES, iter, relres1, absres, factor);
102
103 // Check convergence
104 if ( relres1 < tol ) break;
105 }
106
107 if ( prtlvl > PRINT_NONE ) {
108 ITS_FINAL(iter, MaxIt, relres1);
109 fasp_gettime(&solve_end);
110 fasp_cputime("AMG solve", solve_end - solve_start);
111 }
112
113#if DEBUG_MODE > 0
114 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
115#endif
116
117 if ( iter > MaxIt )
118 return ERROR_SOLVER_MAXIT;
119 else
120 return iter;}
121
142INT fasp_amg_solve_amli (AMG_data *mgl,
143 AMG_param *param)
144{
145 dCSRmat *ptrA = &mgl[0].A;
146 dvector *b = &mgl[0].b, *x = &mgl[0].x, *r = &mgl[0].w;
147
148 const INT MaxIt = param->maxit;
149 const SHORT prtlvl = param->print_level;
150 const REAL tol = param->tol;
151 const REAL sumb = fasp_blas_dvec_norm2(b); // L2norm(b)
152
153 // local variables
154 REAL solve_start, solve_end;
155 REAL relres1 = 1.0, absres0 = sumb, absres, factor;
156 INT iter = 0;
157
158#if DEBUG_MODE > 0
159 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
160 printf("### DEBUG: nrow = %d, ncol = %d, nnz = %d\n",
161 mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
162#endif
163
164 fasp_gettime(&solve_start);
165
166 // Print iteration information if needed
167 fasp_itinfo(prtlvl, STOP_REL_RES, iter, relres1, sumb, 0.0);
168
169 // If b = 0, set x = 0 to be a trivial solution
170 if ( sumb <= SMALLREAL ) fasp_dvec_set(x->row, x, 0.0);
171
172 // MG solver here
173 while ( (iter++ < MaxIt) & (sumb > SMALLREAL) ) {
174
175 // Call one AMLI cycle
176 fasp_solver_amli(mgl, param, 0);
177
178 // Form residual r = b-A*x
179 fasp_dvec_cp(b, r);
180 fasp_blas_dcsr_aAxpy(-1.0, ptrA, x->val, r->val);
181
182 // Compute norms of r and convergence factor
183 absres = fasp_blas_dvec_norm2(r); // residual ||r||
184 relres1 = absres/MAX(SMALLREAL,sumb); // relative residual ||r||/||b||
185 factor = absres/absres0; // contraction factor
186 absres0 = absres; // prepare for next iteration
187
188 // Print iteration information if needed
189 fasp_itinfo(prtlvl, STOP_REL_RES, iter, relres1, absres, factor);
190
191 // Check convergence
192 if ( relres1 < tol ) break;
193 }
194
195 if ( prtlvl > PRINT_NONE ) {
196 ITS_FINAL(iter, MaxIt, relres1);
197 fasp_gettime(&solve_end);
198 fasp_cputime("AMLI solve", solve_end - solve_start);
199 }
200
201#if DEBUG_MODE > 0
202 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
203#endif
204
205 if ( iter > MaxIt )
206 return ERROR_SOLVER_MAXIT;
207 else
208 return iter;
209}
210
230INT fasp_amg_solve_namli (AMG_data *mgl,
231 AMG_param *param)
232{
233 dCSRmat *ptrA = &mgl[0].A;
234 dvector *b = &mgl[0].b, *x = &mgl[0].x, *r = &mgl[0].w;
235
236 const INT MaxIt = param->maxit;
237 const SHORT prtlvl = param->print_level;
238 const REAL tol = param->tol;
239 const REAL sumb = fasp_blas_dvec_norm2(b); // L2norm(b)
240
241 // local variables
242 REAL solve_start, solve_end;
243 REAL relres1 = 1.0, absres0 = sumb, absres, factor;
244 INT iter = 0;
245
246#if DEBUG_MODE > 0
247 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
248 printf("### DEBUG: nrow = %d, ncol = %d, nnz = %d\n",
249 mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
250#endif
251
252 fasp_gettime(&solve_start);
253
254 // Print iteration information if needed
255 fasp_itinfo(prtlvl, STOP_REL_RES, iter, relres1, sumb, 0.0);
256
257 // If b = 0, set x = 0 to be a trivial solution
258 if ( sumb <= SMALLREAL ) fasp_dvec_set(x->row, x, 0.0);
259
260 while ( (iter++ < MaxIt) & (sumb > SMALLREAL) ) // MG solver here
261 {
262 // one multigrid cycle
263 fasp_solver_namli(mgl, param, 0, mgl[0].num_levels);
264
265 // r = b-A*x
266 fasp_dvec_cp(b, r);
267 fasp_blas_dcsr_aAxpy(-1.0, ptrA, x->val, r->val);
268
269 absres = fasp_blas_dvec_norm2(r); // residual ||r||
270 relres1 = absres/MAX(SMALLREAL,sumb); // relative residual ||r||/||b||
271 factor = absres/absres0; // contraction factor
272
273 // output iteration information if needed
274 fasp_itinfo(prtlvl, STOP_REL_RES, iter, relres1, absres, factor);
275
276 if ( relres1 < tol ) break; // early exit condition
277
278 absres0 = absres;
279 }
280
281 if ( prtlvl > PRINT_NONE ) {
282 ITS_FINAL(iter, MaxIt, relres1);
283 fasp_gettime(&solve_end);
284 fasp_cputime("Nonlinear AMLI solve", solve_end - solve_start);
285 }
286
287#if DEBUG_MODE > 0
288 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
289#endif
290
291 if ( iter > MaxIt )
292 return ERROR_SOLVER_MAXIT;
293 else
294 return iter;
295}
296
308void fasp_famg_solve (AMG_data *mgl,
309 AMG_param *param)
310{
311 dCSRmat *ptrA = &mgl[0].A;
312 dvector *b = &mgl[0].b, *x = &mgl[0].x, *r = &mgl[0].w;
313
314 const SHORT prtlvl = param->print_level;
315 const REAL sumb = fasp_blas_dvec_norm2(b); // L2norm(b)
316
317 // local variables
318 REAL solve_start, solve_end;
319 REAL relres1 = 1.0, absres;
320
321#if DEBUG_MODE > 0
322 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
323 printf("### DEBUG: nrow = %d, ncol = %d, nnz = %d\n",
324 mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
325#endif
326
327 fasp_gettime(&solve_start);
328
329 // If b = 0, set x = 0 to be a trivial solution
330 if ( sumb <= SMALLREAL ) fasp_dvec_set(x->row, x, 0.0);
331
332 // Call one full multigrid cycle
333 fasp_solver_fmgcycle(mgl, param);
334
335 // Form residual r = b-A*x
336 fasp_dvec_cp(b, r);
337 fasp_blas_dcsr_aAxpy(-1.0, ptrA, x->val, r->val);
338
339 // Compute norms of r and convergence factor
340 absres = fasp_blas_dvec_norm2(r); // residual ||r||
341 relres1 = absres/MAX(SMALLREAL,sumb); // relative residual ||r||/||b||
342
343 if ( prtlvl > PRINT_NONE ) {
344 printf("FMG finishes with relative residual %e.\n", relres1);
345 fasp_gettime(&solve_end);
346 fasp_cputime("FMG solve", solve_end - solve_start);
347 }
348
349#if DEBUG_MODE > 0
350 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
351#endif
352
353 return;
354}
355
356/*---------------------------------*/
357/*-- End of File --*/
358/*---------------------------------*/
fasp_cputime
void fasp_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: AuxMessage.c:179
fasp_itinfo
void fasp_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: AuxMessage.c:41
fasp_gettime
void fasp_gettime(REAL *time)
Get system time.
Definition: AuxTiming.c:37
fasp_dvec_set
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
fasp_dvec_cp
void fasp_dvec_cp(const dvector *x, dvector *y)
Copy dvector x to dvector y.
Definition: AuxVector.c:334
fasp_blas_dcsr_aAxpy
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:494
fasp_blas_dvec_norm2
REAL fasp_blas_dvec_norm2(const dvector *x)
L2 norm of dvector x.
Definition: BlaVector.c:170
fasp_solver_fmgcycle
void fasp_solver_fmgcycle(AMG_data *mgl, AMG_param *param)
Solve Ax=b with non-recursive full multigrid K-cycle.
Definition: PreMGCycleFull.c:47
fasp_solver_mgcycle
void fasp_solver_mgcycle(AMG_data *mgl, AMG_param *param)
Solve Ax=b with non-recursive multigrid cycle.
Definition: PreMGCycle.c:48
fasp_solver_namli
void fasp_solver_namli(AMG_data *mgl, AMG_param *param, INT l, INT num_levels)
Solve Ax=b with recursive nonlinear AMLI-cycle.
Definition: PreMGRecurAMLI.c:291
fasp_solver_amli
void fasp_solver_amli(AMG_data *mgl, AMG_param *param, INT l)
Solve Ax=b with recursive AMLI-cycle.
Definition: PreMGRecurAMLI.c:58
fasp_solver_mgrecur
void fasp_solver_mgrecur(AMG_data *mgl, AMG_param *param, INT level)
Solve Ax=b with recursive multigrid K-cycle.
Definition: PreMGRecur.c:47
fasp_amg_solve_amli
INT fasp_amg_solve_amli(AMG_data *mgl, AMG_param *param)
AMLI – SOLVE phase.
Definition: PreMGSolve.c:142
fasp_amg_solve_namli
INT fasp_amg_solve_namli(AMG_data *mgl, AMG_param *param)
Nonlinear AMLI – SOLVE phase.
Definition: PreMGSolve.c:230
fasp_amg_solve
INT fasp_amg_solve(AMG_data *mgl, AMG_param *param)
AMG – SOLVE phase.
Definition: PreMGSolve.c:49
fasp_famg_solve
void fasp_famg_solve(AMG_data *mgl, AMG_param *param)
FMG – SOLVE phase.
Definition: PreMGSolve.c:308
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
ERROR_SOLVER_MAXIT
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:48
STOP_REL_RES
#define STOP_REL_RES
Definition of iterative solver stopping criteria types.
Definition: fasp_const.h:132
PRINT_NONE
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
SMALLREAL
#define SMALLREAL
Definition: fasp_const.h:256
AMG_data
Data for AMG methods.
Definition: fasp.h:804
AMG_data::A
dCSRmat A
pointer to the matrix at level level_num
Definition: fasp.h:817
AMG_data::b
dvector b
pointer to the right-hand side at level level_num
Definition: fasp.h:826
AMG_data::x
dvector x
pointer to the iterative solution at level level_num
Definition: fasp.h:829
AMG_data::w
dvector w
temporary work space
Definition: fasp.h:863
AMG_param
Parameters for AMG methods.
Definition: fasp.h:455
AMG_param::print_level
SHORT print_level
print level for AMG
Definition: fasp.h:461
AMG_param::tol
REAL tol
stopping tolerance for AMG solver
Definition: fasp.h:467
AMG_param::maxit
INT maxit
max number of iterations of AMG
Definition: fasp.h:464
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:157
dCSRmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354