Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
PreAMGSetupUABSR.c
Go to the documentation of this file.
1
22#include <math.h>
23#include <time.h>
24
25#include "fasp.h"
26#include "fasp_functs.h"
27
28/*---------------------------------*/
29/*-- Declare Private Functions --*/
30/*---------------------------------*/
31
32#include "PreAMGAggregation.inl"
33#include "PreAMGAggregationBSR.inl"
34#include "PreAMGAggregationUA.inl"
35
36static SHORT amg_setup_unsmoothP_unsmoothR_bsr (AMG_data_bsr *, AMG_param *);
37
38/*---------------------------------*/
39/*-- Public Functions --*/
40/*---------------------------------*/
41
56 AMG_param *param)
57{
58#if DEBUG_MODE > 0
59 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
60#endif
61
62 SHORT status = amg_setup_unsmoothP_unsmoothR_bsr(mgl, param);
63
64#if DEBUG_MODE > 0
65 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
66#endif
67
68 return status;
69}
70
71/*---------------------------------*/
72/*-- Private Functions --*/
73/*---------------------------------*/
74
92static SHORT amg_setup_unsmoothP_unsmoothR_bsr (AMG_data_bsr *mgl,
93 AMG_param *param)
94{
95 const SHORT CondType = 1; // Condensation method used for AMG
96
97 const SHORT prtlvl = param->print_level;
98 const SHORT csolver = param->coarse_solver;
99 const SHORT min_cdof = MAX(param->coarse_dof,50);
100 const INT m = mgl[0].A.ROW;
101 const INT nb = mgl[0].A.nb;
102
103 SHORT max_levels = param->max_levels;
104 SHORT i, lvl = 0, status = FASP_SUCCESS;
105 REAL setup_start, setup_end;
106
107 AMG_data *mgl_csr = fasp_amg_data_create(max_levels);
108
109 dCSRmat temp1, temp2;
110
111#if DEBUG_MODE > 0
112 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
113 printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n",
114 mgl[0].A.ROW, mgl[0].A.COL, mgl[0].A.NNZ);
115#endif
116
117 fasp_gettime(&setup_start);
118
119 /*-----------------------*/
120 /*--local working array--*/
121 /*-----------------------*/
122
123 // level info (fine: 0; coarse: 1)
124 ivector *vertices = (ivector *)fasp_mem_calloc(max_levels, sizeof(ivector));
125
126 //each elvel stores the information of the number of aggregations
127 INT *num_aggs = (INT *)fasp_mem_calloc(max_levels, sizeof(INT));
128
129 // each level stores the information of the strongly coupled neighborhoods
130 dCSRmat *Neighbor = (dCSRmat *)fasp_mem_calloc(max_levels, sizeof(dCSRmat));
131
132 for ( i=0; i<max_levels; ++i ) num_aggs[i] = 0;
133
134 /*------------------------------------------*/
135 /*-- setup null spaces for whole Jacobian --*/
136 /*------------------------------------------*/
137
138 /*
139 mgl[0].near_kernel_dim = 1;
140 mgl[0].near_kernel_basis = (REAL **)fasp_mem_calloc(mgl->near_kernel_dim, sizeof(REAL*));
141
142 for ( i=0; i < mgl->near_kernel_dim; ++i ) mgl[0].near_kernel_basis[i] = NULL;
143 */
144
145 /*-----------------------*/
146 /*-- setup ILU param --*/
147 /*-----------------------*/
148
149 // initialize ILU parameters
150 mgl->ILU_levels = param->ILU_levels;
151 ILU_param iluparam;
152
153 if ( param->ILU_levels > 0 ) {
154 iluparam.print_level = param->print_level;
155 iluparam.ILU_lfil = param->ILU_lfil;
156 iluparam.ILU_droptol = param->ILU_droptol;
157 iluparam.ILU_relax = param->ILU_relax;
158 iluparam.ILU_type = param->ILU_type;
159 }
160
161 /*----------------------------*/
162 /*--- checking aggregation ---*/
163 /*----------------------------*/
164 if (param->aggregation_type == PAIRWISE)
165 param->pair_number = MIN(param->pair_number, max_levels);
166
167 // Main AMG setup loop
168 while ( (mgl[lvl].A.ROW > min_cdof) && (lvl < max_levels-1) ) {
169
170 /*-- setup ILU decomposition if necessary */
171 if ( lvl < param->ILU_levels ) {
172 status = fasp_ilu_dbsr_setup(&mgl[lvl].A, &mgl[lvl].LU, &iluparam);
173 if ( status < 0 ) {
174 if ( prtlvl > PRINT_MIN ) {
175 printf("### WARNING: ILU setup on level-%d failed!\n", lvl);
176 printf("### WARNING: Disable ILU for level >= %d.\n", lvl);
177 }
178 param->ILU_levels = lvl;
179 }
180 }
181
182 /*-- get the diagonal inverse --*/
183 mgl[lvl].diaginv = fasp_dbsr_getdiaginv(&mgl[lvl].A);
184
185 switch ( CondType ) {
186 case 2:
187 mgl[lvl].PP = condenseBSR(&mgl[lvl].A); break;
188 default:
189 mgl[lvl].PP = condenseBSRLinf(&mgl[lvl].A); break;
190 }
191
192 /*-- Aggregation --*/
193 switch ( param->aggregation_type ) {
194
195 case VMB: // VMB aggregation
196
197 status = aggregation_vmb (&mgl[lvl].PP, &vertices[lvl], param,
198 lvl+1, &Neighbor[lvl], &num_aggs[lvl]);
199
200 /*-- Choose strength threshold adaptively --*/
201 if ( num_aggs[lvl]*4 > mgl[lvl].PP.row )
202 param->strong_coupled /= 4;
203 else if ( num_aggs[lvl]*1.25 < mgl[lvl].PP.row )
204 param->strong_coupled *= 1.5;
205
206 break;
207
208 case NPAIR: // non-symmetric pairwise matching aggregation
209
210 mgl_csr[lvl].A = mgl[lvl].PP;
211 status = aggregation_nsympair (mgl_csr, param, lvl, vertices,
212 &num_aggs[lvl]);
213
214 break;
215
216 default: // symmetric pairwise matching aggregation
217
218 mgl_csr[lvl].A = mgl[lvl].PP;
219 status = aggregation_symmpair (mgl_csr, param, lvl, vertices,
220 &num_aggs[lvl]);
221
222 // TODO: Need to design better algorithm for pairwise BSR -- Xiaozhe
223 // TODO: Unsymmetric pairwise aggregation not finished -- Chensong
224 // TODO: Check why this fails for BSR --Chensong
225
226 break;
227 }
228
229 if ( status < 0 ) {
230 // When error happens, force solver to use the current multigrid levels!
231 if ( prtlvl > PRINT_MIN ) {
232 printf("### WARNING: Forming aggregates on level-%d failed!\n", lvl);
233 }
234 status = FASP_SUCCESS; break;
235 }
236
237 /* -- Form Prolongation --*/
238 if ( lvl == 0 && mgl[0].near_kernel_dim >0 ) {
239 form_tentative_p_bsr1(&vertices[lvl], &mgl[lvl].P, &mgl[0],
240 num_aggs[lvl], mgl[0].near_kernel_dim,
241 mgl[0].near_kernel_basis);
242 }
243 else {
244 form_boolean_p_bsr(&vertices[lvl], &mgl[lvl].P, &mgl[0], num_aggs[lvl]);
245 }
246
247 /*-- Form resitriction --*/
248 fasp_dbsr_trans(&mgl[lvl].P, &mgl[lvl].R);
249
250 /*-- Form coarse level stiffness matrix --*/
251 fasp_blas_dbsr_rap(&mgl[lvl].R, &mgl[lvl].A, &mgl[lvl].P, &mgl[lvl+1].A);
252
253 /* -- Form extra near kernal space if needed --*/
254 if (mgl[lvl].A_nk != NULL){
255
256 mgl[lvl+1].A_nk = (dCSRmat *)fasp_mem_calloc(1, sizeof(dCSRmat));
257 mgl[lvl+1].P_nk = (dCSRmat *)fasp_mem_calloc(1, sizeof(dCSRmat));
258 mgl[lvl+1].R_nk = (dCSRmat *)fasp_mem_calloc(1, sizeof(dCSRmat));
259
260 temp1 = fasp_format_dbsr_dcsr(&mgl[lvl].R);
261 fasp_blas_dcsr_mxm(&temp1, mgl[lvl].P_nk, mgl[lvl+1].P_nk);
262 fasp_dcsr_trans(mgl[lvl+1].P_nk, mgl[lvl+1].R_nk);
263 temp2 = fasp_format_dbsr_dcsr(&mgl[lvl+1].A);
264 fasp_blas_dcsr_rap(mgl[lvl+1].R_nk, &temp2, mgl[lvl+1].P_nk, mgl[lvl+1].A_nk);
265 fasp_dcsr_free(&temp1);
266 fasp_dcsr_free(&temp2);
267
268 }
269
270 fasp_dcsr_free(&Neighbor[lvl]);
271 fasp_ivec_free(&vertices[lvl]);
272
273 ++lvl;
274 }
275
276 // Setup coarse level systems for direct solvers (BSR version)
277 switch (csolver) {
278
279#if WITH_MUMPS
280 case SOLVER_MUMPS: {
281 // Setup MUMPS direct solver on the coarsest level
282 mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
283 mgl[lvl].mumps.job = 1;
284 fasp_solver_mumps_steps(&mgl[lvl].Ac, &mgl[lvl].b, &mgl[lvl].x, &mgl[lvl].mumps);
285 break;
286 }
287#endif
288
289#if WITH_UMFPACK
290 case SOLVER_UMFPACK: {
291 // Need to sort the matrix A for UMFPACK to work
292 mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
293 dCSRmat Ac_tran;
294 fasp_dcsr_trans(&mgl[lvl].Ac, &Ac_tran);
295 fasp_dcsr_sort(&Ac_tran);
296 fasp_dcsr_cp(&Ac_tran, &mgl[lvl].Ac);
297 fasp_dcsr_free(&Ac_tran);
298 mgl[lvl].Numeric = fasp_umfpack_factorize(&mgl[lvl].Ac, 0);
299 break;
300 }
301#endif
302
303#if WITH_SuperLU
304 case SOLVER_SUPERLU: {
305 /* Setup SuperLU direct solver on the coarsest level */
306 mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
307 }
308#endif
309
310#if WITH_PARDISO
311 case SOLVER_PARDISO: {
312 mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
313 fasp_dcsr_sort(&mgl[lvl].Ac);
314 fasp_pardiso_factorize(&mgl[lvl].Ac, &mgl[lvl].pdata, prtlvl);
315 break;
316 }
317#endif
318
319 default:
320 // Do nothing!
321 break;
322 }
323
324
325 // setup total level number and current level
326 mgl[0].num_levels = max_levels = lvl+1;
327 mgl[0].w = fasp_dvec_create(3*m*nb);
328
329 if (mgl[0].A_nk != NULL){
330
331#if WITH_UMFPACK
332 // Need to sort the matrix A_nk for UMFPACK
333 fasp_dcsr_trans(mgl[0].A_nk, &temp1);
334 fasp_dcsr_sort(&temp1);
335 fasp_dcsr_cp(&temp1, mgl[0].A_nk);
336 fasp_dcsr_free(&temp1);
337#endif
338
339 }
340
341 for ( lvl = 1; lvl < max_levels; lvl++ ) {
342 const INT mm = mgl[lvl].A.ROW*nb;
343 mgl[lvl].num_levels = max_levels;
344 mgl[lvl].b = fasp_dvec_create(mm);
345 mgl[lvl].x = fasp_dvec_create(mm);
346 mgl[lvl].w = fasp_dvec_create(3*mm);
347 mgl[lvl].ILU_levels = param->ILU_levels - lvl; // initialize ILU levels!
348
349 if (mgl[lvl].A_nk != NULL){
350
351#if WITH_UMFPACK
352 // Need to sort the matrix A_nk for UMFPACK
353 fasp_dcsr_trans(mgl[lvl].A_nk, &temp1);
354 fasp_dcsr_sort(&temp1);
355 fasp_dcsr_cp(&temp1, mgl[lvl].A_nk);
356 fasp_dcsr_free(&temp1);
357#endif
358
359 }
360
361 }
362
363 if ( prtlvl > PRINT_NONE ) {
364 fasp_gettime(&setup_end);
365 fasp_amgcomplexity_bsr(mgl,prtlvl);
366 fasp_cputime("Unsmoothed aggregation (BSR) setup", setup_end - setup_start);
367 }
368
369 fasp_mem_free(vertices); vertices = NULL;
370 fasp_mem_free(num_aggs); num_aggs = NULL;
371 fasp_mem_free(Neighbor); Neighbor = NULL;
372
373#if DEBUG_MODE > 0
374 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
375#endif
376
377 return status;
378}
379
380/*---------------------------------*/
381/*-- End of File --*/
382/*---------------------------------*/
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:155
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
void fasp_cputime(const char *message, const REAL cputime)
Print CPU walltime.
Definition: AuxMessage.c:179
void fasp_amgcomplexity_bsr(const AMG_data_bsr *mgl, const SHORT prtlvl)
Print complexities of AMG method for BSR matrices.
Definition: AuxMessage.c:136
void fasp_gettime(REAL *time)
Get system time.
Definition: AuxTiming.c:36
void fasp_ivec_free(ivector *u)
Free vector data space of INT type.
Definition: AuxVector.c:164
dvector fasp_dvec_create(const INT m)
Create dvector data space of REAL type.
Definition: AuxVector.c:62
dCSRmat fasp_format_dbsr_dcsr(const dBSRmat *B)
Transfer a 'dBSRmat' type matrix into a dCSRmat.
Definition: BlaFormat.c:497
SHORT fasp_ilu_dbsr_setup(dBSRmat *A, ILU_data *iludata, ILU_param *iluparam)
Get ILU decoposition of a BSR matrix A.
INT fasp_dbsr_trans(const dBSRmat *A, dBSRmat *AT)
Find A^T from given dBSRmat matrix A.
Definition: BlaSparseBSR.c:199
dvector fasp_dbsr_getdiaginv(const dBSRmat *A)
Get D^{-1} of matrix A.
Definition: BlaSparseBSR.c:486
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:184
void fasp_dcsr_sort(dCSRmat *A)
Sort each row of A in ascending order w.r.t. column indices.
Definition: BlaSparseCSR.c:385
void fasp_dcsr_cp(const dCSRmat *A, dCSRmat *B)
copy a dCSRmat to a new one B=A
Definition: BlaSparseCSR.c:851
INT fasp_dcsr_trans(const dCSRmat *A, dCSRmat *AT)
Find transpose of dCSRmat matrix A.
Definition: BlaSparseCSR.c:952
void fasp_blas_dbsr_rap(const dBSRmat *R, const dBSRmat *A, const dBSRmat *P, dBSRmat *B)
dBSRmat sparse matrix multiplication B=R*A*P
Definition: BlaSpmvBSR.c:4961
void fasp_blas_dcsr_mxm(const dCSRmat *A, const dCSRmat *B, dCSRmat *C)
Sparse matrix multiplication C=A*B.
Definition: BlaSpmvCSR.c:888
void fasp_blas_dcsr_rap(const dCSRmat *R, const dCSRmat *A, const dCSRmat *P, dCSRmat *RAP)
Triple sparse matrix multiplication B=R*A*P.
Definition: BlaSpmvCSR.c:994
SHORT fasp_amg_setup_ua_bsr(AMG_data_bsr *mgl, AMG_param *param)
Set up phase of unsmoothed aggregation AMG (BSR format)
AMG_data * fasp_amg_data_create(SHORT max_levels)
Create and initialize AMG_data for classical and SA AMG.
Definition: PreDataInit.c:64
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:188
Main header file for the FASP project.
#define MIN(a, b)
Definition: fasp.h:75
#define REAL
Definition: fasp.h:67
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:74
#define INT
Definition: fasp.h:64
#define SOLVER_PARDISO
Definition: fasp_const.h:126
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
#define NPAIR
Definition: fasp_const.h:171
#define SOLVER_MUMPS
Definition: fasp_const.h:125
#define PAIRWISE
Definition of aggregation types.
Definition: fasp_const.h:169
#define SOLVER_SUPERLU
Definition: fasp_const.h:123
#define SOLVER_UMFPACK
Definition: fasp_const.h:124
#define VMB
Definition: fasp_const.h:170
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
#define PRINT_MIN
Definition: fasp_const.h:74
Data for multigrid levels in dBSRmat format.
Definition: fasp_block.h:146
dCSRmat PP
pointer to the pressure block (only for reservoir simulation)
Definition: fasp_block.h:182
dBSRmat A
pointer to the matrix at level level_num
Definition: fasp_block.h:155
dCSRmat * A_nk
Matrix data for near kernal.
Definition: fasp_block.h:218
dCSRmat Ac
pointer to the matrix at level level_num (csr format)
Definition: fasp_block.h:173
Mumps_data mumps
data for MUMPS
Definition: fasp_block.h:231
dvector b
pointer to the right-hand side at level level_num
Definition: fasp_block.h:164
void * Numeric
pointer to the numerical dactorization from UMFPACK
Definition: fasp_block.h:176
INT ILU_levels
number of levels use ILU smoother
Definition: fasp_block.h:203
INT num_levels
number of levels in use <= max_levels
Definition: fasp_block.h:152
dCSRmat * R_nk
Resriction for near kernal.
Definition: fasp_block.h:224
dvector x
pointer to the iterative solution at level level_num
Definition: fasp_block.h:167
dCSRmat * P_nk
Prolongation for near kernal.
Definition: fasp_block.h:221
dvector diaginv
pointer to the diagonal inverse at level level_num
Definition: fasp_block.h:170
dvector w
temporary work space
Definition: fasp_block.h:228
Data for AMG methods.
Definition: fasp.h:790
dCSRmat A
pointer to the matrix at level level_num
Definition: fasp.h:803
Parameters for AMG methods.
Definition: fasp.h:447
INT ILU_lfil
level of fill-in for ILUs and ILUk
Definition: fasp.h:555
REAL ILU_relax
relaxation for ILUs
Definition: fasp.h:561
SHORT print_level
print level for AMG
Definition: fasp.h:453
SHORT aggregation_type
aggregation type
Definition: fasp.h:510
SHORT ILU_levels
number of levels use ILU smoother
Definition: fasp.h:549
SHORT coarse_solver
coarse solver type
Definition: fasp.h:492
REAL strong_coupled
strong coupled threshold for aggregate
Definition: fasp.h:534
SHORT ILU_type
ILU type for smoothing.
Definition: fasp.h:552
INT coarse_dof
max number of coarsest level DOF
Definition: fasp.h:465
REAL ILU_droptol
drop tolerance for ILUt
Definition: fasp.h:558
INT pair_number
number of pairwise matchings
Definition: fasp.h:531
SHORT max_levels
max number of levels of AMG
Definition: fasp.h:462
Parameters for ILU.
Definition: fasp.h:396
INT ILU_lfil
level of fill-in for ILUk
Definition: fasp.h:405
REAL ILU_relax
add the sum of dropped elements to diagonal element in proportion relax
Definition: fasp.h:411
SHORT print_level
print level
Definition: fasp.h:399
SHORT ILU_type
ILU type for decomposition.
Definition: fasp.h:402
REAL ILU_droptol
drop tolerance for ILUt
Definition: fasp.h:408
INT job
work for MUMPS
Definition: fasp.h:601
INT nb
dimension of each sub-block
Definition: fasp_block.h:46
INT ROW
number of rows of sub-blocks in matrix A, M
Definition: fasp_block.h:37
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:143
Vector with n entries of INT type.
Definition: fasp.h:361