Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
amg_setup_rs.c
Go to the documentation of this file.
1 
10 #include <time.h>
11 
12 #ifdef _OPENMP
13 #include <omp.h>
14 #endif
15 
16 #include "fasp.h"
17 #include "fasp_functs.h"
18 
19 /*---------------------------------*/
20 /*-- Public Functions --*/
21 /*---------------------------------*/
22 
48  AMG_param *param)
49 {
50  const SHORT prtlvl = param->print_level;
51  const SHORT cycle_type = param->cycle_type;
52  const SHORT csolver = param->coarse_solver;
53  const SHORT min_cdof = MAX(param->coarse_dof,MIN_CDOF);
54  const INT m = mgl[0].A.row;
55 
56  // local variables
57  SHORT status = FASP_SUCCESS;
58  INT lvl = 0, max_lvls = param->max_levels;
59  REAL setup_start, setup_end;
60  ILU_param iluparam;
61  Schwarz_param swzparam;
62  iCSRmat Scouple; // strong n-couplings
63 
64  // level info (fine: 0; coarse: 1)
65  ivector vertices = fasp_ivec_create(m);
66 
67 #if DEBUG_MODE > 0
68  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
69  printf("### DEBUG: n = %d, nnz = %d\n", mgl[0].A.row, mgl[0].A.nnz);
70 #endif
71 
72  fasp_gettime(&setup_start);
73 
74  // Make sure classical AMG will not call fasp_blas_dcsr_mxv_agg!
75  param->tentative_smooth = 1.0;
76 
77  // If user want to use aggressive coarsening but did not specify number of
78  // levels use aggressive coarsening, make sure apply aggressive coarsening
79  // on the finest level only !!!
80  if ( param->coarsening_type == COARSE_AC ) {
81  param->aggressive_level = MAX(param->aggressive_level, 1);
82  }
83 
84  // Initialize AMLI coefficients
85  if ( cycle_type == AMLI_CYCLE ) {
86  const INT amlideg = param->amli_degree;
87  param->amli_coef = (REAL *)fasp_mem_calloc(amlideg+1,sizeof(REAL));
88  fasp_amg_amli_coef(2.0, 0.5, amlideg, param->amli_coef);
89  }
90 
91  // Initialize ILU parameters
92  mgl->ILU_levels = param->ILU_levels;
93  if ( param->ILU_levels > 0 ) {
94  iluparam.print_level = param->print_level;
95  iluparam.ILU_lfil = param->ILU_lfil;
96  iluparam.ILU_droptol = param->ILU_droptol;
97  iluparam.ILU_relax = param->ILU_relax;
98  iluparam.ILU_type = param->ILU_type;
99  }
100 
101  // Initialize Schwarz parameters
102  mgl->Schwarz_levels = param->Schwarz_levels;
103  if ( param->Schwarz_levels > 0 ) {
104  swzparam.Schwarz_mmsize = param->Schwarz_mmsize;
105  swzparam.Schwarz_maxlvl = param->Schwarz_maxlvl;
106  swzparam.Schwarz_type = param->Schwarz_type;
107  swzparam.Schwarz_blksolver = param->Schwarz_blksolver;
108  }
109 
110 #if DIAGONAL_PREF
111  // Reorder each row to keep the diagonal entries appear first !!!
112  fasp_dcsr_diagpref(&mgl[0].A);
113 #endif
114 
115  // Main AMG setup loop
116  while ( (mgl[lvl].A.row > min_cdof) && (lvl < max_lvls-1) ) {
117 
118 #if DEBUG_MODE > 2
119  printf("### DEBUG: level = %d, row = %d, nnz = %d\n",
120  lvl, mgl[lvl].A.row, mgl[lvl].A.nnz);
121 #endif
122 
123  /*-- Setup ILU decomposition if needed --*/
124  if ( lvl < param->ILU_levels ) {
125  status = fasp_ilu_dcsr_setup(&mgl[lvl].A, &mgl[lvl].LU, &iluparam);
126  if ( status < 0 ) {
127  if ( prtlvl > PRINT_MIN ) {
128  printf("### WARNING: ILU setup on level-%d failed!\n", lvl);
129  printf("### WARNING: Disable ILU for level >= %d.\n", lvl);
130  }
131  param->ILU_levels = lvl;
132  }
133  }
134 
135  /*-- Setup Schwarz smoother if needed --*/
136  if ( lvl < param->Schwarz_levels ) {
137  mgl[lvl].Schwarz.A = fasp_dcsr_sympat(&mgl[lvl].A);
138  fasp_dcsr_shift(&(mgl[lvl].Schwarz.A), 1);
139  fasp_Schwarz_setup(&mgl[lvl].Schwarz, &swzparam);
140  }
141 
142  /*-- Coarsening and form the structure of interpolation --*/
143  status = fasp_amg_coarsening_rs(&mgl[lvl].A, &vertices, &mgl[lvl].P,
144  &Scouple, param);
145 
146  // Check 1: Did coarsening step succeeded?
147  if ( status < 0 ) {
148  // When error happens, stop at the current multigrid level!
149  if ( prtlvl > PRINT_MIN ) {
150  printf("### WARNING: Could not find any C-variables!\n");
151  printf("### WARNING: RS coarsening on level-%d failed!\n", lvl);
152  }
153  status = FASP_SUCCESS; break;
154  }
155 
156  // Check 2: Is coarse sparse too small?
157  if ( mgl[lvl].P.col < MIN_CDOF ) break;
158 
159  // Check 3: Does this coarsening step too aggressive?
160  if ( mgl[lvl].P.row > mgl[lvl].P.col * 10.0 ) {
161  if ( prtlvl > PRINT_MIN ) {
162  printf("### WARNING: Coarsening might be too aggressive!\n");
163  printf("### WARNING: Fine level = %d, coarse level = %d. Discard!\n",
164  mgl[lvl].P.row, mgl[lvl].P.col);
165  }
166  break;
167  }
168 
169  /*-- Perform aggressive coarsening only up to the specified level --*/
170  if ( mgl[lvl].P.col*1.5 > mgl[lvl].A.row ) param->coarsening_type = COARSE_RS;
171  if ( lvl == param->aggressive_level ) param->coarsening_type = COARSE_RS;
172 
173  /*-- Store the C/F marker --*/
174  {
175  INT size = mgl[lvl].A.row;
176  mgl[lvl].cfmark = fasp_ivec_create(size);
177  memcpy(mgl[lvl].cfmark.val, vertices.val, size*sizeof(INT));
178  }
179 
180  /*-- Form interpolation --*/
181  fasp_amg_interp(&mgl[lvl].A, &vertices, &mgl[lvl].P, &Scouple, param);
182 
183  /*-- Form coarse level matrix: two RAP routines available! --*/
184  fasp_dcsr_trans(&mgl[lvl].P, &mgl[lvl].R);
185 
186  fasp_blas_dcsr_rap(&mgl[lvl].R, &mgl[lvl].A, &mgl[lvl].P, &mgl[lvl+1].A);
187 
188  /*-- Clean up Scouple generated in coarsening --*/
189  fasp_mem_free(Scouple.IA);
190  fasp_mem_free(Scouple.JA);
191 
192  // Check 4: Is the coarse matrix too dense?
193  if ( mgl[lvl].A.nnz / mgl[lvl].A.row > mgl[lvl].A.col * 0.2 ) {
194  if ( prtlvl > PRINT_MIN ) {
195  printf("### WARNING: Coarse matrix is too dense!\n");
196  printf("### WARNING: m = n = %d, nnz = %d!\n",
197  mgl[lvl].A.col, mgl[lvl].A.nnz);
198  }
199  break;
200  }
201 
202  ++lvl;
203 
204 #if DIAGONAL_PREF
205  // reorder each row to make diagonal appear first
206  fasp_dcsr_diagpref(&mgl[lvl].A);
207 #endif
208 
209  }
210 
211  // Setup coarse level systems for direct solvers
212  switch (csolver) {
213 
214 #if WITH_MUMPS
215  case SOLVER_MUMPS: {
216  // Setup MUMPS direct solver on the coarsest level
217  mgl[lvl].mumps.job = 1;
218  fasp_solver_mumps_steps(&mgl[lvl].A, &mgl[lvl].b, &mgl[lvl].x, &mgl[lvl].mumps);
219  break;
220  }
221 #endif
222 
223 #if WITH_UMFPACK
224  case SOLVER_UMFPACK: {
225  // Need to sort the matrix A for UMFPACK to work
226  dCSRmat Ac_tran;
227  Ac_tran = fasp_dcsr_create(mgl[lvl].A.row, mgl[lvl].A.col, mgl[lvl].A.nnz);
228  fasp_dcsr_transz(&mgl[lvl].A, NULL, &Ac_tran);
229  // It is equivalent to do transpose and then sort
230  // fasp_dcsr_trans(&mgl[lvl].A, &Ac_tran);
231  // fasp_dcsr_sort(&Ac_tran);
232  fasp_dcsr_cp(&Ac_tran, &mgl[lvl].A);
233  fasp_dcsr_free(&Ac_tran);
234  mgl[lvl].Numeric = fasp_umfpack_factorize(&mgl[lvl].A, 0);
235  break;
236  }
237 #endif
238 
239 #if WITH_PARDISO
240  case SOLVER_PARDISO: {
241  fasp_dcsr_sort(&mgl[lvl].A);
242  fasp_pardiso_factorize(&mgl[lvl].A, &mgl[lvl].pdata, 0);
243  break;
244  }
245 #endif
246 
247  default:
248  // Do nothing!
249  break;
250  }
251 
252  // setup total level number and current level
253  mgl[0].num_levels = max_lvls = lvl+1;
254  mgl[0].w = fasp_dvec_create(m);
255 
256  for ( lvl = 1; lvl < max_lvls; ++lvl ) {
257  INT mm = mgl[lvl].A.row;
258  mgl[lvl].num_levels = max_lvls;
259  mgl[lvl].b = fasp_dvec_create(mm);
260  mgl[lvl].x = fasp_dvec_create(mm);
261 
262  mgl[lvl].cycle_type = cycle_type; // initialize cycle type!
263  mgl[lvl].ILU_levels = param->ILU_levels - lvl; // initialize ILU levels!
264  mgl[lvl].Schwarz_levels = param->Schwarz_levels -lvl; // initialize Schwarz!
265 
266  // allocate work arrays for the solve phase
267  if ( cycle_type == NL_AMLI_CYCLE )
268  mgl[lvl].w = fasp_dvec_create(3*mm);
269  else
270  mgl[lvl].w = fasp_dvec_create(2*mm);
271  }
272 
273  fasp_ivec_free(&vertices);
274 
275  if ( prtlvl > PRINT_NONE ) {
276  fasp_gettime(&setup_end);
277  print_amgcomplexity(mgl, prtlvl);
278  print_cputime("Classical AMG setup", setup_end - setup_start);
279  }
280 
281 #if DEBUG_MODE > 0
282  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
283 #endif
284 
285  return status;
286 }
287 
288 /*---------------------------------*/
289 /*-- End of File --*/
290 /*---------------------------------*/
dvector fasp_dvec_create(const INT m)
Create dvector data space of REAL type.
Definition: vec.c:56
void fasp_blas_dcsr_rap(dCSRmat *R, dCSRmat *A, dCSRmat *P, dCSRmat *RAP)
Triple sparse matrix multiplication B=R*A*P.
Definition: blas_csr.c:866
dCSRmat fasp_dcsr_create(const INT m, const INT n, const INT nnz)
Create CSR sparse matrix data memory space.
Definition: sparse_csr.c:34
Parameters for AMG solver.
Definition: fasp.h:583
void fasp_dcsr_diagpref(dCSRmat *A)
Re-order the column and data arrays of a CSR matrix, so that the first entry in each row is the diago...
Definition: sparse_csr.c:553
REAL tentative_smooth
relaxation parameter for smoothing the tentative prolongation
Definition: fasp.h:676
SHORT coarsening_type
coarsening type
Definition: fasp.h:643
Schwarz_data Schwarz
data of Schwarz smoother
Definition: fasp.h:778
INT coarse_dof
max number of coarsest level DOF
Definition: fasp.h:601
#define REAL
Definition: fasp.h:67
SHORT amli_degree
degree of the polynomial used by AMLI cycle
Definition: fasp.h:634
Vector with n entries of INT type.
Definition: fasp.h:356
void fasp_amg_amli_coef(const REAL lambda_max, const REAL lambda_min, const INT degree, REAL *coef)
Compute the coefficients of the polynomial used by AMLI-cycle.
Definition: amlirecur.c:706
SHORT num_levels
number of levels in use <= max_levels
Definition: fasp.h:730
SHORT ILU_type
ILU type for decomposition.
Definition: fasp.h:380
#define COARSE_AC
Definition: fasp_const.h:200
INT fasp_Schwarz_setup(Schwarz_data *Schwarz, Schwarz_param *param)
Setup phase for the Schwarz methods.
INT fasp_dcsr_trans(dCSRmat *A, dCSRmat *AT)
Find transpose of dCSRmat matrix A.
Definition: sparse_csr.c:826
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
INT Schwarz_levels
number of levels use Schwarz smoother
Definition: fasp.h:775
REAL ILU_relax
add the sum of dropped elements to diagonal element in proportion relax
Definition: fasp.h:389
SHORT fasp_amg_setup_rs(AMG_data *mgl, AMG_param *param)
Setup phase of Ruge and Stuben's classic AMG.
Definition: amg_setup_rs.c:47
INT Schwarz_maxlvl
maximal level for constructing the blocks
Definition: fasp.h:443
INT ILU_lfil
level of fill-in for ILUk
Definition: fasp.h:383
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
INT aggressive_level
number of levels use aggressive coarsening
Definition: fasp.h:661
#define INT
Definition: fasp.h:64
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:193
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
INT Schwarz_levels
number of levels use Schwarz smoother
Definition: fasp.h:700
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void fasp_ivec_free(ivector *u)
Free vector data space of INT type.
Definition: vec.c:159
dvector x
pointer to the iterative solution at level level_num
Definition: fasp.h:747
void fasp_dcsr_cp(dCSRmat *A, dCSRmat *B)
copy a dCSRmat to a new one B=A
Definition: sparse_csr.c:723
void fasp_dcsr_transz(dCSRmat *A, INT *p, dCSRmat *AT)
Generalized transpose of A: (n x m) matrix given in dCSRmat format.
Definition: sparse_csr.c:1366
void fasp_dcsr_shift(dCSRmat *A, INT offset)
Re-index a REAL matrix in CSR format to make the index starting from 0 or 1.
Definition: sparse_csr.c:1085
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 fasp_amg_coarsening_rs(dCSRmat *A, ivector *vertices, dCSRmat *P, iCSRmat *S, AMG_param *param)
Standard and aggressive coarsening schemes.
Definition: coarsening_rs.c:61
SHORT ILU_levels
number of levels use ILU smoother
Definition: fasp.h:682
SHORT Schwarz_type
type for Schwarz method
Definition: fasp.h:440
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:190
SHORT max_levels
max number of levels of AMG
Definition: fasp.h:598
REAL * amli_coef
coefficients of the polynomial used by AMLI cycle
Definition: fasp.h:637
void fasp_gettime(REAL *time)
Get system time.
Definition: timing.c:28
#define SOLVER_MUMPS
Definition: fasp_const.h:125
void print_amgcomplexity(AMG_data *mgl, const SHORT prtlvl)
Print complexities of AMG method.
Definition: message.c:79
INT Schwarz_mmsize
maximal size of blocks
Definition: fasp.h:446
REAL ILU_relax
relaxation for ILUs
Definition: fasp.h:694
void fasp_amg_interp(dCSRmat *A, ivector *vertices, dCSRmat *P, iCSRmat *S, AMG_param *param)
Generate interpolation operator P.
Definition: interpolation.c:48
INT * val
actual vector entries
Definition: fasp.h:362
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
REAL ILU_droptol
drop tolerance for ILUt
Definition: fasp.h:386
INT row
row number of matrix A, m
Definition: fasp.h:151
void fasp_dcsr_sort(dCSRmat *A)
Sort each row of A in ascending order w.r.t. column indices.
Definition: sparse_csr.c:358
INT Schwarz_blksolver
type of Schwarz block solver
Definition: fasp.h:449
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
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
dCSRmat P
prolongation operator at level level_num
Definition: fasp.h:741
INT ILU_levels
number of levels use ILU smoother
Definition: fasp.h:761
SHORT ILU_type
ILU type for smoothing.
Definition: fasp.h:685
SHORT cycle_type
type of AMG cycle
Definition: fasp.h:604
dCSRmat A
pointer to the matrix
Definition: fasp.h:510
Parameters for Schwarz method.
Definition: fasp.h:434
INT cycle_type
cycle type
Definition: fasp.h:787
SHORT fasp_ilu_dcsr_setup(dCSRmat *A, ILU_data *iludata, ILU_param *iluparam)
Get ILU decomposition of a CSR matrix A.
Definition: ilu_setup_csr.c:50
Parameters for ILU.
Definition: fasp.h:374
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define AMLI_CYCLE
Definition: fasp_const.h:177
Mumps_data mumps
data for MUMPS
Definition: fasp.h:784
INT Schwarz_mmsize
maximal block size
Definition: fasp.h:703
#define COARSE_RS
Definition of coarsening types.
Definition: fasp_const.h:197
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
SHORT print_level
print level for AMG
Definition: fasp.h:589
Data for AMG solvers.
Definition: fasp.h:722
#define SOLVER_UMFPACK
Definition: fasp_const.h:124
void * Numeric
pointer to the numerical factorization from UMFPACK
Definition: fasp.h:752
INT Schwarz_maxlvl
maximal levels
Definition: fasp.h:706
#define NL_AMLI_CYCLE
Definition: fasp_const.h:178
REAL ILU_droptol
drop tolerance for ILUt
Definition: fasp.h:691
ivector fasp_ivec_create(const INT m)
Create vector data space of INT type.
Definition: vec.c:78
INT job
work for MUMPS
Definition: fasp.h:467
INT Schwarz_type
type of Schwarz method
Definition: fasp.h:709
INT ILU_lfil
level of fill-in for ILUs and ILUk
Definition: fasp.h:688
#define MIN_CDOF
Definition: fasp_const.h:242
dvector w
Temporary work space.
Definition: fasp.h:781
Sparse matrix of INT type in CSR format.
Definition: fasp.h:178
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: sparse_csr.c:166
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72
#define PRINT_MIN
Definition: fasp_const.h:80
dCSRmat fasp_dcsr_sympat(dCSRmat *A)
Get symmetric part of a dCSRmat matrix.
Definition: sparse_csr.c:1232
SHORT print_level
print level
Definition: fasp.h:377