Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
amg_setup_sa.c
Go to the documentation of this file.
1 
10 #include <math.h>
11 #include <time.h>
12 
13 #ifdef _OPENMP
14 #include <omp.h>
15 #endif
16 
17 #include "fasp.h"
18 #include "fasp_functs.h"
19 #include "aggregation_csr.inl"
20 #include "aggregation_bsr.inl"
21 
22 static SHORT amg_setup_smoothP_smoothR (AMG_data *, AMG_param *);
23 static SHORT amg_setup_smoothP_unsmoothR (AMG_data *, AMG_param *);
24 static SHORT amg_setup_smoothP_smoothR_bsr (AMG_data_bsr *mgl, AMG_param *param);
25 
26 /*---------------------------------*/
27 /*-- Public Functions --*/
28 /*---------------------------------*/
29 
49  AMG_param *param)
50 {
51  SHORT status = FASP_SUCCESS;
52 
53 #if DEBUG_MODE > 0
54  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
55  printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n",
56  mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
57 #endif
58 
59 #if TRUE
60  status = amg_setup_smoothP_smoothR(mgl, param);
61 #else // smoothed P, unsmoothed R
62  status = amg_setup_smoothP_unsmoothR(mgl, param);
63 #endif
64 
65 #if DEBUG_MODE > 0
66  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
67 #endif
68 
69  return status;
70 }
71 
86  AMG_param *param)
87 {
88 #if DEBUG_MODE > 0
89  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
90 #endif
91 
92  SHORT status = amg_setup_smoothP_smoothR_bsr(mgl, param);
93 
94 #if DEBUG_MODE > 0
95  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
96 #endif
97 
98  return status;
99 }
100 
101 /*---------------------------------*/
102 /*-- Private Functions --*/
103 /*---------------------------------*/
104 
120 static SHORT amg_setup_smoothP_smoothR (AMG_data *mgl,
121  AMG_param *param)
122 {
123  const SHORT prtlvl = param->print_level;
124  const SHORT cycle_type = param->cycle_type;
125  const SHORT csolver = param->coarse_solver;
126  const SHORT min_cdof = MAX(param->coarse_dof,50);
127  const INT m = mgl[0].A.row;
128 
129  // local variables
130  SHORT max_levels = param->max_levels, lvl = 0, status = FASP_SUCCESS;
131  INT i, j;
132  REAL setup_start, setup_end;
133  ILU_param iluparam;
134  Schwarz_param swzparam;
135 
136 #if DEBUG_MODE > 0
137  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
138 #endif
139 
140  fasp_gettime(&setup_start);
141 
142  // level info (fine: 0; coarse: 1)
143  ivector *vertices = (ivector *)fasp_mem_calloc(max_levels,sizeof(ivector));
144 
145  // each elvel stores the information of the number of aggregations
146  INT *num_aggregations = (INT *)fasp_mem_calloc(max_levels,sizeof(INT));
147 
148  // each level stores the information of the strongly coupled neighbourhood
149  dCSRmat *Neighbor = (dCSRmat *)fasp_mem_calloc(max_levels,sizeof(dCSRmat));
150 
151  // each level stores the information of the tentative prolongations
152  dCSRmat *tentp = (dCSRmat *)fasp_mem_calloc(max_levels,sizeof(dCSRmat));
153 
154  // Initialize level information
155  for ( i = 0; i < max_levels; ++i ) num_aggregations[i] = 0;
156 
157  mgl[0].near_kernel_dim = 1;
158  mgl[0].near_kernel_basis = (REAL **)fasp_mem_calloc(mgl->near_kernel_dim,sizeof(REAL*));
159 
160  for ( i = 0; i < mgl->near_kernel_dim; ++i ) {
161  mgl[0].near_kernel_basis[i] = (REAL *)fasp_mem_calloc(m,sizeof(REAL));
162  for ( j = 0; j < m; ++j ) mgl[0].near_kernel_basis[i][j] = 1.0;
163  }
164 
165  // Initialize ILU parameters
166  mgl->ILU_levels = param->ILU_levels;
167  if ( param->ILU_levels > 0 ) {
168  iluparam.print_level = param->print_level;
169  iluparam.ILU_lfil = param->ILU_lfil;
170  iluparam.ILU_droptol = param->ILU_droptol;
171  iluparam.ILU_relax = param->ILU_relax;
172  iluparam.ILU_type = param->ILU_type;
173  }
174 
175  // Initialize Schwarz parameters
176  mgl->Schwarz_levels = param->Schwarz_levels;
177  if ( param->Schwarz_levels > 0 ) {
178  swzparam.Schwarz_mmsize = param->Schwarz_mmsize;
179  swzparam.Schwarz_maxlvl = param->Schwarz_maxlvl;
180  swzparam.Schwarz_type = param->Schwarz_type;
181  swzparam.Schwarz_blksolver = param->Schwarz_blksolver;
182  }
183 
184  // Initialize AMLI coefficients
185  if ( cycle_type == AMLI_CYCLE ) {
186  const INT amlideg = param->amli_degree;
187  param->amli_coef = (REAL *)fasp_mem_calloc(amlideg+1,sizeof(REAL));
188  REAL lambda_max = 2.0, lambda_min = lambda_max/4;
189  fasp_amg_amli_coef(lambda_max, lambda_min, amlideg, param->amli_coef);
190  }
191 
192 #if DIAGONAL_PREF
193  fasp_dcsr_diagpref(&mgl[0].A); // reorder each row to make diagonal appear first
194 #endif
195 
196  /*----------------------------*/
197  /*--- checking aggregation ---*/
198  /*----------------------------*/
199  if ( param->aggregation_type == PAIRWISE )
200  param->pair_number = MIN(param->pair_number, max_levels);
201 
202  // Main AMG setup loop
203  while ( (mgl[lvl].A.row > min_cdof) && (lvl < max_levels-1) ) {
204 
205 #if DEBUG_MODE > 2
206  printf("### DEBUG: level = %d, row = %d, nnz = %d\n",
207  lvl, mgl[lvl].A.row, mgl[lvl].A.nnz);
208 #endif
209 
210  /*-- setup ILU decomposition if necessary */
211  if ( lvl < param->ILU_levels ) {
212  status = fasp_ilu_dcsr_setup(&mgl[lvl].A, &mgl[lvl].LU, &iluparam);
213  if ( status < 0 ) {
214  if ( prtlvl > PRINT_MIN ) {
215  printf("### WARNING: ILU setup on level-%d failed!\n", lvl);
216  printf("### WARNING: Disable ILU for level >= %d.\n", lvl);
217  }
218  param->ILU_levels = lvl;
219  }
220  }
221 
222  /* -- setup Schwarz smoother if necessary */
223  if ( lvl < param->Schwarz_levels ) {
224  mgl[lvl].Schwarz.A = fasp_dcsr_sympat(&mgl[lvl].A);
225  fasp_dcsr_shift(&(mgl[lvl].Schwarz.A), 1);
226  fasp_Schwarz_setup(&mgl[lvl].Schwarz, &swzparam);
227  }
228 
229  /*-- Aggregation --*/
230  status = aggregation_vmb(&mgl[lvl].A, &vertices[lvl], param, lvl+1,
231  &Neighbor[lvl], &num_aggregations[lvl]);
232 
233  // Check 1: Did coarsening step succeed?
234  if ( status < 0 ) {
235  // When error happens, stop at the current multigrid level!
236  if ( prtlvl > PRINT_MIN ) {
237  printf("### WARNING: Forming aggregates on level-%d failed!\n", lvl);
238  }
239  status = FASP_SUCCESS; break;
240  }
241 
242  /* -- Form Tentative prolongation --*/
243  form_tentative_p(&vertices[lvl], &tentp[lvl], mgl[0].near_kernel_basis,
244  lvl+1, num_aggregations[lvl]);
245 
246  /* -- Form smoothed prolongation -- */
247  smooth_agg(&mgl[lvl].A, &tentp[lvl], &mgl[lvl].P, param, lvl+1,
248  &Neighbor[lvl]);
249 
250  // Check 2: Is coarse sparse too small?
251  if ( mgl[lvl].P.col < MIN_CDOF ) break;
252 
253  // Check 3: Does this coarsening step too aggressive?
254  if ( mgl[lvl].P.row > mgl[lvl].P.col * MAX_CRATE ) {
255  if ( prtlvl > PRINT_MIN ) {
256  printf("### WARNING: Coarsening might be too aggressive!\n");
257  printf("### WARNING: Fine level = %d, coarse level = %d. Discard!\n",
258  mgl[lvl].P.row, mgl[lvl].P.col);
259  }
260  break;
261  }
262 
263  // Check 4: Is this coarsening ratio too small?
264  if ( (REAL)mgl[lvl].P.col > mgl[lvl].P.row * MIN_CRATE ) {
265  if ( prtlvl > PRINT_MIN ) {
266  printf("### WARNING: Coarsening rate is too small!\n");
267  printf("### WARNING: Fine level = %d, coarse level = %d. Discard!\n",
268  mgl[lvl].P.row, mgl[lvl].P.col);
269  }
270  break;
271  }
272 
273  /*-- Form restriction --*/
274  fasp_dcsr_trans(&mgl[lvl].P, &mgl[lvl].R);
275 
276  /*-- Form coarse level stiffness matrix --*/
277  fasp_blas_dcsr_rap(&mgl[lvl].R, &mgl[lvl].A, &mgl[lvl].P, &mgl[lvl+1].A);
278 
279  fasp_dcsr_free(&Neighbor[lvl]);
280  fasp_dcsr_free(&tentp[lvl]);
281  fasp_ivec_free(&vertices[lvl]);
282 
283  ++lvl;
284 
285 #if DIAGONAL_PREF
286  // reorder each row to make diagonal appear first
287  fasp_dcsr_diagpref(&mgl[lvl].A);
288 #endif
289 
290  }
291 
292  // Setup coarse level systems for direct solvers
293  switch (csolver) {
294 
295 #if WITH_MUMPS
296  case SOLVER_MUMPS: {
297  // Setup MUMPS direct solver on the coarsest level
298  mgl[lvl].mumps.job = 1;
299  fasp_solver_mumps_steps(&mgl[lvl].A, &mgl[lvl].b, &mgl[lvl].x, &mgl[lvl].mumps);
300  break;
301  }
302 #endif
303 
304 #if WITH_UMFPACK
305  case SOLVER_UMFPACK: {
306  // Need to sort the matrix A for UMFPACK to work
307  dCSRmat Ac_tran;
308  Ac_tran = fasp_dcsr_create(mgl[lvl].A.row, mgl[lvl].A.col, mgl[lvl].A.nnz);
309  fasp_dcsr_transz(&mgl[lvl].A, NULL, &Ac_tran);
310  // It is equivalent to do transpose and then sort
311  // fasp_dcsr_trans(&mgl[lvl].A, &Ac_tran);
312  // fasp_dcsr_sort(&Ac_tran);
313  fasp_dcsr_cp(&Ac_tran, &mgl[lvl].A);
314  fasp_dcsr_free(&Ac_tran);
315  mgl[lvl].Numeric = fasp_umfpack_factorize(&mgl[lvl].A, 0);
316  break;
317  }
318 #endif
319 
320 #if WITH_PARDISO
321  case SOLVER_PARDISO: {
322  fasp_dcsr_sort(&mgl[lvl].A);
323  fasp_pardiso_factorize(&mgl[lvl].A, &mgl[lvl].pdata, 0);
324  break;
325  }
326 #endif
327 
328  default:
329  // Do nothing!
330  break;
331  }
332 
333  // setup total level number and current level
334  mgl[0].num_levels = max_levels = lvl+1;
335  mgl[0].w = fasp_dvec_create(m);
336 
337  for ( lvl = 1; lvl < max_levels; ++lvl) {
338  INT mm = mgl[lvl].A.row;
339  mgl[lvl].num_levels = max_levels;
340  mgl[lvl].b = fasp_dvec_create(mm);
341  mgl[lvl].x = fasp_dvec_create(mm);
342 
343  mgl[lvl].cycle_type = cycle_type; // initialize cycle type!
344  mgl[lvl].ILU_levels = param->ILU_levels - lvl; // initialize ILU levels!
345  mgl[lvl].Schwarz_levels = param->Schwarz_levels -lvl; // initialize Schwarz!
346 
347  if ( cycle_type == NL_AMLI_CYCLE )
348  mgl[lvl].w = fasp_dvec_create(3*mm);
349  else
350  mgl[lvl].w = fasp_dvec_create(2*mm);
351  }
352 
353  if ( prtlvl > PRINT_NONE ) {
354  fasp_gettime(&setup_end);
355  print_amgcomplexity(mgl,prtlvl);
356  print_cputime("Smoothed aggregation setup", setup_end - setup_start);
357  }
358 
359  fasp_mem_free(vertices);
360  fasp_mem_free(num_aggregations);
361  fasp_mem_free(Neighbor);
362  fasp_mem_free(tentp);
363 
364 #if DEBUG_MODE > 0
365  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
366 #endif
367 
368  return status;
369 }
370 
386 static SHORT amg_setup_smoothP_unsmoothR (AMG_data *mgl,
387  AMG_param *param)
388 {
389  const SHORT prtlvl = param->print_level;
390  const SHORT cycle_type = param->cycle_type;
391  const SHORT csolver = param->coarse_solver;
392  const SHORT min_cdof = MAX(param->coarse_dof,50);
393  const INT m = mgl[0].A.row;
394 
395  // local variables
396  SHORT max_levels = param->max_levels, lvl = 0, status = FASP_SUCCESS;
397  INT i, j;
398  REAL setup_start, setup_end;
399  ILU_param iluparam;
400  Schwarz_param swzparam;
401 
402 #if DEBUG_MODE > 0
403  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
404 #endif
405 
406  fasp_gettime(&setup_start);
407 
408  // level info (fine: 0; coarse: 1)
409  ivector *vertices = (ivector *)fasp_mem_calloc(max_levels,sizeof(ivector));
410 
411  // each level stores the information of the number of aggregations
412  INT *num_aggregations = (INT *)fasp_mem_calloc(max_levels,sizeof(INT));
413 
414  // each level stores the information of the strongly coupled neighbourhood
415  dCSRmat *Neighbor = (dCSRmat *)fasp_mem_calloc(max_levels,sizeof(dCSRmat));
416 
417  // each level stores the information of the tentative prolongations
418  dCSRmat *tentp = (dCSRmat *)fasp_mem_calloc(max_levels,sizeof(dCSRmat));
419  dCSRmat *tentr = (dCSRmat *)fasp_mem_calloc(max_levels,sizeof(dCSRmat));
420 
421  for ( i = 0; i < max_levels; ++i ) num_aggregations[i] = 0;
422 
423  mgl[0].near_kernel_dim = 1;
424 
425  mgl[0].near_kernel_basis = (REAL **)fasp_mem_calloc(mgl->near_kernel_dim,sizeof(REAL*));
426 
427  for ( i = 0; i < mgl->near_kernel_dim; ++i ) {
428  mgl[0].near_kernel_basis[i] = (REAL *)fasp_mem_calloc(m,sizeof(REAL));
429  for ( j = 0; j < m; ++j ) mgl[0].near_kernel_basis[i][j] = 1.0;
430  }
431 
432  // Initialize ILU parameters
433  if ( param->ILU_levels > 0 ) {
434  iluparam.print_level = param->print_level;
435  iluparam.ILU_lfil = param->ILU_lfil;
436  iluparam.ILU_droptol = param->ILU_droptol;
437  iluparam.ILU_relax = param->ILU_relax;
438  iluparam.ILU_type = param->ILU_type;
439  }
440 
441  // Initialize Schwarz parameters
442  mgl->Schwarz_levels = param->Schwarz_levels;
443  if ( param->Schwarz_levels > 0 ) {
444  swzparam.Schwarz_mmsize = param->Schwarz_mmsize;
445  swzparam.Schwarz_maxlvl = param->Schwarz_maxlvl;
446  swzparam.Schwarz_type = param->Schwarz_type;
447  swzparam.Schwarz_blksolver = param->Schwarz_blksolver;
448  }
449 
450  // Initialize AMLI coefficients
451  if ( cycle_type == AMLI_CYCLE ) {
452  const INT amlideg = param->amli_degree;
453  param->amli_coef = (REAL *)fasp_mem_calloc(amlideg+1,sizeof(REAL));
454  REAL lambda_max = 2.0, lambda_min = lambda_max/4;
455  fasp_amg_amli_coef(lambda_max, lambda_min, amlideg, param->amli_coef);
456  }
457 
458  // Main AMG setup loop
459  while ( (mgl[lvl].A.row > min_cdof) && (lvl < max_levels-1) ) {
460 
461  /*-- setup ILU decomposition if necessary */
462  if ( lvl < param->ILU_levels ) {
463  status = fasp_ilu_dcsr_setup(&mgl[lvl].A, &mgl[lvl].LU, &iluparam);
464  if ( status < 0 ) {
465  if ( prtlvl > PRINT_MIN ) {
466  printf("### WARNING: ILU setup on level-%d failed!\n", lvl);
467  printf("### WARNING: Disable ILU for level >= %d.\n", lvl);
468  }
469  param->ILU_levels = lvl;
470  }
471  }
472 
473  /* -- setup Schwarz smoother if necessary */
474  if ( lvl < param->Schwarz_levels ) {
475  mgl[lvl].Schwarz.A = fasp_dcsr_sympat(&mgl[lvl].A);
476  fasp_dcsr_shift(&(mgl[lvl].Schwarz.A), 1);
477  fasp_Schwarz_setup(&mgl[lvl].Schwarz, &swzparam);
478  }
479 
480  /*-- Aggregation --*/
481  status = aggregation_vmb(&mgl[lvl].A, &vertices[lvl], param, lvl+1,
482  &Neighbor[lvl], &num_aggregations[lvl]);
483 
484  // Check 1: Did coarsening step succeeded?
485  if ( status < 0 ) {
486  // When error happens, stop at the current multigrid level!
487  if ( prtlvl > PRINT_MIN ) {
488  printf("### WARNING: Forming aggregates on level-%d failed!\n", lvl);
489  }
490  status = FASP_SUCCESS; break;
491  }
492 
493  /* -- Form Tentative prolongation --*/
494  form_tentative_p(&vertices[lvl], &tentp[lvl], mgl[0].near_kernel_basis,
495  lvl+1, num_aggregations[lvl]);
496 
497  /* -- Form smoothed prolongation -- */
498  smooth_agg(&mgl[lvl].A, &tentp[lvl], &mgl[lvl].P, param, lvl+1,
499  &Neighbor[lvl]);
500 
501  // Check 2: Is coarse sparse too small?
502  if ( mgl[lvl].P.col < MIN_CDOF ) break;
503 
504  // Check 3: Does this coarsening step too aggressive?
505  if ( mgl[lvl].P.row > mgl[lvl].P.col * MAX_CRATE ) {
506  if ( prtlvl > PRINT_MIN ) {
507  printf("### WARNING: Coarsening might be too aggressive!\n");
508  printf("### WARNING: Fine level = %d, coarse level = %d. Discard!\n",
509  mgl[lvl].P.row, mgl[lvl].P.col);
510  }
511  break;
512  }
513 
514  // Check 4: Is this coarsening ratio too small?
515  if ( (REAL)mgl[lvl].P.col > mgl[lvl].P.row * MIN_CRATE ) {
516  if ( prtlvl > PRINT_MIN ) {
517  printf("### WARNING: Coarsening rate is too small!\n");
518  printf("### WARNING: Fine level = %d, coarse level = %d. Discard!\n",
519  mgl[lvl].P.row, mgl[lvl].P.col);
520  }
521  break;
522  }
523 
524  /*-- Form restriction --*/
525  fasp_dcsr_trans(&mgl[lvl].P, &mgl[lvl].R);
526  fasp_dcsr_trans(&tentp[lvl], &tentr[lvl]);
527 
528  /*-- Form coarse level stiffness matrix --*/
529  fasp_blas_dcsr_rap_agg(&tentr[lvl], &mgl[lvl].A, &tentp[lvl], &mgl[lvl+1].A);
530 
531  fasp_dcsr_free(&Neighbor[lvl]);
532  fasp_dcsr_free(&tentp[lvl]);
533  fasp_ivec_free(&vertices[lvl]);
534 
535  ++lvl;
536  }
537 
538  // Setup coarse level systems for direct solvers
539  switch (csolver) {
540 
541 #if WITH_MUMPS
542  case SOLVER_MUMPS: {
543  // Setup MUMPS direct solver on the coarsest level
544  mgl[lvl].mumps.job = 1;
545  fasp_solver_mumps_steps(&mgl[lvl].A, &mgl[lvl].b, &mgl[lvl].x, &mgl[lvl].mumps);
546  break;
547  }
548 #endif
549 
550 #if WITH_UMFPACK
551  case SOLVER_UMFPACK: {
552  // Need to sort the matrix A for UMFPACK to work
553  dCSRmat Ac_tran;
554  Ac_tran = fasp_dcsr_create(mgl[lvl].A.row, mgl[lvl].A.col, mgl[lvl].A.nnz);
555  fasp_dcsr_transz(&mgl[lvl].A, NULL, &Ac_tran);
556  // It is equivalent to do transpose and then sort
557  // fasp_dcsr_trans(&mgl[lvl].A, &Ac_tran);
558  // fasp_dcsr_sort(&Ac_tran);
559  fasp_dcsr_cp(&Ac_tran, &mgl[lvl].A);
560  fasp_dcsr_free(&Ac_tran);
561  mgl[lvl].Numeric = fasp_umfpack_factorize(&mgl[lvl].A, 0);
562  break;
563  }
564 #endif
565 
566 #if WITH_PARDISO
567  case SOLVER_PARDISO: {
568  fasp_dcsr_sort(&mgl[lvl].A);
569  fasp_pardiso_factorize(&mgl[lvl].A, &mgl[lvl].pdata, 0);
570  break;
571  }
572 #endif
573 
574  default:
575  // Do nothing!
576  break;
577  }
578 
579  // setup total level number and current level
580  mgl[0].num_levels = max_levels = lvl+1;
581  mgl[0].w = fasp_dvec_create(m);
582 
583  for ( lvl = 1; lvl < max_levels; ++lvl) {
584  INT mm = mgl[lvl].A.row;
585  mgl[lvl].num_levels = max_levels;
586  mgl[lvl].b = fasp_dvec_create(mm);
587  mgl[lvl].x = fasp_dvec_create(mm);
588 
589  mgl[lvl].cycle_type = cycle_type; // initialize cycle type!
590  mgl[lvl].ILU_levels = param->ILU_levels - lvl; // initialize ILU levels!
591  mgl[lvl].Schwarz_levels = param->Schwarz_levels -lvl; // initialize Schwarz!
592 
593  if ( cycle_type == NL_AMLI_CYCLE )
594  mgl[lvl].w = fasp_dvec_create(3*mm);
595  else
596  mgl[lvl].w = fasp_dvec_create(2*mm);
597  }
598 
599  if ( prtlvl > PRINT_NONE ) {
600  fasp_gettime(&setup_end);
601  print_amgcomplexity(mgl,prtlvl);
602  print_cputime("Smoothed aggregation 1/2 setup", setup_end - setup_start);
603  }
604 
605  fasp_mem_free(vertices);
606  fasp_mem_free(num_aggregations);
607  fasp_mem_free(Neighbor);
608  fasp_mem_free(tentp);
609  fasp_mem_free(tentr);
610 
611 #if DEBUG_MODE > 0
612  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
613 #endif
614 
615  return status;
616 }
617 
618 
635 static SHORT amg_setup_smoothP_smoothR_bsr (AMG_data_bsr *mgl,
636  AMG_param *param)
637 {
638  const SHORT prtlvl = param->print_level;
639  const SHORT csolver = param->coarse_solver;
640  const SHORT min_cdof = MAX(param->coarse_dof,50);
641  const INT m = mgl[0].A.ROW;
642  const INT nb = mgl[0].A.nb;
643 
644  ILU_param iluparam;
645  SHORT max_levels=param->max_levels;
646  SHORT i, lvl=0, status=FASP_SUCCESS;
647  REAL setup_start, setup_end;
648 
649  dCSRmat temp1, temp2;
650 
651 #if DEBUG_MODE > 0
652  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
653  printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n",
654  mgl[0].A.ROW, mgl[0].A.COL, mgl[0].A.NNZ);
655 #endif
656 
657  fasp_gettime(&setup_start);
658 
659  /*-----------------------*/
660  /*--local working array--*/
661  /*-----------------------*/
662 
663  // level info (fine: 0; coarse: 1)
664  ivector *vertices = (ivector *)fasp_mem_calloc(max_levels, sizeof(ivector));
665 
666  // each level stores the information of the number of aggregations
667  INT *num_aggs = (INT *)fasp_mem_calloc(max_levels, sizeof(INT));
668 
669  // each level stores the information of the strongly coupled neighbourhood
670  dCSRmat *Neighbor = (dCSRmat *)fasp_mem_calloc(max_levels, sizeof(dCSRmat));
671 
672  // each level stores the information of the tentative prolongations
673  dBSRmat *tentp = (dBSRmat *)fasp_mem_calloc(max_levels,sizeof(dBSRmat));
674 
675  for ( i=0; i<max_levels; ++i ) num_aggs[i] = 0;
676 
677  /*-----------------------*/
678  /*-- setup null spaces --*/
679  /*-----------------------*/
680 
681  // null space for whole Jacobian
682  //mgl[0].near_kernel_dim = 1;
683  //mgl[0].near_kernel_basis = (REAL **)fasp_mem_calloc(mgl->near_kernel_dim, sizeof(REAL*));
684 
685  //for ( i=0; i < mgl->near_kernel_dim; ++i ) mgl[0].near_kernel_basis[i] = NULL;
686 
687  /*-----------------------*/
688  /*-- setup ILU param --*/
689  /*-----------------------*/
690 
691  // initialize ILU parameters
692  mgl->ILU_levels = param->ILU_levels;
693  if ( param->ILU_levels > 0 ) {
694  iluparam.print_level = param->print_level;
695  iluparam.ILU_lfil = param->ILU_lfil;
696  iluparam.ILU_droptol = param->ILU_droptol;
697  iluparam.ILU_relax = param->ILU_relax;
698  iluparam.ILU_type = param->ILU_type;
699  }
700 
701  /*----------------------------*/
702  /*--- checking aggregation ---*/
703  /*----------------------------*/
704 
705  if (param->aggregation_type == PAIRWISE)
706  param->pair_number = MIN(param->pair_number, max_levels);
707 
708  // Main AMG setup loop
709  while ( (mgl[lvl].A.ROW > min_cdof) && (lvl < max_levels-1) ) {
710 
711  /*-- setup ILU decomposition if necessary */
712  if ( lvl < param->ILU_levels ) {
713  status = fasp_ilu_dbsr_setup(&mgl[lvl].A, &mgl[lvl].LU, &iluparam);
714  if ( status < 0 ) {
715  if ( prtlvl > PRINT_MIN ) {
716  printf("### WARNING: ILU setup on level-%d failed!\n", lvl);
717  printf("### WARNING: Disable ILU for level >= %d.\n", lvl);
718  }
719  param->ILU_levels = lvl;
720  }
721  }
722 
723  /*-- get the diagonal inverse --*/
724  mgl[lvl].diaginv = fasp_dbsr_getdiaginv(&mgl[lvl].A);
725 
726  /*-- Aggregation --*/
727  //mgl[lvl].PP = fasp_dbsr_getblk_dcsr(&mgl[lvl].A);
728  mgl[lvl].PP = fasp_dbsr_Linfinity_dcsr(&mgl[lvl].A);
729 
730  switch ( param->aggregation_type ) {
731 
732  case VMB: // VMB aggregation
733  // Same as default
734 
735  default: // only one aggregation is tested!!! --Chensong
736 
737  status = aggregation_vmb(&mgl[lvl].PP, &vertices[lvl], param, lvl+1,
738  &Neighbor[lvl], &num_aggs[lvl]);
739 
740  /*-- Choose strength threshold adaptively --*/
741  if ( num_aggs[lvl]*4 > mgl[lvl].PP.row )
742  param->strong_coupled /= 4;
743  else if ( num_aggs[lvl]*1.25 < mgl[lvl].PP.row )
744  param->strong_coupled *= 1.5;
745 
746  break;
747  }
748 
749  if ( status < 0 ) {
750  // When error happens, force solver to use the current multigrid levels!
751  if ( prtlvl > PRINT_MIN ) {
752  printf("### WARNING: Aggregation on level-%d failed!\n", lvl);
753  }
754  status = FASP_SUCCESS; break;
755  }
756 
757  /* -- Form Tentative prolongation --*/
758  if (lvl == 0 && mgl[0].near_kernel_dim >0 ){
759  form_tentative_p_bsr1(&vertices[lvl], &tentp[lvl], &mgl[0], lvl+1,
760  num_aggs[lvl], mgl[0].near_kernel_dim,
761  mgl[0].near_kernel_basis);
762  }
763  else{
764  form_boolean_p_bsr(&vertices[lvl], &tentp[lvl], &mgl[0], lvl+1, num_aggs[lvl]);
765  }
766 
767  /* -- Smoothing -- */
768  smooth_agg_bsr(&mgl[lvl].A, &tentp[lvl], &mgl[lvl].P, param, lvl+1,
769  &Neighbor[lvl]);
770 
771  /*-- Form restriction --*/
772  fasp_dbsr_trans(&mgl[lvl].P, &mgl[lvl].R);
773 
774  /*-- Form coarse level stiffness matrix --*/
775  fasp_blas_dbsr_rap(&mgl[lvl].R, &mgl[lvl].A, &mgl[lvl].P, &mgl[lvl+1].A);
776 
777  /*-- Form extra near kernel space if needed --*/
778  if (mgl[lvl].A_nk != NULL){
779 
780  mgl[lvl+1].A_nk = (dCSRmat *)fasp_mem_calloc(1, sizeof(dCSRmat));
781  mgl[lvl+1].P_nk = (dCSRmat *)fasp_mem_calloc(1, sizeof(dCSRmat));
782  mgl[lvl+1].R_nk = (dCSRmat *)fasp_mem_calloc(1, sizeof(dCSRmat));
783 
784  temp1 = fasp_format_dbsr_dcsr(&mgl[lvl].R);
785  fasp_blas_dcsr_mxm(&temp1, mgl[lvl].P_nk, mgl[lvl+1].P_nk);
786  fasp_dcsr_trans(mgl[lvl+1].P_nk, mgl[lvl+1].R_nk);
787  temp2 = fasp_format_dbsr_dcsr(&mgl[lvl+1].A);
788  fasp_blas_dcsr_rap(mgl[lvl+1].R_nk, &temp2, mgl[lvl+1].P_nk, mgl[lvl+1].A_nk);
789  fasp_dcsr_free(&temp1);
790  fasp_dcsr_free(&temp2);
791 
792  }
793 
794  fasp_dcsr_free(&Neighbor[lvl]);
795  fasp_ivec_free(&vertices[lvl]);
796  fasp_dbsr_free(&tentp[lvl]);
797 
798  ++lvl;
799  }
800 
801  // Setup coarse level systems for direct solvers (BSR version)
802  switch (csolver) {
803 
804 #if WITH_MUMPS
805  case SOLVER_MUMPS: {
806  // Setup MUMPS direct solver on the coarsest level
807  mgl[lvl].mumps.job = 1;
808  mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
809  fasp_solver_mumps_steps(&mgl[lvl].Ac, &mgl[lvl].b, &mgl[lvl].x, &mgl[lvl].mumps);
810  break;
811  }
812 #endif
813 
814 #if WITH_UMFPACK
815  case SOLVER_UMFPACK: {
816  // Need to sort the matrix A for UMFPACK to work
817  mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
818  dCSRmat Ac_tran;
819  fasp_dcsr_trans(&mgl[lvl].Ac, &Ac_tran);
820  fasp_dcsr_sort(&Ac_tran);
821  fasp_dcsr_cp(&Ac_tran, &mgl[lvl].Ac);
822  fasp_dcsr_free(&Ac_tran);
823  mgl[lvl].Numeric = fasp_umfpack_factorize(&mgl[lvl].Ac, 0);
824  break;
825  }
826 #endif
827 
828 #if WITH_SuperLU
829  case SOLVER_SUPERLU: {
830  /* Setup SuperLU direct solver on the coarsest level */
831  mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
832  }
833 #endif
834 
835 #if WITH_PARDISO
836  case SOLVER_PARDISO: {
837  mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
838  fasp_dcsr_sort(&mgl[lvl].Ac);
839  fasp_pardiso_factorize(&mgl[lvl].Ac, &mgl[lvl].pdata, 0);
840  break;
841  }
842 #endif
843 
844  default:
845  // Do nothing!
846  break;
847  }
848 
849  // setup total level number and current level
850  mgl[0].num_levels = max_levels = lvl+1;
851  mgl[0].w = fasp_dvec_create(3*m*nb);
852 
853  if (mgl[0].A_nk != NULL){
854 
855 #if WITH_UMFPACK
856  // Need to sort the matrix A_nk for UMFPACK
857  fasp_dcsr_trans(mgl[0].A_nk, &temp1);
858  fasp_dcsr_sort(&temp1);
859  fasp_dcsr_cp(&temp1, mgl[0].A_nk);
860  fasp_dcsr_free(&temp1);
861  mgl[0].Numeric = fasp_umfpack_factorize(mgl[0].A_nk, 0);
862 #endif
863 
864  }
865 
866  for ( lvl = 1; lvl < max_levels; lvl++ ) {
867  const INT mm = mgl[lvl].A.ROW*nb;
868  mgl[lvl].num_levels = max_levels;
869  mgl[lvl].b = fasp_dvec_create(mm);
870  mgl[lvl].x = fasp_dvec_create(mm);
871  mgl[lvl].w = fasp_dvec_create(3*mm);
872  mgl[lvl].ILU_levels = param->ILU_levels - lvl; // initialize ILU levels!
873 
874  if (mgl[lvl].A_nk != NULL){
875 
876 #if WITH_UMFPACK
877  // Need to sort the matrix A_nk for UMFPACK
878  fasp_dcsr_trans(mgl[lvl].A_nk, &temp1);
879  fasp_dcsr_sort(&temp1);
880  fasp_dcsr_cp(&temp1, mgl[lvl].A_nk);
881  fasp_dcsr_free(&temp1);
882  mgl[lvl].Numeric = fasp_umfpack_factorize(mgl[lvl].A_nk, 0);
883 #endif
884 
885  }
886 
887  }
888 
889  if ( prtlvl > PRINT_NONE ) {
890  fasp_gettime(&setup_end);
891  print_amgcomplexity_bsr(mgl,prtlvl);
892  print_cputime("Smoothed aggregation (BSR) setup", setup_end - setup_start);
893  }
894 
895  fasp_mem_free(vertices);
896  fasp_mem_free(num_aggs);
897  fasp_mem_free(Neighbor);
898 
899 #if DEBUG_MODE > 0
900  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
901 #endif
902 
903  return status;
904 }
905 
906 /*---------------------------------*/
907 /*-- End of File --*/
908 /*---------------------------------*/
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
dvector x
pointer to the iterative solution at level level_num
Definition: fasp_block.h:219
Parameters for AMG solver.
Definition: fasp.h:583
REAL strong_coupled
strong coupled threshold for aggregate
Definition: fasp.h:670
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
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
dBSRmat A
pointer to the matrix at level level_num
Definition: fasp_block.h:207
#define PAIRWISE
Definition of aggregation types.
Definition: fasp_const.h:169
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
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 nb
dimension of each sub-block
Definition: fasp_block.h:56
#define SOLVER_SUPERLU
Definition: fasp_const.h:123
Mumps_data mumps
data for MUMPS
Definition: fasp_block.h:283
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
void print_amgcomplexity_bsr(AMG_data_bsr *mgl, const SHORT prtlvl)
Print complexities of AMG method for BSR matrices.
Definition: message.c:122
REAL ILU_relax
add the sum of dropped elements to diagonal element in proportion relax
Definition: fasp.h:389
INT Schwarz_maxlvl
maximal level for constructing the blocks
Definition: fasp.h:443
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:44
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
dvector w
temporary work space
Definition: fasp_block.h:280
#define INT
Definition: fasp.h:64
void fasp_blas_dcsr_rap_agg(dCSRmat *R, dCSRmat *A, dCSRmat *P, dCSRmat *RAP)
Triple sparse matrix multiplication B=R*A*P.
Definition: blas_csr.c:1148
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
INT near_kernel_dim
dimension of the near kernel for SAMG
Definition: fasp.h:767
dCSRmat fasp_dbsr_Linfinity_dcsr(dBSRmat *A)
get dCSRmat from a dBSRmat matrix using L_infinity norm of each small block
Definition: sparse_block.c:312
REAL ** near_kernel_basis
basis of near kernel space for SAMG
Definition: fasp.h:770
Data for multigrid levels. (BSR format)
Definition: fasp_block.h:198
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
SHORT fasp_ilu_dbsr_setup(dBSRmat *A, ILU_data *iludata, ILU_param *iluparam)
Get ILU decoposition of a BSR matrix A.
Definition: ilu_setup_bsr.c:45
SHORT fasp_amg_setup_sa_bsr(AMG_data_bsr *mgl, AMG_param *param)
Set up phase of smoothed aggregation AMG (BSR format)
Definition: amg_setup_sa.c:85
INT NNZ
number of nonzero sub-blocks in matrix A, NNZ
Definition: fasp_block.h:53
SHORT aggregation_type
aggregation type
Definition: fasp.h:646
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
#define MIN_CRATE
Definition: fasp_const.h:243
SHORT fasp_amg_setup_sa(AMG_data *mgl, AMG_param *param)
Set up phase of smoothed aggregation AMG.
Definition: amg_setup_sa.c:48
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
dvector b
pointer to the right-hand side at level level_num
Definition: fasp_block.h:216
SHORT coarse_solver
coarse solver type
Definition: fasp.h:628
#define SOLVER_PARDISO
Definition: fasp_const.h:126
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.
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
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
dCSRmat * P_nk
Prolongation for near kernal.
Definition: fasp_block.h:273
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
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
dvector diaginv
pointer to the diagonal inverse at level level_num
Definition: fasp_block.h:222
dCSRmat fasp_format_dbsr_dcsr(dBSRmat *B)
Transfer a 'dBSRmat' type matrix into a dCSRmat.
Definition: formats.c:495
SHORT cycle_type
type of AMG cycle
Definition: fasp.h:604
dCSRmat A
pointer to the matrix
Definition: fasp.h:510
INT COL
number of cols of sub-blocks in matrix A, N
Definition: fasp_block.h:50
dCSRmat Ac
pointer to the matrix at level level_num (csr format)
Definition: fasp_block.h:225
dCSRmat * R_nk
Resriction for near kernal.
Definition: fasp_block.h:276
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
INT ILU_levels
number of levels use ILU smoother
Definition: fasp_block.h:255
#define MIN(a, b)
Definition: fasp.h:73
#define AMLI_CYCLE
Definition: fasp_const.h:177
Mumps_data mumps
data for MUMPS
Definition: fasp.h:784
#define VMB
Definition: fasp_const.h:170
INT Schwarz_mmsize
maximal block size
Definition: fasp.h:703
INT pair_number
number of pairwise matchings
Definition: fasp.h:667
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
void fasp_blas_dcsr_mxm(dCSRmat *A, dCSRmat *B, dCSRmat *C)
Sparse matrix multiplication C=A*B.
Definition: blas_csr.c:759
INT num_levels
number of levels in use <= max_levels
Definition: fasp_block.h:204
dvector fasp_dbsr_getdiaginv(dBSRmat *A)
Get D^{-1} of matrix A.
Definition: sparse_bsr.c:392
Data for AMG solvers.
Definition: fasp.h:722
void fasp_blas_dbsr_rap(dBSRmat *R, dBSRmat *A, dBSRmat *P, dBSRmat *B)
dBSRmat sparse matrix multiplication B=R*A*P
Definition: blas_bsr.c:4895
#define SOLVER_UMFPACK
Definition: fasp_const.h:124
void * Numeric
pointer to the numerical factorization from UMFPACK
Definition: fasp.h:752
INT ROW
number of rows of sub-blocks in matrix A, M
Definition: fasp_block.h:47
INT Schwarz_maxlvl
maximal levels
Definition: fasp.h:706
INT fasp_dbsr_trans(dBSRmat *A, dBSRmat *AT)
Find A^T from given dBSRmat matrix A.
Definition: sparse_bsr.c:208
#define NL_AMLI_CYCLE
Definition: fasp_const.h:178
dCSRmat PP
pointer to the pressure block (only for reservoir simulation)
Definition: fasp_block.h:234
REAL ILU_droptol
drop tolerance for ILUt
Definition: fasp.h:691
dCSRmat * A_nk
Matrix data for near kernal.
Definition: fasp_block.h:270
INT job
work for MUMPS
Definition: fasp.h:467
INT Schwarz_type
type of Schwarz method
Definition: fasp.h:709
#define MAX_CRATE
Definition: fasp_const.h:244
INT ILU_lfil
level of fill-in for ILUs and ILUk
Definition: fasp.h:688
#define MIN_CDOF
Definition: fasp_const.h:242
void fasp_dbsr_free(dBSRmat *A)
Free memory space for BSR format sparse matrix.
Definition: sparse_bsr.c:133
dvector w
Temporary work space.
Definition: fasp.h:781
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
void * Numeric
pointer to the numerical dactorization from UMFPACK
Definition: fasp_block.h:228
#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