Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
PreMGCycleFull.c
Go to the documentation of this file.
1
16#include <math.h>
17#include <time.h>
18
19#include "fasp.h"
20#include "fasp_functs.h"
21
22/*---------------------------------*/
23/*-- Declare Private Functions --*/
24/*---------------------------------*/
25
26#include "PreMGUtil.inl"
27#include "PreMGSmoother.inl"
28
29/*---------------------------------*/
30/*-- Public Functions --*/
31/*---------------------------------*/
32
47void fasp_solver_fmgcycle (AMG_data *mgl,
48 AMG_param *param)
49{
50 const SHORT maxit = 3; // Max allowed V-cycles in each level
51 const SHORT amg_type = param->AMG_type;
52 const SHORT prtlvl = param->print_level;
53 const SHORT nl = mgl[0].num_levels;
54 const SHORT smoother = param->smoother;
55 const SHORT smooth_order = param->smooth_order;
56 const SHORT coarse_solver = param->coarse_solver;
57
58 const REAL relax = param->relaxation;
59 const SHORT ndeg = param->polynomial_degree;
60 const REAL tol = param->tol*1e-4;
61
62 // local variables
63 INT l, i, lvl, num_cycle;
64 REAL alpha = 1.0, relerr;
65
66 // Schwarz parameters
67 SWZ_param swzparam;
68 if ( param->SWZ_levels > 0 ) {
69 swzparam.SWZ_blksolver = param->SWZ_blksolver;
70 }
71
72#if DEBUG_MODE > 0
73 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
74 printf("### DEBUG: n=%d, nnz=%d\n", mgl[0].A.row, mgl[0].A.nnz);
75#endif
76
77 if ( prtlvl >= PRINT_MOST )
78 printf("FMG_level = %d, ILU_level = %d\n", nl, param->ILU_levels);
79
80 // restriction r1 = R*r0
81 switch (amg_type) {
82
83 case UA_AMG:
84 for (l=0;l<nl-1;l++)
85 fasp_blas_dcsr_mxv_agg(&mgl[l].R, mgl[l].b.val, mgl[l+1].b.val);
86 break;
87
88 default:
89 for (l=0;l<nl-1;l++)
90 fasp_blas_dcsr_mxv(&mgl[l].R, mgl[l].b.val, mgl[l+1].b.val);
91 break;
92
93 }
94
95 fasp_dvec_set(mgl[l].A.row, &mgl[l].x, 0.0); // initial guess
96
97 // If only one level, just direct solver
98 if ( nl==1 ) {
99
100 switch (coarse_solver) {
101
102#if WITH_PARDISO
103 case SOLVER_PARDISO: {
104 /* use Intel MKL PARDISO direct solver on the coarsest level */
105 fasp_pardiso_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].pdata, 0);
106 break;
107 }
108#endif
109
110#if WITH_SuperLU
111 case SOLVER_SUPERLU:
112 /* use SuperLU direct solver on the coarsest level */
113 fasp_solver_superlu(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, 0);
114 break;
115#endif
116
117#if WITH_UMFPACK
118 case SOLVER_UMFPACK:
119 /* use UMFPACK direct solver on the coarsest level */
120 fasp_umfpack_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, mgl[nl-1].Numeric, 0);
121 break;
122#endif
123
124#if WITH_MUMPS
125 case SOLVER_MUMPS:
126 /* use MUMPS direct solver on the coarsest level */
127 mgl[nl-1].mumps.job = 2;
128 fasp_solver_mumps_steps(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].mumps);
129 break;
130#endif
131
132 default:
133 /* use iterative solver on the coarest level */
134 fasp_coarse_itsolver(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, tol, prtlvl);
135
136 }
137
138 return;
139
140 }
141
142 for ( i=1; i<nl; i++ ) {
143
144 // Coarse Space Solver:
145 switch (coarse_solver) {
146
147#if WITH_PARDISO
148 case SOLVER_PARDISO: {
149 /* use Intel MKL PARDISO direct solver on the coarsest level */
150 fasp_pardiso_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].pdata, 0);
151 break;
152 }
153#endif
154
155#if WITH_SuperLU
156 case SOLVER_SUPERLU:
157 /* use SuperLU direct solver on the coarsest level */
158 fasp_solver_superlu(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, 0);
159 break;
160#endif
161
162#if WITH_UMFPACK
163 case SOLVER_UMFPACK:
164 /* use UMFPACK direct solver on the coarsest level */
165 fasp_umfpack_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, mgl[nl-1].Numeric, 0);
166 break;
167#endif
168
169#if WITH_MUMPS
170 case SOLVER_MUMPS:
171 /* use MUMPS direct solver on the coarsest level */
172 mgl[nl-1].mumps.job = 2;
173 fasp_solver_mumps_steps(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].mumps);
174 break;
175#endif
176
177 default:
178 /* use iterative solver on the coarest level */
179 fasp_coarse_itsolver(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, tol, prtlvl);
180
181 }
182
183 // Slash part: /-cycle
184 {
185 --l; // go back to finer level
186
187 // find the optimal scaling factor alpha
188 if ( param->coarse_scaling == ON ) {
189 alpha = fasp_blas_darray_dotprod(mgl[l+1].A.row, mgl[l+1].x.val, mgl[l+1].b.val)
190 / fasp_blas_dcsr_vmv(&mgl[l+1].A, mgl[l+1].x.val, mgl[l+1].x.val);
191 alpha = MIN(alpha, 1.0); // Add this for safty! --Chensong on 10/04/2014
192 }
193
194 // prolongation u = u + alpha*P*e1
195 switch (amg_type) {
196 case UA_AMG:
197 fasp_blas_dcsr_aAxpy_agg(alpha, &mgl[l].P, mgl[l+1].x.val, mgl[l].x.val); break;
198 default:
199 fasp_blas_dcsr_aAxpy(alpha, &mgl[l].P, mgl[l+1].x.val, mgl[l].x.val); break;
200 }
201 }
202
203 // initialzie rel error
204 num_cycle = 0; relerr = BIGREAL;
205
206 while ( relerr > param->tol && num_cycle < maxit) {
207
208 ++num_cycle;
209
210 // form residual r = b - A x
211 fasp_darray_cp(mgl[l].A.row, mgl[l].b.val, mgl[l].w.val);
212 fasp_blas_dcsr_aAxpy(-1.0,&mgl[l].A, mgl[l].x.val, mgl[l].w.val);
213 relerr = fasp_blas_dvec_norm2(&mgl[l].w) / fasp_blas_dvec_norm2(&mgl[l].b);
214
215 // Forward Sweep
216 for ( lvl=0; lvl<i; lvl++ ) {
217
218 // pre smoothing
219 if (l<param->ILU_levels) {
220 fasp_smoother_dcsr_ilu(&mgl[l].A, &mgl[l].b, &mgl[l].x, &mgl[l].LU);
221 }
222 else if (l<mgl->SWZ_levels) {
223 switch (mgl[l].Schwarz.SWZ_type) {
224 case SCHWARZ_SYMMETRIC:
225 fasp_dcsr_swz_forward(&mgl[l].Schwarz, &swzparam, &mgl[l].x, &mgl[l].b);
226 fasp_dcsr_swz_backward(&mgl[l].Schwarz, &swzparam,&mgl[l].x, &mgl[l].b);
227 break;
228 default:
229 fasp_dcsr_swz_forward(&mgl[l].Schwarz, &swzparam, &mgl[l].x, &mgl[l].b);
230 break;
231 }
232 }
233
234 else {
235 fasp_dcsr_presmoothing(smoother,&mgl[l].A,&mgl[l].b,&mgl[l].x,param->presmooth_iter,
236 0,mgl[l].A.row-1,1,relax,ndeg,smooth_order,mgl[l].cfmark.val);
237 }
238
239 // form residual r = b - A x
240 fasp_darray_cp(mgl[l].A.row, mgl[l].b.val, mgl[l].w.val);
241 fasp_blas_dcsr_aAxpy(-1.0,&mgl[l].A, mgl[l].x.val, mgl[l].w.val);
242
243 // restriction r1 = R*r0
244 switch (amg_type) {
245 case UA_AMG:
246 fasp_blas_dcsr_mxv_agg(&mgl[l].R, mgl[l].w.val, mgl[l+1].b.val);
247 break;
248 default:
249 fasp_blas_dcsr_mxv(&mgl[l].R, mgl[l].w.val, mgl[l+1].b.val);
250 break;
251 }
252
253 ++l;
254
255 // prepare for the next level
256 fasp_dvec_set(mgl[l].A.row, &mgl[l].x, 0.0);
257
258 } // end for lvl
259
260 // CoarseSpaceSolver:
261 switch (coarse_solver) {
262
263#if WITH_PARDISO
264 case SOLVER_PARDISO: {
265 /* use Intel MKL PARDISO direct solver on the coarsest level */
266 fasp_pardiso_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].pdata, 0);
267 break;
268 }
269#endif
270
271#if WITH_SuperLU
272 case SOLVER_SUPERLU:
273 /* use SuperLU direct solver on the coarsest level */
274 fasp_solver_superlu(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, 0);
275 break;
276#endif
277
278#if WITH_UMFPACK
279 case SOLVER_UMFPACK:
280 /* use UMFPACK direct solver on the coarsest level */
281 fasp_umfpack_solve(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, mgl[nl-1].Numeric, 0);
282 break;
283#endif
284
285#if WITH_MUMPS
286 case SOLVER_MUMPS:
287 /* use MUMPS direct solver on the coarsest level */
288 mgl[nl-1].mumps.job = 2;
289 fasp_solver_mumps_steps(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, &mgl[nl-1].mumps);
290 break;
291#endif
292
293 default:
294 /* use iterative solver on the coarest level */
295 fasp_coarse_itsolver(&mgl[nl-1].A, &mgl[nl-1].b, &mgl[nl-1].x, tol, prtlvl);
296
297 }
298
299 // Backward Sweep
300 for ( lvl=0; lvl<i; lvl++ ) {
301
302 --l;
303
304 // find the optimal scaling factor alpha
305 if ( param->coarse_scaling == ON ) {
306 alpha = fasp_blas_darray_dotprod(mgl[l+1].A.row, mgl[l+1].x.val, mgl[l+1].b.val)
307 / fasp_blas_dcsr_vmv(&mgl[l+1].A, mgl[l+1].x.val, mgl[l+1].x.val);
308 alpha = MIN(alpha, 1.0); // Add this for safty! --Chensong on 10/04/2014
309 }
310
311 // prolongation u = u + alpha*P*e1
312 switch (amg_type)
313 {
314 case UA_AMG:
315 fasp_blas_dcsr_aAxpy_agg(alpha, &mgl[l].P, mgl[l+1].x.val, mgl[l].x.val);
316 break;
317 default:
318 fasp_blas_dcsr_aAxpy(alpha, &mgl[l].P, mgl[l+1].x.val, mgl[l].x.val);
319 break;
320 }
321
322 // post-smoothing
323 if (l<param->ILU_levels) {
324 fasp_smoother_dcsr_ilu(&mgl[l].A, &mgl[l].b, &mgl[l].x, &mgl[l].LU);
325 }
326 else if (l<mgl->SWZ_levels) {
327 switch (mgl[l].Schwarz.SWZ_type) {
328 case SCHWARZ_SYMMETRIC:
329 fasp_dcsr_swz_backward(&mgl[l].Schwarz, &swzparam,&mgl[l].x, &mgl[l].b);
330 fasp_dcsr_swz_forward(&mgl[l].Schwarz, &swzparam, &mgl[l].x, &mgl[l].b);
331 break;
332 default:
333 fasp_dcsr_swz_backward(&mgl[l].Schwarz, &swzparam,&mgl[l].x, &mgl[l].b);
334 break;
335 }
336 }
337
338 else {
339 fasp_dcsr_postsmoothing(smoother,&mgl[l].A,&mgl[l].b,&mgl[l].x,param->postsmooth_iter,
340 0,mgl[l].A.row-1,-1,relax,ndeg,smooth_order,mgl[l].cfmark.val);
341 }
342
343 } // end while
344
345 } //end while
346
347 } // end for
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_darray_cp
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:210
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_blas_darray_dotprod
REAL fasp_blas_darray_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies x and y.
Definition: BlaArray.c:771
fasp_dcsr_swz_forward
void fasp_dcsr_swz_forward(SWZ_data *swzdata, SWZ_param *swzparam, dvector *x, dvector *b)
Schwarz smoother: forward sweep.
Definition: BlaSchwarzSetup.c:218
fasp_dcsr_swz_backward
void fasp_dcsr_swz_backward(SWZ_data *swzdata, SWZ_param *swzparam, dvector *x, dvector *b)
Schwarz smoother: backward sweep.
Definition: BlaSchwarzSetup.c:328
fasp_blas_dcsr_mxv_agg
void fasp_blas_dcsr_mxv_agg(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x (nonzeros of A = 1)
Definition: BlaSpmvCSR.c:438
fasp_blas_dcsr_mxv
void fasp_blas_dcsr_mxv(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvCSR.c:242
fasp_blas_dcsr_aAxpy_agg
void fasp_blas_dcsr_aAxpy_agg(const REAL alpha, const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y (nonzeros of A = 1)
Definition: BlaSpmvCSR.c:727
fasp_blas_dcsr_vmv
REAL fasp_blas_dcsr_vmv(const dCSRmat *A, const REAL *x, const REAL *y)
vector-Matrix-vector multiplication alpha = y'*A*x
Definition: BlaSpmvCSR.c:839
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_smoother_dcsr_ilu
void fasp_smoother_dcsr_ilu(dCSRmat *A, dvector *b, dvector *x, void *data)
ILU method as a smoother.
Definition: ItrSmootherCSR.c:1279
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_mumps_steps
int fasp_solver_mumps_steps(dCSRmat *ptrA, dvector *b, dvector *u, Mumps_data *mumps)
Solve Ax=b by MUMPS in three steps.
Definition: XtrMumps.c:196
fasp_solver_superlu
int fasp_solver_superlu(dCSRmat *ptrA, dvector *b, dvector *u, const SHORT prtlvl)
Solve Au=b by SuperLU.
Definition: XtrSuperlu.c:47
fasp.h
Main header file for the FASP project.
MIN
#define MIN(a, b)
Definition: fasp.h:83
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
INT
#define INT
Definition: fasp.h:72
SOLVER_PARDISO
#define SOLVER_PARDISO
Definition: fasp_const.h:126
PRINT_MOST
#define PRINT_MOST
Definition: fasp_const.h:77
SOLVER_MUMPS
#define SOLVER_MUMPS
Definition: fasp_const.h:125
SOLVER_SUPERLU
#define SOLVER_SUPERLU
Definition: fasp_const.h:123
SCHWARZ_SYMMETRIC
#define SCHWARZ_SYMMETRIC
Definition: fasp_const.h:158
BIGREAL
#define BIGREAL
Some global constants.
Definition: fasp_const.h:255
SOLVER_UMFPACK
#define SOLVER_UMFPACK
Definition: fasp_const.h:124
ON
#define ON
Definition of switch.
Definition: fasp_const.h:67
UA_AMG
#define UA_AMG
Definition: fasp_const.h:165
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::mumps
Mumps_data mumps
data for MUMPS
Definition: fasp.h:866
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::num_levels
SHORT num_levels
number of levels in use <= max_levels
Definition: fasp.h:812
AMG_data::w
dvector w
temporary work space
Definition: fasp.h:863
AMG_data::cfmark
ivector cfmark
pointer to the CF marker at level level_num
Definition: fasp.h:840
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::polynomial_degree
SHORT polynomial_degree
degree of the polynomial smoother
Definition: fasp.h:497
AMG_param::coarse_scaling
SHORT coarse_scaling
switch of scaling of the coarse grid correction
Definition: fasp.h:503
AMG_param::ILU_levels
SHORT ILU_levels
number of levels use ILU smoother
Definition: fasp.h:560
AMG_param::AMG_type
SHORT AMG_type
type of AMG method
Definition: fasp.h:458
AMG_param::tol
REAL tol
stopping tolerance for AMG solver
Definition: fasp.h:467
AMG_param::coarse_solver
SHORT coarse_solver
coarse solver type
Definition: fasp.h:500
AMG_param::relaxation
REAL relaxation
relaxation parameter for Jacobi and SOR smoother
Definition: fasp.h:494
AMG_param::smoother
SHORT smoother
smoother type
Definition: fasp.h:482
AMG_param::SWZ_blksolver
INT SWZ_blksolver
type of Schwarz block solver
Definition: fasp.h:590
AMG_param::postsmooth_iter
SHORT postsmooth_iter
number of postsmoothers
Definition: fasp.h:491
AMG_param::SWZ_levels
INT SWZ_levels
number of levels use Schwarz smoother
Definition: fasp.h:578
AMG_param::presmooth_iter
SHORT presmooth_iter
number of presmoothers
Definition: fasp.h:488
AMG_param::smooth_order
SHORT smooth_order
smoother order
Definition: fasp.h:485
Mumps_data::job
INT job
work for MUMPS
Definition: fasp.h:615
SWZ_param
Parameters for Schwarz method.
Definition: fasp.h:430
SWZ_param::SWZ_blksolver
INT SWZ_blksolver
type of Schwarz block solver
Definition: fasp.h:445
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360
ivector::val
INT * val
actual vector entries
Definition: fasp.h:374