Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
amg_setup_ua.c
Go to the documentation of this file.
1 
10 #include <math.h>
11 #include <time.h>
12 
13 #include "fasp.h"
14 #include "fasp_functs.h"
15 #include "aggregation_csr.inl"
16 #include "aggregation_bsr.inl"
17 
18 static SHORT amg_setup_unsmoothP_unsmoothR(AMG_data *, AMG_param *);
19 static SHORT amg_setup_unsmoothP_unsmoothR_bsr(AMG_data_bsr *, AMG_param *);
20 
21 /*---------------------------------*/
22 /*-- Public Functions --*/
23 /*---------------------------------*/
24 
39  AMG_param *param)
40 {
41 #if DEBUG_MODE > 0
42  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
43  printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n",
44  mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
45 #endif
46 
47  SHORT status = amg_setup_unsmoothP_unsmoothR(mgl, param);
48 
49 #if DEBUG_MODE > 0
50  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
51 #endif
52 
53  return status;
54 }
55 
70  AMG_param *param)
71 {
72 #if DEBUG_MODE > 0
73  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
74 #endif
75 
76  SHORT status = amg_setup_unsmoothP_unsmoothR_bsr(mgl, param);
77 
78 #if DEBUG_MODE > 0
79  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
80 #endif
81 
82  return status;
83 }
84 
85 /*---------------------------------*/
86 /*-- Private Functions --*/
87 /*---------------------------------*/
88 
109 static SHORT amg_setup_unsmoothP_unsmoothR (AMG_data *mgl,
110  AMG_param *param)
111 {
112  const SHORT prtlvl = param->print_level;
113  const SHORT cycle_type = param->cycle_type;
114  const SHORT csolver = param->coarse_solver;
115  const SHORT min_cdof = MAX(param->coarse_dof,50);
116  const INT m = mgl[0].A.row;
117 
118  // empiric value
119  const REAL cplxmax = 3.0;
120  const REAL xsi = 0.6;
121  INT icum = 1.0;
122  REAL eta, fracratio;
123 
124  // local variables
125  SHORT max_levels = param->max_levels, lvl = 0, status = FASP_SUCCESS;
126  INT i;
127  REAL setup_start, setup_end;
128  ILU_param iluparam;
129  Schwarz_param swzparam;
130 
131 #if DEBUG_MODE > 0
132  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
133  printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n",
134  mgl[0].A.row, mgl[0].A.col, mgl[0].A.nnz);
135 #endif
136 
137  fasp_gettime(&setup_start);
138 
139  // level info (fine: 0; coarse: 1)
140  ivector *vertices = (ivector *)fasp_mem_calloc(max_levels,sizeof(ivector));
141 
142  // each level stores the information of the number of aggregations
143  INT *num_aggs = (INT *)fasp_mem_calloc(max_levels,sizeof(INT));
144 
145  // each level stores the information of the strongly coupled neighborhoods
146  dCSRmat *Neighbor = (dCSRmat *)fasp_mem_calloc(max_levels,sizeof(dCSRmat));
147 
148  // Initialize level information
149  for ( i = 0; i < max_levels; ++i ) num_aggs[i] = 0;
150 
151  mgl[0].near_kernel_dim = 1;
152  mgl[0].near_kernel_basis = (REAL **)fasp_mem_calloc(mgl->near_kernel_dim,sizeof(REAL*));
153 
154  for ( i = 0; i < mgl->near_kernel_dim; ++i ) {
155  mgl[0].near_kernel_basis[i] = (REAL *)fasp_mem_calloc(m,sizeof(REAL));
156  fasp_array_set(m, mgl[0].near_kernel_basis[i], 1.0);
157  }
158 
159  // Initialize ILU parameters
160  mgl->ILU_levels = param->ILU_levels;
161  if ( param->ILU_levels > 0 ) {
162  iluparam.print_level = param->print_level;
163  iluparam.ILU_lfil = param->ILU_lfil;
164  iluparam.ILU_droptol = param->ILU_droptol;
165  iluparam.ILU_relax = param->ILU_relax;
166  iluparam.ILU_type = param->ILU_type;
167  }
168 
169  // Initialize Schwarz parameters
170  mgl->Schwarz_levels = param->Schwarz_levels;
171  if ( param->Schwarz_levels > 0 ) {
172  swzparam.Schwarz_mmsize = param->Schwarz_mmsize;
173  swzparam.Schwarz_maxlvl = param->Schwarz_maxlvl;
174  swzparam.Schwarz_type = param->Schwarz_type;
175  swzparam.Schwarz_blksolver = param->Schwarz_blksolver;
176  }
177 
178  // Initialize AMLI coefficients
179  if ( cycle_type == AMLI_CYCLE ) {
180  const INT amlideg = param->amli_degree;
181  param->amli_coef = (REAL *)fasp_mem_calloc(amlideg+1,sizeof(REAL));
182  REAL lambda_max = 2.0, lambda_min = lambda_max/4;
183  fasp_amg_amli_coef(lambda_max, lambda_min, amlideg, param->amli_coef);
184  }
185 
186 #if DIAGONAL_PREF
187  fasp_dcsr_diagpref(&mgl[0].A); // reorder each row to make diagonal appear first
188 #endif
189 
190  /*----------------------------*/
191  /*--- checking aggregation ---*/
192  /*----------------------------*/
193 
194  // Pairwise matching algorithm requires diagonal preference ordering
195  if ( param->aggregation_type == PAIRWISE ) {
196  param->pair_number = MIN(param->pair_number, max_levels);
197  fasp_dcsr_diagpref(&mgl[0].A);
198  }
199 
200  // Main AMG setup loop
201  while ( (mgl[lvl].A.row > min_cdof) && (lvl < max_levels-1) ) {
202 
203 #if DEBUG_MODE > 2
204  printf("### DEBUG: level = %d, row = %d, nnz = %d\n",
205  lvl, mgl[lvl].A.row, mgl[lvl].A.nnz);
206 #endif
207 
208  /*-- Setup ILU decomposition if necessary */
209  if ( lvl < param->ILU_levels ) {
210  status = fasp_ilu_dcsr_setup(&mgl[lvl].A, &mgl[lvl].LU, &iluparam);
211  if ( status < 0 ) {
212  if ( prtlvl > PRINT_MIN ) {
213  printf("### WARNING: ILU setup on level-%d failed!\n", lvl);
214  printf("### WARNING: Disable ILU for level >= %d.\n", lvl);
215  }
216  param->ILU_levels = lvl;
217  }
218  }
219 
220  /*-- Setup Schwarz smoother if necessary */
221  if ( lvl < param->Schwarz_levels ) {
222  mgl[lvl].Schwarz.A=fasp_dcsr_sympat(&mgl[lvl].A);
223  fasp_dcsr_shift(&(mgl[lvl].Schwarz.A), 1);
224  fasp_Schwarz_setup(&mgl[lvl].Schwarz, &swzparam);
225  }
226 
227  /*-- Aggregation --*/
228  switch ( param->aggregation_type ) {
229 
230  case VMB: // VMB aggregation
231 
232  status = aggregation_vmb(&mgl[lvl].A, &vertices[lvl], param,
233  lvl+1, &Neighbor[lvl], &num_aggs[lvl]);
234 
235  /*-- Choose strength threshold adaptively --*/
236  if ( num_aggs[lvl]*4 > mgl[lvl].A.row )
237  param->strong_coupled /= 2;
238  else if ( num_aggs[lvl]*1.25 < mgl[lvl].A.row )
239  param->strong_coupled *= 2;
240 
241  break;
242 
243  default: // pairwise matching aggregation
244 
245  status = aggregation_pairwise(mgl, param, lvl, vertices,
246  &num_aggs[lvl]);
247 
248  break;
249  }
250 
251  // Check 1: Did coarsening step succeed?
252  if ( status < 0 ) {
253  // When error happens, stop at the current multigrid level!
254  if ( prtlvl > PRINT_MIN ) {
255  printf("### WARNING: Forming aggregates on level-%d failed!\n", lvl);
256  }
257  status = FASP_SUCCESS; break;
258  }
259 
260  /*-- Form Prolongation --*/
261  form_tentative_p(&vertices[lvl], &mgl[lvl].P, mgl[0].near_kernel_basis,
262  lvl+1, num_aggs[lvl]);
263 
264  // Check 2: Is coarse sparse too small?
265  if ( mgl[lvl].P.col < MIN_CDOF ) break;
266 
267  // Check 3: Does this coarsening step too aggressive?
268  if ( mgl[lvl].P.row > mgl[lvl].P.col * MAX_CRATE ) {
269  if ( prtlvl > PRINT_MIN ) {
270  printf("### WARNING: Coarsening might be too aggressive!\n");
271  printf("### WARNING: Fine level = %d, coarse level = %d. Discard!\n",
272  mgl[lvl].P.row, mgl[lvl].P.col);
273  }
274  break;
275  }
276 
277  // Check 4: Is this coarsening ratio too small?
278 #if 0
279  if ( (REAL)mgl[lvl].P.col > mgl[lvl].P.row * MIN_CRATE ) {
280  if ( prtlvl > PRINT_MIN ) {
281  printf("### WARNING: Coarsening rate is too small!\n");
282  printf("### WARNING: Fine level = %d, coarse level = %d. Discard!\n",
283  mgl[lvl].P.row, mgl[lvl].P.col);
284  }
285  break;
286  }
287 #endif
288  if ( (REAL)mgl[lvl].P.col > mgl[lvl].P.row * MIN_CRATE ) param->quality_bound *= 2.0;
289 
290  /*-- Form restriction --*/
291  fasp_dcsr_trans(&mgl[lvl].P, &mgl[lvl].R);
292 
293  /*-- Form coarse level stiffness matrix --*/
294  fasp_blas_dcsr_rap_agg(&mgl[lvl].R, &mgl[lvl].A, &mgl[lvl].P,
295  &mgl[lvl+1].A);
296 
297  fasp_dcsr_free(&Neighbor[lvl]);
298  fasp_ivec_free(&vertices[lvl]);
299 
300  ++lvl;
301 
302 #if DIAGONAL_PREF
303  fasp_dcsr_diagpref(&mgl[lvl].A); // reorder each row to make diagonal appear first
304 #endif
305 
306  }
307 
308  // Setup coarse level systems for direct solvers
309  switch (csolver) {
310 
311 #if WITH_MUMPS
312  case SOLVER_MUMPS: {
313  // Setup MUMPS direct solver on the coarsest level
314  mgl[lvl].mumps.job = 1;
315  fasp_solver_mumps_steps(&mgl[lvl].A, &mgl[lvl].b, &mgl[lvl].x, &mgl[lvl].mumps);
316  break;
317  }
318 #endif
319 
320 #if WITH_UMFPACK
321  case SOLVER_UMFPACK: {
322  // Need to sort the matrix A for UMFPACK to work
323  dCSRmat Ac_tran;
324  Ac_tran = fasp_dcsr_create(mgl[lvl].A.row, mgl[lvl].A.col, mgl[lvl].A.nnz);
325  fasp_dcsr_transz(&mgl[lvl].A, NULL, &Ac_tran);
326  // It is equivalent to do transpose and then sort
327  // fasp_dcsr_trans(&mgl[lvl].A, &Ac_tran);
328  // fasp_dcsr_sort(&Ac_tran);
329  fasp_dcsr_cp(&Ac_tran, &mgl[lvl].A);
330  fasp_dcsr_free(&Ac_tran);
331  mgl[lvl].Numeric = fasp_umfpack_factorize(&mgl[lvl].A, 0);
332  break;
333  }
334 #endif
335 
336 #if WITH_PARDISO
337  case SOLVER_PARDISO: {
338  fasp_dcsr_sort(&mgl[lvl].A);
339  fasp_pardiso_factorize(&mgl[lvl].A, &mgl[lvl].pdata, 0);
340  break;
341  }
342 #endif
343 
344  default:
345  // Do nothing!
346  break;
347  }
348 
349  // setup total level number and current level
350  mgl[0].num_levels = max_levels = lvl+1;
351  mgl[0].w = fasp_dvec_create(m);
352 
353  for ( lvl = 1; lvl < max_levels; ++lvl) {
354  INT mm = mgl[lvl].A.row;
355  mgl[lvl].num_levels = max_levels;
356  mgl[lvl].b = fasp_dvec_create(mm);
357  mgl[lvl].x = fasp_dvec_create(mm);
358 
359  mgl[lvl].cycle_type = cycle_type; // initialize cycle type!
360  mgl[lvl].ILU_levels = param->ILU_levels - lvl; // initialize ILU levels!
361  mgl[lvl].Schwarz_levels = param->Schwarz_levels -lvl; // initialize Schwarz!
362 
363  if ( cycle_type == NL_AMLI_CYCLE )
364  mgl[lvl].w = fasp_dvec_create(3*mm);
365  else
366  mgl[lvl].w = fasp_dvec_create(2*mm);
367  }
368 
369  // setup for cycle type of unsmoothed aggregation
370  eta = xsi/((1-xsi)*(cplxmax-1));
371  mgl[0].cycle_type = 1;
372  mgl[max_levels-1].cycle_type = 0;
373 
374  for (lvl = 1; lvl < max_levels-1; ++lvl) {
375  fracratio = (REAL)mgl[lvl].A.nnz/mgl[0].A.nnz;
376  mgl[lvl].cycle_type = (INT)(pow((REAL)xsi,(REAL)lvl)/(eta*fracratio*icum));
377  // safe-guard step: make cycle type >= 1 and <= 2
378  mgl[lvl].cycle_type = MAX(1, MIN(2, mgl[lvl].cycle_type));
379  icum = icum * mgl[lvl].cycle_type;
380  }
381 
382  if ( prtlvl > PRINT_NONE ) {
383  fasp_gettime(&setup_end);
384  print_amgcomplexity(mgl,prtlvl);
385  print_cputime("Unsmoothed aggregation setup", setup_end - setup_start);
386  }
387 
388  fasp_mem_free(Neighbor);
389  fasp_mem_free(vertices);
390  fasp_mem_free(num_aggs);
391 
392 #if DEBUG_MODE > 0
393  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
394 #endif
395 
396  return status;
397 }
398 
416 static SHORT amg_setup_unsmoothP_unsmoothR_bsr (AMG_data_bsr *mgl,
417  AMG_param *param)
418 {
419  const SHORT prtlvl = param->print_level;
420  const SHORT csolver = param->coarse_solver;
421  const SHORT min_cdof = MAX(param->coarse_dof,50);
422  const INT m = mgl[0].A.ROW;
423  const INT nb = mgl[0].A.nb;
424 
425  SHORT max_levels=param->max_levels;
426  SHORT i, lvl=0, status=FASP_SUCCESS;
427  REAL setup_start, setup_end;
428 
429  AMG_data *mgl_csr=fasp_amg_data_create(max_levels);
430 
431  dCSRmat temp1, temp2;
432 
433 #if DEBUG_MODE > 0
434  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
435  printf("### DEBUG: nr=%d, nc=%d, nnz=%d\n",
436  mgl[0].A.ROW, mgl[0].A.COL, mgl[0].A.NNZ);
437 #endif
438 
439  fasp_gettime(&setup_start);
440 
441  /*-----------------------*/
442  /*--local working array--*/
443  /*-----------------------*/
444 
445  // level info (fine: 0; coarse: 1)
446  ivector *vertices = (ivector *)fasp_mem_calloc(max_levels, sizeof(ivector));
447 
448  //each elvel stores the information of the number of aggregations
449  INT *num_aggs = (INT *)fasp_mem_calloc(max_levels, sizeof(INT));
450 
451  // each level stores the information of the strongly coupled neighborhoods
452  dCSRmat *Neighbor = (dCSRmat *)fasp_mem_calloc(max_levels, sizeof(dCSRmat));
453 
454  for ( i=0; i<max_levels; ++i ) num_aggs[i] = 0;
455 
456  /*-----------------------*/
457  /*-- setup null spaces --*/
458  /*-----------------------*/
459 
460  // null space for whole Jacobian
461  /*
462  mgl[0].near_kernel_dim = 1;
463  mgl[0].near_kernel_basis = (REAL **)fasp_mem_calloc(mgl->near_kernel_dim, sizeof(REAL*));
464 
465  for ( i=0; i < mgl->near_kernel_dim; ++i ) mgl[0].near_kernel_basis[i] = NULL;
466  */
467 
468  /*-----------------------*/
469  /*-- setup ILU param --*/
470  /*-----------------------*/
471 
472  // initialize ILU parameters
473  mgl->ILU_levels = param->ILU_levels;
474  ILU_param iluparam;
475 
476  if ( param->ILU_levels > 0 ) {
477  iluparam.print_level = param->print_level;
478  iluparam.ILU_lfil = param->ILU_lfil;
479  iluparam.ILU_droptol = param->ILU_droptol;
480  iluparam.ILU_relax = param->ILU_relax;
481  iluparam.ILU_type = param->ILU_type;
482  }
483 
484  /*----------------------------*/
485  /*--- checking aggregation ---*/
486  /*----------------------------*/
487  if (param->aggregation_type == PAIRWISE)
488  param->pair_number = MIN(param->pair_number, max_levels);
489 
490  // Main AMG setup loop
491  while ( (mgl[lvl].A.ROW > min_cdof) && (lvl < max_levels-1) ) {
492 
493 
494  /*-- setup ILU decomposition if necessary */
495  if ( lvl < param->ILU_levels ) {
496  status = fasp_ilu_dbsr_setup(&mgl[lvl].A, &mgl[lvl].LU, &iluparam);
497  if ( status < 0 ) {
498  if ( prtlvl > PRINT_MIN ) {
499  printf("### WARNING: ILU setup on level-%d failed!\n", lvl);
500  printf("### WARNING: Disable ILU for level >= %d.\n", lvl);
501  }
502  param->ILU_levels = lvl;
503  }
504  }
505 
506  /*-- get the diagonal inverse --*/
507  mgl[lvl].diaginv = fasp_dbsr_getdiaginv(&mgl[lvl].A);
508 
509  /*-- Aggregation --*/
510  //mgl[lvl].PP = fasp_dbsr_getblk_dcsr(&mgl[lvl].A);
511  mgl[lvl].PP = fasp_dbsr_Linfinity_dcsr(&mgl[lvl].A); // TODO: Try different way to form the scalar block!! -- Xiaozhe
512 
513  switch ( param->aggregation_type ) {
514 
515  case VMB: // VMB aggregation
516 
517  status = aggregation_vmb(&mgl[lvl].PP, &vertices[lvl], param,
518  lvl+1, &Neighbor[lvl], &num_aggs[lvl]);
519 
520  /*-- Choose strength threshold adaptively --*/
521  if ( num_aggs[lvl]*4 > mgl[lvl].PP.row )
522  param->strong_coupled /= 4;
523  else if ( num_aggs[lvl]*1.25 < mgl[lvl].PP.row )
524  param->strong_coupled *= 1.5;
525 
526  break;
527 
528  default: // pairwise matching aggregation
529 
530  //mgl_csr[lvl].A=fasp_dcsr_create(mgl[lvl].PP.row,mgl[lvl].PP.col,mgl[lvl].PP.nnz);
531  //fasp_dcsr_cp(&mgl[lvl].PP, &mgl_csr[lvl].A);
532  mgl_csr[lvl].A = mgl[lvl].PP;
533 
534  status = aggregation_pairwise(mgl_csr, param, lvl, vertices, &num_aggs[lvl]);
535 
536  // TODO: Need to design better algorithm for this part -- Xiaozhe
537 
538  break;
539  }
540 
541  if ( status < 0 ) {
542  // When error happens, force solver to use the current multigrid levels!
543  if ( prtlvl > PRINT_MIN ) {
544  printf("### WARNING: Forming aggregates on level-%d failed!\n", lvl);
545  }
546  status = FASP_SUCCESS; break;
547  }
548 
549  /* -- Form Prolongation --*/
550  //form_tentative_p_bsr(&vertices[lvl], &mgl[lvl].P, &mgl[0], lvl+1, num_aggs[lvl]);
551 
552  if ( lvl == 0 && mgl[0].near_kernel_dim >0 ) {
553  form_tentative_p_bsr1(&vertices[lvl], &mgl[lvl].P, &mgl[0], lvl+1,
554  num_aggs[lvl], mgl[0].near_kernel_dim,
555  mgl[0].near_kernel_basis);
556  }
557  else {
558 
559  form_boolean_p_bsr(&vertices[lvl], &mgl[lvl].P, &mgl[0], lvl+1,
560  num_aggs[lvl]);
561  }
562 
563  /*-- Form resitriction --*/
564  fasp_dbsr_trans(&mgl[lvl].P, &mgl[lvl].R);
565 
566  /*-- Form coarse level stiffness matrix --*/
567  fasp_blas_dbsr_rap(&mgl[lvl].R, &mgl[lvl].A, &mgl[lvl].P, &mgl[lvl+1].A);
568 
569  /* -- Form extra near kernal space if needed --*/
570  if (mgl[lvl].A_nk != NULL){
571 
572  mgl[lvl+1].A_nk = (dCSRmat *)fasp_mem_calloc(1, sizeof(dCSRmat));
573  mgl[lvl+1].P_nk = (dCSRmat *)fasp_mem_calloc(1, sizeof(dCSRmat));
574  mgl[lvl+1].R_nk = (dCSRmat *)fasp_mem_calloc(1, sizeof(dCSRmat));
575 
576  temp1 = fasp_format_dbsr_dcsr(&mgl[lvl].R);
577  fasp_blas_dcsr_mxm(&temp1, mgl[lvl].P_nk, mgl[lvl+1].P_nk);
578  fasp_dcsr_trans(mgl[lvl+1].P_nk, mgl[lvl+1].R_nk);
579  temp2 = fasp_format_dbsr_dcsr(&mgl[lvl+1].A);
580  fasp_blas_dcsr_rap(mgl[lvl+1].R_nk, &temp2, mgl[lvl+1].P_nk, mgl[lvl+1].A_nk);
581  fasp_dcsr_free(&temp1);
582  fasp_dcsr_free(&temp2);
583 
584  }
585 
586  fasp_dcsr_free(&Neighbor[lvl]);
587  fasp_ivec_free(&vertices[lvl]);
588 
589  ++lvl;
590  }
591 
592  // Setup coarse level systems for direct solvers (BSR version)
593  switch (csolver) {
594 
595 #if WITH_MUMPS
596  case SOLVER_MUMPS: {
597  // Setup MUMPS direct solver on the coarsest level
598  mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
599  mgl[lvl].mumps.job = 1;
600  fasp_solver_mumps_steps(&mgl[lvl].Ac, &mgl[lvl].b, &mgl[lvl].x, &mgl[lvl].mumps);
601  break;
602  }
603 #endif
604 
605 #if WITH_UMFPACK
606  case SOLVER_UMFPACK: {
607  // Need to sort the matrix A for UMFPACK to work
608  mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
609  dCSRmat Ac_tran;
610  fasp_dcsr_trans(&mgl[lvl].Ac, &Ac_tran);
611  fasp_dcsr_sort(&Ac_tran);
612  fasp_dcsr_cp(&Ac_tran, &mgl[lvl].Ac);
613  fasp_dcsr_free(&Ac_tran);
614  mgl[lvl].Numeric = fasp_umfpack_factorize(&mgl[lvl].Ac, 0);
615  break;
616  }
617 #endif
618 
619 #if WITH_SuperLU
620  case SOLVER_SUPERLU: {
621  /* Setup SuperLU direct solver on the coarsest level */
622  mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
623  }
624 #endif
625 
626 #if WITH_PARDISO
627  case SOLVER_PARDISO: {
628  mgl[lvl].Ac = fasp_format_dbsr_dcsr(&mgl[lvl].A);
629  fasp_dcsr_sort(&mgl[lvl].Ac);
630  fasp_pardiso_factorize(&mgl[lvl].Ac, &mgl[lvl].pdata, 0);
631  break;
632  }
633 #endif
634 
635  default:
636  // Do nothing!
637  break;
638  }
639 
640  // setup total level number and current level
641  mgl[0].num_levels = max_levels = lvl+1;
642  mgl[0].w = fasp_dvec_create(3*m*nb);
643 
644  if (mgl[0].A_nk != NULL){
645 
646 #if WITH_UMFPACK
647  // Need to sort the matrix A_nk for UMFPACK
648  fasp_dcsr_trans(mgl[0].A_nk, &temp1);
649  fasp_dcsr_sort(&temp1);
650  fasp_dcsr_cp(&temp1, mgl[0].A_nk);
651  fasp_dcsr_free(&temp1);
652 #endif
653 
654  }
655 
656  for ( lvl = 1; lvl < max_levels; lvl++ ) {
657  const INT mm = mgl[lvl].A.ROW*nb;
658  mgl[lvl].num_levels = max_levels;
659  mgl[lvl].b = fasp_dvec_create(mm);
660  mgl[lvl].x = fasp_dvec_create(mm);
661  mgl[lvl].w = fasp_dvec_create(3*mm);
662  mgl[lvl].ILU_levels = param->ILU_levels - lvl; // initialize ILU levels!
663 
664  if (mgl[lvl].A_nk != NULL){
665 
666 #if WITH_UMFPACK
667  // Need to sort the matrix A_nk for UMFPACK
668  fasp_dcsr_trans(mgl[lvl].A_nk, &temp1);
669  fasp_dcsr_sort(&temp1);
670  fasp_dcsr_cp(&temp1, mgl[lvl].A_nk);
671  fasp_dcsr_free(&temp1);
672 #endif
673 
674  }
675 
676  }
677 
678  if ( prtlvl > PRINT_NONE ) {
679  fasp_gettime(&setup_end);
680  print_amgcomplexity_bsr(mgl,prtlvl);
681  print_cputime("Unsmoothed aggregation (BSR) setup", setup_end - setup_start);
682  }
683 
684  fasp_mem_free(vertices);
685  fasp_mem_free(num_aggs);
686  fasp_mem_free(Neighbor);
687 
688 #if DEBUG_MODE > 0
689  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
690 #endif
691 
692  return status;
693 }
694 
695 /*---------------------------------*/
696 /*-- End of File --*/
697 /*---------------------------------*/
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
SHORT fasp_amg_setup_ua(AMG_data *mgl, AMG_param *param)
Set up phase of unsmoothed aggregation AMG.
Definition: amg_setup_ua.c:38
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
REAL quality_bound
quality threshold for pairwise aggregation
Definition: fasp.h:607
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
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
void fasp_array_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: array.c:48
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
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
AMG_data * fasp_amg_data_create(SHORT max_levels)
Create and initialize AMG_data for classical and SA AMG.
Definition: init.c:56
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
SHORT fasp_amg_setup_ua_bsr(AMG_data_bsr *mgl, AMG_param *param)
Set up phase of unsmoothed aggregation AMG (BSR format)
Definition: amg_setup_ua.c:69
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
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