Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
coarsening_rs.c
Go to the documentation of this file.
1 
12 #ifdef _OPENMP
13 #include <omp.h>
14 #endif
15 
16 #include "fasp.h"
17 #include "fasp_functs.h"
18 
19 // Private routines for RS coarsening
20 static INT cfsplitting_cls (dCSRmat *, iCSRmat *, ivector *);
21 static INT cfsplitting_clsp (dCSRmat *, iCSRmat *, ivector *);
22 static INT cfsplitting_agg (dCSRmat *, iCSRmat *, ivector *, INT);
23 static INT cfsplitting_mis (iCSRmat *, ivector *, ivector *);
24 static INT clean_ff_couplings (iCSRmat *, ivector *, INT, INT);
25 static INT compress_S (iCSRmat *);
26 static void strong_couplings (dCSRmat *, iCSRmat *, AMG_param *);
27 static void form_P_pattern_dir (dCSRmat *, iCSRmat *, ivector *, INT, INT);
28 static void form_P_pattern_std (dCSRmat *, iCSRmat *, ivector *, INT, INT);
29 static void ordering1 (iCSRmat *, ivector *);
30 
31 #include "linklist.inl"
32 
33 /*---------------------------------*/
34 /*-- Public Functions --*/
35 /*---------------------------------*/
36 
62  ivector *vertices,
63  dCSRmat *P,
64  iCSRmat *S,
65  AMG_param *param)
66 {
67  const SHORT coarse_type = param->coarsening_type;
68  const INT agg_path = param->aggressive_path;
69  const INT row = A->row;
70 
71  // local variables
72  SHORT interp_type = param->interpolation_type;
73  INT col = 0;
74 
75 #if DEBUG_MODE > 0
76  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
77 #endif
78 
79 #if DEBUG_MODE > 1
80  printf("### DEBUG: Step 1. Find strong connections ......\n");
81 #endif
82 
83  // make sure standard interp is used for aggressive coarsening
84  if ( coarse_type == COARSE_AC ) interp_type = INTERP_STD;
85 
86  // find strong couplings and return them in S
87  strong_couplings(A, S, param);
88 
89 #if DEBUG_MODE > 1
90  printf("### DEBUG: Step 2. C/F splitting ......\n");
91 #endif
92 
93  switch ( coarse_type ) {
94 
95  case COARSE_RS: // Classical coarsening
96  col = cfsplitting_cls(A, S, vertices); break;
97 
98  case COARSE_RSP: // Classical coarsening with positive connections
99  col = cfsplitting_clsp(A, S, vertices); break;
100 
101  case COARSE_AC: // Aggressive coarsening
102  col = cfsplitting_agg(A, S, vertices, agg_path); break;
103 
104  case COARSE_CR: // Compatible relaxation (Need to be fixed --Chensong)
105  col = fasp_amg_coarsening_cr(0, row-1, A, vertices, param); break;
106 
107 #if 0
108  case COARSE_MIS:
109  ivector order = fasp_ivec_create(row);
110  compress_S(S);
111  ordering1(S, &order);
112  col = cfsplitting_mis(S, vertices, &order);
113  fasp_ivec_free(&order);
114  break;
115 #endif
116 
117  default:
118  fasp_chkerr(ERROR_AMG_COARSE_TYPE, __FUNCTION__);
119 
120  }
121 
122  if ( col <= 0 ) return ERROR_UNKNOWN;
123 
124 #if DEBUG_MODE > 1
125  printf("### DEBUG: Step 3. Find support of C points ......\n");
126 #endif
127 
128  switch ( interp_type ) {
129 
130  case INTERP_DIR: // Direct interpolation or ...
131  case INTERP_ENG: // Energy-min interpolation
132  col = clean_ff_couplings(S, vertices, row, col);
133  form_P_pattern_dir(P, S, vertices, row, col);
134  break;
135 
136  case INTERP_STD: // Standard interpolation
137  form_P_pattern_std(P, S, vertices, row, col); break;
138 
139  default:
140  fasp_chkerr(ERROR_AMG_INTERP_TYPE, __FUNCTION__);
141 
142  }
143 
144 #if DEBUG_MODE > 0
145  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
146 #endif
147 
148  return FASP_SUCCESS;
149 }
150 
151 
152 /*---------------------------------*/
153 /*-- Private Functions --*/
154 /*---------------------------------*/
155 
173 static void strong_couplings (dCSRmat *A,
174  iCSRmat *S,
175  AMG_param *param )
176 {
177  const REAL max_row_sum = param->max_row_sum;
178  const REAL epsilon_str = param->strong_threshold;
179  const INT row = A->row, col = A->col, row1 = row+1;
180  const INT nnz = A->nnz;
181 
182  INT *ia = A->IA, *ja = A->JA;
183  REAL *aj = A->val;
184 
185  // local variables
186  INT i, j, begin_row, end_row;
187  REAL row_scl, row_sum;
188 
189  INT nthreads = 1, use_openmp = FALSE;
190 
191 #ifdef _OPENMP
192  if ( row > OPENMP_HOLDS ) {
193  use_openmp = TRUE;
194  nthreads = FASP_GET_NUM_THREADS();
195  }
196 #endif
197 
198  // get the diagonal entry of A: assume all connections are strong
199  dvector diag; fasp_dcsr_getdiag(0, A, &diag);
200 
201  // copy the structure of A to S
202  S->row = row; S->col = col; S->nnz = nnz; S->val = NULL;
203  S->IA = (INT *)fasp_mem_calloc(row1, sizeof(INT));
204  S->JA = (INT *)fasp_mem_calloc(nnz, sizeof(INT));
205  fasp_iarray_cp(row1, ia, S->IA);
206  fasp_iarray_cp(nnz, ja, S->JA);
207 
208  if ( use_openmp ) {
209 
210  INT mybegin, myend, myid;
211 #ifdef _OPENMP
212 #pragma omp parallel for private(myid, mybegin,myend,i,row_scl,row_sum,begin_row,end_row,j)
213 #endif
214  for ( myid = 0; myid < nthreads; myid++ ) {
215  FASP_GET_START_END(myid, nthreads, row, &mybegin, &myend);
216  for ( i = mybegin; i < myend; i++) {
217 
218  // Compute most negative entry in each row and row sum
219  row_scl = row_sum = 0.0;
220  begin_row = ia[i]; end_row = ia[i+1];
221  for ( j = begin_row; j < end_row; j++ ) {
222  row_scl = MIN(row_scl, aj[j]);
223  row_sum += aj[j];
224  }
225 
226  // Find diagonal entries of S and remove them later
227  for ( j = begin_row; j < end_row; j++ ) {
228  if ( ja[j] == i ) { S->JA[j] = -1; break; }
229  }
230 
231  // Mark entire row as weak couplings if strongly diagonal-dominant
232  if ( ABS(row_sum) > max_row_sum * ABS(diag.val[i]) ) {
233  for ( j = begin_row; j < end_row; j++ ) S->JA[j] = -1;
234  }
235  else {
236  for ( j = begin_row; j < end_row; j++) {
237  // If a_{ij} >= \epsilon_{str} * \min a_{ij}, the connection
238  // j->i is set to be weak; positive entries result in weak
239  // connections
240  if ( A->val[j] >= epsilon_str*row_scl ) S->JA[j] = -1;
241  }
242  }
243 
244  } // end for i
245  } // end for myid
246 
247  }
248 
249  else {
250 
251  for ( i = 0; i < row; ++i ) {
252 
253  // Compute row scale and row sum
254  row_scl = row_sum = 0.0;
255  begin_row = ia[i]; end_row = ia[i+1];
256 
257  for ( j = begin_row; j < end_row; j++ ) {
258 
259  // Originally: Not consider positive entries
260  // row_sum += aj[j];
261  // Now changed to --Chensong 05/17/2013
262  row_sum += ABS(aj[j]);
263 
264  // Originally: Not consider positive entries
265  // row_scl = MAX(row_scl, -aj[j]); // smallest negative
266  // Now changed to --Chensong 06/01/2013
267  if ( ja[j] != i ) row_scl = MAX(row_scl, ABS(aj[j])); // largest abs
268 
269  }
270 
271  // Multiply by the strength threshold
272  row_scl *= epsilon_str;
273 
274  // Find diagonal entries of S and remove them later
275  for ( j = begin_row; j < end_row; j++ ) {
276  if ( ja[j] == i ) { S->JA[j] = -1; break; }
277  }
278 
279  // Mark entire row as weak couplings if strongly diagonal-dominant
280  // Originally: Not consider positive entries
281  // if ( ABS(row_sum) > max_row_sum * ABS(diag.val[i]) ) {
282  // Now changed to --Chensong 05/17/2013
283  if ( row_sum < (2 - max_row_sum) * ABS(diag.val[i]) ) {
284  for ( j = begin_row; j < end_row; j++ ) S->JA[j] = -1;
285  }
286  else {
287  for ( j = begin_row; j < end_row; j++ ) {
288  if ( -A->val[j] <= row_scl ) S->JA[j] = -1; // only n-couplings
289  }
290  }
291  } // end for i
292 
293  } // end if openmp
294 
295  fasp_dvec_free(&diag);
296 }
297 
312 static INT compress_S (iCSRmat *S)
313 {
314  const INT row = S->row;
315  INT * ia = S->IA;
316 
317  // local variables
318  INT index, i, j, begin_row, end_row;
319 
320  // compress S: remove weak connections and form strong coupling matrix
321  for ( index = i = 0; i < row; ++i ) {
322 
323  begin_row = ia[i]; end_row = ia[i+1];
324 
325  ia[i] = index;
326  for ( j = begin_row; j < end_row; j++ ) {
327  if ( S->JA[j] > -1 ) S->JA[index++] = S->JA[j]; // strong couplings
328  }
329 
330  }
331 
332  S->nnz = S->IA[row] = index;
333 
334  if ( S->nnz <= 0 ) {
335  return ERROR_UNKNOWN;
336  }
337  else {
338  return FASP_SUCCESS;
339  }
340 }
341 
356 static void rem_positive_ff (dCSRmat *A,
357  iCSRmat *Stemp,
358  ivector *vertices)
359 {
360  const INT row = A->row;
361  INT *ia = A->IA, *vec = vertices->val;
362 
363  REAL row_scl, max_entry;
364  INT i, j, ji, max_index;
365 
366  for ( i = 0; i < row; ++i ) {
367 
368  if ( vec[i] != FGPT ) continue; // skip non F-variables
369 
370  row_scl = 0.0;
371  for ( ji = ia[i]; ji < ia[i+1]; ++ji ) {
372  j = A->JA[ji];
373  if ( j == i ) continue; // skip diagonal
374  row_scl = MAX(row_scl, ABS(A->val[ji])); // max abs entry
375  } // end for ji
376  row_scl *= 0.75;
377 
378  // looking for strong F-F connections
379  max_index = -1; max_entry = 0.0;
380  for ( ji = ia[i]; ji < ia[i+1]; ++ji ) {
381  j = A->JA[ji];
382  if ( j == i ) continue; // skip diagonal
383  if ( vec[j] != FGPT ) continue; // skip F-C connections
384  if ( A->val[ji] > row_scl ) {
385  Stemp->JA[ji] = j;
386  if ( A->val[ji] > max_entry ) {
387  max_entry = A->val[ji];
388  max_index = j; // max positive entry
389  }
390  }
391  } // end for ji
392 
393  // mark max positive entry as C-point
394  if ( max_index != -1 ) vec[max_index] = CGPT;
395 
396  } // end for i
397 
398 }
399 
422 static INT cfsplitting_cls (dCSRmat *A,
423  iCSRmat *S,
424  ivector *vertices)
425 {
426  const INT row = A->row;
427 
428  // local variables
429  INT col = 0;
430  INT maxmeas, maxnode, num_left = 0;
431  INT measure, newmeas;
432  INT i, j, k, l;
433  INT myid, mybegin, myend;
434  INT *vec = vertices->val;
435  INT *work = (INT*)fasp_mem_calloc(3*row,sizeof(INT));
436  INT *lists = work, *where = lists+row, *lambda = where+row;
437 
438 #if RS_C1
439  INT set_empty = 1;
440  INT jkeep = 0, cnt, index;
441  INT row_end_S, ji, row_end_S_nabor, jj;
442  INT *graph_array = lambda;
443 #else
444  INT *ia = A->IA;
445 #endif
446 
447  LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
448 
449  INT nthreads = 1, use_openmp = FALSE;
450 
451 #if DEBUG_MODE > 0
452  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
453 #endif
454 
455 #ifdef _OPENMP
456  if ( row > OPENMP_HOLDS ) {
457  use_openmp = TRUE;
458  nthreads = FASP_GET_NUM_THREADS();
459  }
460 #endif
461 
462  // 0. Compress S and form S_transpose
463  col = compress_S(S);
464  if ( col < 0 ) goto FINISHED; // compression failed!!!
465 
466  iCSRmat ST; fasp_icsr_trans(S, &ST);
467 
468  // 1. Initialize lambda
469  if ( use_openmp ) {
470 #ifdef _OPENMP
471 #pragma omp parallel for private(myid, mybegin,myend,i)
472 #endif
473  for ( myid = 0; myid < nthreads; myid++ ) {
474  FASP_GET_START_END(myid, nthreads, row, &mybegin, &myend);
475  for ( i = mybegin; i < myend; i++ ) lambda[i] = ST.IA[i+1] - ST.IA[i];
476  }
477  }
478  else {
479  for ( i = 0; i < row; ++i ) lambda[i] = ST.IA[i+1] - ST.IA[i];
480  }
481 
482  // 2. Before C/F splitting algorithm starts, filter out the variables which
483  // have no connections at all and mark them as special F-variables.
484  if ( use_openmp ) {
485 
486 #ifdef _OPENMP
487 #pragma omp parallel for reduction(+:num_left) private(myid, mybegin, myend, i)
488 #endif
489  for ( myid = 0; myid < nthreads; myid++ ) {
490  FASP_GET_START_END(myid, nthreads, row, &mybegin, &myend);
491  for ( i = mybegin; i < myend; i++ ) {
492 #if RS_C1
493  if ( S->IA[i+1] == S->IA[i] ) {
494 #else
495  if ( (ia[i+1]-ia[i]) <= 1 ) {
496 #endif
497  vec[i] = ISPT; // set i as an ISOLATED fine node
498  lambda[i] = 0;
499  }
500  else {
501  vec[i] = UNPT; // set i as a undecided node
502  num_left++;
503  }
504  }
505  } // end for myid
506 
507  }
508 
509  else {
510 
511  for ( i = 0; i < row; ++i ) {
512 
513 #if RS_C1
514  if ( S->IA[i+1] == S->IA[i] ) {
515 #else
516  if ( (ia[i+1]-ia[i]) <= 1 ) {
517 #endif
518  vec[i] = ISPT; // set i as an ISOLATED fine node
519  lambda[i] = 0;
520  }
521  else {
522  vec[i] = UNPT; // set i as a undecided node
523  num_left++;
524  }
525  } // end for i
526 
527  }
528 
529  // 3. Form linked list for lambda (max to min)
530  for ( i = 0; i < row; ++i ) {
531 
532  if ( vec[i] == ISPT ) continue; // skip isolated variables
533 
534  measure = lambda[i];
535 
536  if ( measure > 0 ) {
537  enter_list(&LoL_head, &LoL_tail, lambda[i], i, lists, where);
538  }
539  else {
540 
541  if ( measure < 0 ) printf("### WARNING: Negative lambda[%d]!\n", i);
542 
543  // Set variables with non-positive measure as F-variables
544  vec[i] = FGPT; // no strong connections, set i as fine node
545  --num_left;
546 
547  // Update lambda and linked list after i->F
548  for ( k = S->IA[i]; k < S->IA[i+1]; ++k ) {
549  j = S->JA[k];
550  if ( vec[j] == ISPT ) continue; // skip isolate variables
551  if ( j < i ) {
552  newmeas = lambda[j];
553  if ( newmeas > 0 ) {
554  remove_node(&LoL_head, &LoL_tail, newmeas, j, lists, where);
555  }
556  newmeas = ++(lambda[j]);
557  enter_list(&LoL_head, &LoL_tail, newmeas, j, lists, where);
558  }
559  else {
560  newmeas = ++(lambda[j]);
561  }
562  }
563 
564  } // end if measure
565 
566  } // end for i
567 
568  // 4. Main loop
569  while ( num_left > 0 ) {
570 
571  // pick $i\in U$ with $\max\lambda_i: C:=C\cup\{i\}, U:=U\\{i\}$
572  maxnode = LoL_head->head;
573  maxmeas = lambda[maxnode];
574  if ( maxmeas == 0 ) printf("### WARNING: Head of the list has measure 0!\n");
575 
576  vec[maxnode] = CGPT; // set maxnode as coarse node
577  lambda[maxnode] = 0;
578  --num_left;
579  remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
580  col++;
581 
582  // for all $j\in S_i^T\cap U: F:=F\cup\{j\}, U:=U\backslash\{j\}$
583  for ( i = ST.IA[maxnode]; i < ST.IA[maxnode+1]; ++i ) {
584 
585  j = ST.JA[i];
586 
587  if ( vec[j] != UNPT ) continue; // skip decided variables
588 
589  vec[j] = FGPT; // set j as fine node
590  remove_node(&LoL_head, &LoL_tail, lambda[j], j, lists, where);
591  --num_left;
592 
593  // Update lambda and linked list after j->F
594  for ( l = S->IA[j]; l < S->IA[j+1]; l++ ) {
595  k = S->JA[l];
596  if ( vec[k] == UNPT ) { // k is unknown
597  remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
598  newmeas = ++(lambda[k]);
599  enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
600  }
601  }
602 
603  } // end for i
604 
605  // Update lambda and linked list after maxnode->C
606  for ( i = S->IA[maxnode]; i < S->IA[maxnode+1]; ++i ) {
607 
608  j = S->JA[i];
609 
610  if ( vec[j] != UNPT ) continue; // skip decided variables
611 
612  measure = lambda[j];
613  remove_node(&LoL_head, &LoL_tail, measure, j, lists, where);
614  lambda[j] = --measure;
615 
616  if ( measure > 0 ) {
617  enter_list(&LoL_head, &LoL_tail, measure, j, lists, where);
618  }
619  else { // j is the only point left, set as fine variable
620  vec[j] = FGPT;
621  --num_left;
622 
623  // Update lambda and linked list after j->F
624  for ( l = S->IA[j]; l < S->IA[j+1]; l++ ) {
625  k = S->JA[l];
626  if ( vec[k] == UNPT ) { // k is unknown
627  remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
628  newmeas = ++(lambda[k]);
629  enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
630  }
631  } // end for l
632  } // end if
633 
634  } // end for
635 
636  } // end while
637 
638 #if RS_C1
639 
640  // C/F splitting of RS coarsening check C1 Criterion
641  fasp_iarray_set(row, graph_array, -1);
642  for (i = 0; i < row; i ++)
643  {
644  if (vec[i] == FGPT)
645  {
646  row_end_S = S->IA[i+1];
647  for (ji = S->IA[i]; ji < row_end_S; ji ++)
648  {
649  j = S->JA[ji];
650  if (vec[j] == CGPT)
651  {
652  graph_array[j] = i;
653  }
654  }
655  cnt = 0;
656  for (ji = S->IA[i]; ji < row_end_S; ji ++)
657  {
658  j = S->JA[ji];
659  if (vec[j] == FGPT)
660  {
661  set_empty = 1;
662  row_end_S_nabor = S->IA[j+1];
663  for (jj = S->IA[j]; jj < row_end_S_nabor; jj ++)
664  {
665  index = S->JA[jj];
666  if (graph_array[index] == i)
667  {
668  set_empty = 0;
669  break;
670  }
671  }
672  if (set_empty)
673  {
674  if (cnt == 0)
675  {
676  vec[j] = CGPT;
677  col++;
678  graph_array[j] = i;
679  jkeep = j;
680  cnt = 1;
681  }
682  else
683  {
684  vec[i] = CGPT;
685  vec[jkeep] = FGPT;
686  break;
687  }
688  }
689  }
690  }
691  }
692  }
693 
694 #endif
695 
696  fasp_icsr_free(&ST);
697 
698  if ( LoL_head ) {
699  list_ptr = LoL_head;
700  LoL_head->prev_node = NULL;
701  LoL_head->next_node = NULL;
702  LoL_head = list_ptr->next_node;
703  fasp_mem_free(list_ptr);
704  }
705 
706 FINISHED:
707  fasp_mem_free(work);
708 
709 #if DEBUG_MODE > 0
710  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
711 #endif
712 
713  return col;
714 }
715 
735 static INT cfsplitting_clsp (dCSRmat *A,
736  iCSRmat *S,
737  ivector *vertices)
738 {
739  const INT row = A->row;
740 
741  // local variables
742  INT col = 0;
743  INT maxmeas, maxnode, num_left = 0;
744  INT measure, newmeas;
745  INT i, j, k, l;
746  INT myid, mybegin, myend;
747 
748  INT *ia = A->IA, *vec = vertices->val;
749  INT *work = (INT*)fasp_mem_calloc(3*row,sizeof(INT));
750  INT *lists = work, *where = lists+row, *lambda = where+row;
751 
752  LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
753 
754  INT nthreads = 1, use_openmp = FALSE;
755 
756 #if DEBUG_MODE > 0
757  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
758 #endif
759 
760 #ifdef _OPENMP
761  if ( row > OPENMP_HOLDS ) {
762  use_openmp = TRUE;
763  nthreads = FASP_GET_NUM_THREADS();
764  }
765 #endif
766 
767  // 0. Compress S and form S_transpose
768  iCSRmat Stemp;
769  Stemp.row = S->row; Stemp.col = S->col; Stemp.nnz = S->nnz;
770  Stemp.IA = (INT *)fasp_mem_calloc(S->row+1, sizeof(INT));
771  Stemp.JA = (INT *)fasp_mem_calloc(S->nnz, sizeof(INT));
772  fasp_iarray_cp (S->row+1, S->IA, Stemp.IA);
773  fasp_iarray_cp (S->nnz, S->JA, Stemp.JA);
774 
775  if ( compress_S(S) < 0 ) goto FINISHED; // compression failed!!!
776 
777  iCSRmat ST; fasp_icsr_trans(S, &ST);
778 
779  // 1. Initialize lambda
780  if ( use_openmp ) {
781 #ifdef _OPENMP
782 #pragma omp parallel for private(myid, mybegin,myend,i)
783 #endif
784  for ( myid = 0; myid < nthreads; myid++ ) {
785  FASP_GET_START_END(myid, nthreads, row, &mybegin, &myend);
786  for ( i = mybegin; i < myend; i++ ) lambda[i] = ST.IA[i+1] - ST.IA[i];
787  }
788  }
789  else {
790  for ( i = 0; i < row; ++i ) lambda[i] = ST.IA[i+1] - ST.IA[i];
791  }
792 
793  // 2. Before C/F splitting algorithm starts, filter out the variables which
794  // have no connections at all and mark them as special F-variables.
795  if ( use_openmp ) {
796 
797 #ifdef _OPENMP
798 #pragma omp parallel for reduction(+:num_left) private(myid, mybegin, myend, i)
799 #endif
800  for ( myid = 0; myid < nthreads; myid++ ) {
801  FASP_GET_START_END(myid, nthreads, row, &mybegin, &myend);
802  for ( i = mybegin; i < myend; i++ ) {
803  if ( (ia[i+1]-ia[i]) <= 1 ) {
804  vec[i] = ISPT; // set i as an ISOLATED fine node
805  lambda[i] = 0;
806  }
807  else {
808  vec[i] = UNPT; // set i as a undecided node
809  num_left++;
810  }
811  }
812  } // end for myid
813 
814  }
815  else {
816 
817  for ( i = 0; i < row; ++i ) {
818  if ( (ia[i+1]-ia[i]) <= 1 ) {
819  vec[i] = ISPT; // set i as an ISOLATED fine node
820  lambda[i] = 0;
821  }
822  else {
823  vec[i] = UNPT; // set i as a undecided node
824  num_left++;
825  }
826  } // end for i
827 
828  }
829 
830  // 3. Form linked list for lambda (max to min)
831  for ( i = 0; i < row; ++i ) {
832 
833  if ( vec[i] == ISPT ) continue; // skip isolated variables
834 
835  measure = lambda[i];
836 
837  if ( measure > 0 ) {
838  enter_list(&LoL_head, &LoL_tail, lambda[i], i, lists, where);
839  }
840  else {
841 
842  if ( measure < 0 ) printf("### WARNING: Negative lambda[%d]!\n", i);
843 
844  // Set variables with non-positive measure as F-variables
845  vec[i] = FGPT; // no strong connections, set i as fine node
846  --num_left;
847 
848  // Update lambda and linked list after i->F
849  for ( k = S->IA[i]; k < S->IA[i+1]; ++k ) {
850 
851  j = S->JA[k];
852  if ( vec[j] == ISPT ) continue; // skip isolate variables
853 
854  if ( j < i ) { // only look at the previous points!!
855  newmeas = lambda[j];
856  if ( newmeas > 0 ) {
857  remove_node(&LoL_head, &LoL_tail, newmeas, j, lists, where);
858  }
859  newmeas = ++(lambda[j]);
860  enter_list(&LoL_head, &LoL_tail, newmeas, j, lists, where);
861  }
862  else { // will be checked later on
863  newmeas = ++(lambda[j]);
864  } // end if
865 
866  } // end for k
867 
868  } // end if measure
869 
870  } // end for i
871 
872  // 4. Main loop
873  while ( num_left > 0 ) {
874 
875  // pick $i\in U$ with $\max\lambda_i: C:=C\cup\{i\}, U:=U\\{i\}$
876  maxnode = LoL_head->head;
877  maxmeas = lambda[maxnode];
878  if ( maxmeas == 0 ) printf("### WARNING: Head of the list has measure 0!\n");
879 
880  vec[maxnode] = CGPT; // set maxnode as coarse node
881  lambda[maxnode] = 0;
882  --num_left;
883  remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
884  col++;
885 
886  // for all $j\in S_i^T\cap U: F:=F\cup\{j\}, U:=U\backslash\{j\}$
887  for ( i = ST.IA[maxnode]; i < ST.IA[maxnode+1]; ++i ) {
888 
889  j = ST.JA[i];
890 
891  if ( vec[j] != UNPT ) continue; // skip decided variables
892 
893  vec[j] = FGPT; // set j as fine node
894  remove_node(&LoL_head, &LoL_tail, lambda[j], j, lists, where);
895  --num_left;
896 
897  // Update lambda and linked list after j->F
898  for ( l = S->IA[j]; l < S->IA[j+1]; l++ ) {
899  k = S->JA[l];
900  if ( vec[k] == UNPT ) { // k is unknown
901  remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
902  newmeas = ++(lambda[k]);
903  enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
904  }
905  }
906 
907  } // end for i
908 
909  // Update lambda and linked list after maxnode->C
910  for ( i = S->IA[maxnode]; i < S->IA[maxnode+1]; ++i ) {
911 
912  j = S->JA[i];
913 
914  if ( vec[j] != UNPT ) continue; // skip decided variables
915 
916  measure = lambda[j];
917  remove_node(&LoL_head, &LoL_tail, measure, j, lists, where);
918  lambda[j] = --measure;
919 
920  if ( measure > 0 ) {
921  enter_list(&LoL_head, &LoL_tail, measure, j, lists, where);
922  }
923  else { // j is the only point left, set as fine variable
924  vec[j] = FGPT;
925  --num_left;
926 
927  // Update lambda and linked list after j->F
928  for ( l = S->IA[j]; l < S->IA[j+1]; l++ ) {
929  k = S->JA[l];
930  if ( vec[k] == UNPT ) { // k is unknown
931  remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
932  newmeas = ++(lambda[k]);
933  enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
934  }
935  } // end for l
936  } // end if
937 
938  } // end for
939 
940  } // end while
941 
942  fasp_icsr_free(&ST);
943 
944  if ( LoL_head ) {
945  list_ptr = LoL_head;
946  LoL_head->prev_node = NULL;
947  LoL_head->next_node = NULL;
948  LoL_head = list_ptr->next_node;
949  fasp_mem_free(list_ptr);
950  }
951 
952  // Enforce F-C connections. Adding this step helps for the ExxonMobil test
953  // problems! Need more tests though --Chensong 06/08/2013
954  // col = clean_ff_couplings(S, vertices, row, col);
955 
956  rem_positive_ff(A, &Stemp, vertices);
957 
958  if ( compress_S(&Stemp) < 0 ) goto FINISHED; // compression failed!!!
959 
960  S->row = Stemp.row;
961  S->col = Stemp.col;
962  S->nnz = Stemp.nnz;
963 
964  fasp_mem_free(S->IA); S->IA = Stemp.IA;
965  fasp_mem_free(S->JA); S->JA = Stemp.JA;
966 
967 FINISHED:
968  fasp_mem_free(work);
969 
970 #if DEBUG_MODE > 0
971  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
972 #endif
973 
974  return col;
975 }
976 
997 static void strong_couplings_agg1 (dCSRmat *A,
998  iCSRmat *S,
999  iCSRmat *Sh,
1000  ivector *vertices,
1001  ivector *CGPT_index,
1002  ivector *CGPT_rindex)
1003 {
1004  const INT row = A->row;
1005 
1006  // local variables
1007  INT i, j, k;
1008  INT num_c, count, ci, cj, ck, fj, cck;
1009  INT *cp_index, *cp_rindex, *visited;
1010  INT *vec = vertices->val;
1011 
1012  // count the number of coarse grid points
1013  for ( num_c = i = 0; i < row; i++ ) {
1014  if ( vec[i] == CGPT ) num_c++;
1015  }
1016 
1017  // for the reverse indexing of coarse grid points
1018  fasp_ivec_alloc(row, CGPT_rindex);
1019  cp_rindex = CGPT_rindex->val;
1020 
1021  // generate coarse grid point index
1022  fasp_ivec_alloc(num_c, CGPT_index);
1023  cp_index = CGPT_index->val;
1024  for ( j = i = 0; i < row; i++ ) {
1025  if ( vec[i] == CGPT ) {
1026  cp_index[j] = i;
1027  cp_rindex[i] = j;
1028  j++;
1029  }
1030  }
1031 
1032  // allocate space for Sh
1033  Sh->row = Sh->col = num_c;
1034  Sh->val = Sh->JA = NULL;
1035  Sh->IA = (INT*)fasp_mem_calloc(Sh->row+1, sizeof(INT));
1036 
1037  // record the number of times some coarse point is visited
1038  visited = (INT*)fasp_mem_calloc(num_c, sizeof(INT));
1039  fasp_iarray_set(num_c, visited, -1);
1040 
1041  /**********************************************/
1042  /* step 1: Find first the structure IA of Sh */
1043  /**********************************************/
1044 
1045  Sh->IA[0] = 0;
1046 
1047  for ( ci = 0; ci < Sh->row; ci++ ) {
1048 
1049  i = cp_index[ci]; // find the index of the ci-th coarse grid point
1050 
1051  // number of coarse point that i is strongly connected to w.r.t. S(p,2)
1052  count = 0;
1053 
1054  // visit all the fine neighbors that ci is strongly connected to
1055  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1056 
1057  fj = S->JA[j];
1058 
1059  if ( vec[fj] == CGPT && fj != i ) {
1060  cj = cp_rindex[fj];
1061  if ( visited[cj] != ci ) {
1062  visited[cj] = ci; // mark as strongly connected from ci
1063  count++;
1064  }
1065 
1066  }
1067 
1068  else if ( vec[fj] == FGPT ) { // fine grid point,
1069 
1070  // find all the coarse neighbors that fj is strongly connected to
1071  for ( k = S->IA[fj]; k < S->IA[fj+1]; k++ ) {
1072  ck = S->JA[k];
1073  if ( vec[ck] == CGPT && ck != i ) { // it is a coarse grid point
1074  if ( cp_rindex[ck] >= num_c ) {
1075  printf("### ERROR: ck=%d, num_c=%d, out of bound!\n",
1076  ck, num_c);
1077  fasp_chkerr(ERROR_AMG_COARSEING, __FUNCTION__);
1078  }
1079  cck = cp_rindex[ck];
1080 
1081  if ( visited[cck] != ci ) {
1082  visited[cck] = ci; // mark as strongly connected from ci
1083  count++;
1084  }
1085  } //end if
1086  } //end for k
1087 
1088  } //end if
1089 
1090  } //end for j
1091 
1092  Sh->IA[ci+1] = Sh->IA[ci] + count;
1093 
1094  } //end for i
1095 
1096  /*************************/
1097  /* step 2: Find JA of Sh */
1098  /*************************/
1099 
1100  fasp_iarray_set(num_c, visited, -1); // reset visited
1101 
1102  Sh->nnz = Sh->IA[Sh->row];
1103  Sh->JA = (INT*)fasp_mem_calloc(Sh->nnz, sizeof(INT));
1104 
1105  for ( ci = 0; ci < Sh->row; ci++ ) {
1106 
1107  i = cp_index[ci]; // find the index of the i-th coarse grid point
1108  count = Sh->IA[ci]; // count for coarse points
1109 
1110  // visit all the fine neighbors that ci is strongly connected to
1111  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1112 
1113  fj = S->JA[j];
1114 
1115  if ( vec[fj] == CGPT && fj != i ) {
1116  cj = cp_rindex[fj];
1117  if ( visited[cj] != ci ) { // not visited yet
1118  visited[cj] = ci;
1119  Sh->JA[count] = cj;
1120  count++;
1121  }
1122  }
1123  else if ( vec[fj] == FGPT ) { // fine grid point,
1124  //find all the coarse neighbors that fj is strongly connected to
1125  for ( k = S->IA[fj]; k < S->IA[fj+1]; k++ ) {
1126  ck = S->JA[k];
1127  if ( vec[ck] == CGPT && ck != i ) { // coarse grid point
1128  cck = cp_rindex[ck];
1129  if ( visited[cck] != ci ) { // not visited yet
1130  visited[cck] = ci;
1131  Sh->JA[count] = cck;
1132  count++;
1133  }
1134  } // end if
1135  } // end for k
1136  } // end if
1137 
1138  } // end for j
1139 
1140  if ( count != Sh->IA[ci+1] ) {
1141  printf("### WARNING: Inconsistent numbers of nonzeros!\n ");
1142  }
1143 
1144  } // end for ci
1145 
1146  fasp_mem_free(visited);
1147 }
1148 
1175 static void strong_couplings_agg2 (dCSRmat *A,
1176  iCSRmat *S,
1177  iCSRmat *Sh,
1178  ivector *vertices,
1179  ivector *CGPT_index,
1180  ivector *CGPT_rindex)
1181 {
1182  const INT row = A->row;
1183 
1184  // local variables
1185  INT i, j, k;
1186  INT num_c, count, ci, cj, ck, fj, cck;
1187  INT *cp_index, *cp_rindex, *visited;
1188  INT *vec = vertices->val;
1189 
1190  // count the number of coarse grid points
1191  for ( num_c = i = 0; i < row; i++ ) {
1192  if ( vec[i] == CGPT ) num_c++;
1193  }
1194 
1195  // for the reverse indexing of coarse grid points
1196  fasp_ivec_alloc(row, CGPT_rindex);
1197  cp_rindex = CGPT_rindex->val;
1198 
1199  // generate coarse grid point index
1200  fasp_ivec_alloc(num_c, CGPT_index);
1201  cp_index = CGPT_index->val;
1202  for ( j = i = 0; i < row; i++ ) {
1203  if ( vec[i] == CGPT ) {
1204  cp_index[j] = i;
1205  cp_rindex[i] = j;
1206  j++;
1207  }
1208  }
1209 
1210  // allocate space for Sh
1211  Sh->row = Sh->col = num_c;
1212  Sh->val = Sh->JA = NULL;
1213  Sh->IA = (INT*)fasp_mem_calloc(Sh->row+1, sizeof(INT));
1214 
1215  // record the number of times some coarse point is visited
1216  visited = (INT*)fasp_mem_calloc(num_c, sizeof(INT));
1217  memset(visited, 0, sizeof(INT)*num_c);
1218 
1219  /**********************************************/
1220  /* step 1: Find first the structure IA of Sh */
1221  /**********************************************/
1222 
1223  Sh->IA[0] = 0;
1224 
1225  for ( ci = 0; ci < Sh->row; ci++ ) {
1226 
1227  i = cp_index[ci]; // find the index of the ci-th coarse grid point
1228 
1229  // number of coarse point that i is strongly connected to w.r.t. S(p,2)
1230  count = 0;
1231 
1232  // visit all the fine neighbors that ci is strongly connected to
1233  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1234 
1235  fj = S->JA[j];
1236 
1237  if ( vec[fj] == CGPT && fj != i ) {
1238  cj = cp_rindex[fj];
1239  if ( visited[cj] != ci+1 ) { // not visited yet
1240  visited[cj] = ci+1; // mark as strongly connected from ci
1241  count++;
1242  }
1243  }
1244 
1245  else if ( vec[fj] == FGPT ) { // fine grid point
1246 
1247  // find all the coarse neighbors that fj is strongly connected to
1248  for ( k = S->IA[fj]; k < S->IA[fj+1]; k++ ) {
1249 
1250  ck = S->JA[k];
1251 
1252  if ( vec[ck] == CGPT && ck != i ) { // coarse grid point
1253  if ( cp_rindex[ck] >= num_c ) {
1254  printf("### ERROR: ck=%d, num_c=%d, out of bound!\n",
1255  ck, num_c);
1256  fasp_chkerr(ERROR_AMG_COARSEING, __FUNCTION__);
1257  }
1258  cck = cp_rindex[ck];
1259 
1260  if ( visited[cck] == ci+1 ) {
1261  // visited already!
1262  }
1263  else if ( visited[cck] == -ci-1 ) {
1264  visited[cck] = ci+1; // mark as strongly connected from ci
1265  count++;
1266  }
1267  else {
1268  visited[cck] = -ci-1; // mark as visited
1269  }
1270 
1271  } //end if vec[ck]
1272 
1273  } // end for k
1274 
1275  } // end if vec[fj]
1276 
1277  } // end for j
1278 
1279  Sh->IA[ci+1] = Sh->IA[ci] + count;
1280 
1281  } //end for i
1282 
1283  /*************************/
1284  /* step 2: Find JA of Sh */
1285  /*************************/
1286 
1287  memset(visited, 0, sizeof(INT)*num_c); // reset visited
1288 
1289  Sh->nnz = Sh->IA[Sh->row];
1290  Sh->JA = (INT*)fasp_mem_calloc(Sh->nnz,sizeof(INT));
1291 
1292  for ( ci = 0; ci < Sh->row; ci++ ) {
1293 
1294  i = cp_index[ci]; // find the index of the i-th coarse grid point
1295  count = Sh->IA[ci]; // count for coarse points
1296 
1297  // visit all the fine neighbors that ci is strongly connected to
1298  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1299 
1300  fj = S->JA[j];
1301 
1302  if ( vec[fj] == CGPT && fj != i ) {
1303  cj = cp_rindex[fj];
1304  if ( visited[cj] != ci+1 ) { // not visited yet
1305  visited[cj] = ci+1;
1306  Sh->JA[count] = cj;
1307  count++;
1308  }
1309  }
1310 
1311  else if ( vec[fj] == FGPT ) { // fine grid point
1312 
1313  // find all the coarse neighbors that fj is strongly connected to
1314  for ( k = S->IA[fj]; k < S->IA[fj+1]; k++ ) {
1315 
1316  ck = S->JA[k];
1317 
1318  if ( vec[ck] == CGPT && ck != i ) { // coarse grid point
1319  cck = cp_rindex[ck];
1320  if ( visited[cck] == ci+1 ) {
1321  // visited before
1322  }
1323  else if ( visited[cck] == -ci-1 ) {
1324  visited[cck] = ci+1;
1325  Sh->JA[count] = cck;
1326  count++;
1327  }
1328  else {
1329  visited[cck] = -ci-1;
1330  }
1331  } // end if vec[ck]
1332 
1333  } // end for k
1334 
1335  } // end if vec[fj]
1336 
1337  } // end for j
1338 
1339  if ( count != Sh->IA[ci+1] ) {
1340  printf("### WARNING: Inconsistent numbers of nonzeros!\n ");
1341  }
1342 
1343  } // end for ci
1344 
1345  fasp_mem_free(visited);
1346 }
1347 
1370 static INT cfsplitting_agg (dCSRmat *A,
1371  iCSRmat *S,
1372  ivector *vertices,
1373  INT aggressive_path)
1374 {
1375  const INT row = A->row;
1376  INT col = 0; // initialize col(P): returning output
1377 
1378  // local variables
1379  INT *vec = vertices->val, *cp_index;
1380  INT maxmeas, maxnode, num_left = 0;
1381  INT measure, newmeas;
1382  INT i, j, k, l, m, ci, cj, ck, cl, num_c;
1383  SHORT IS_CNEIGH;
1384 
1385  INT *work = (INT*)fasp_mem_calloc(3*row,sizeof(INT));
1386  INT *lists = work, *where = lists+row, *lambda = where+row;
1387 
1388  ivector CGPT_index, CGPT_rindex;
1389  LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
1390 
1391  // Sh is for the strong coupling matrix between temporary CGPTs
1392  // ShT is the transpose of Sh
1393  // Snew is for combining the information from S and Sh
1394  iCSRmat ST, Sh, ShT;
1395 
1396 #if DEBUG_MODE > 0
1397  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
1398 #endif
1399 
1400 #ifdef _OPENMP
1401  // variables for OpenMP
1402  INT myid, mybegin, myend;
1403  INT nthreads = FASP_GET_NUM_THREADS();
1404 #endif
1405 
1406  /************************************************************/
1407  /* Coarsening Phase ONE: find temporary coarse level points */
1408  /************************************************************/
1409 
1410  num_c = cfsplitting_cls(A, S, vertices);
1411  fasp_icsr_trans(S, &ST);
1412 
1413  /************************************************************/
1414  /* Coarsening Phase TWO: find real coarse level points */
1415  /************************************************************/
1416 
1417  // find Sh, the strong coupling between coarse grid points S(path,2)
1418  if ( aggressive_path < 2 )
1419  strong_couplings_agg1(A, S, &Sh, vertices, &CGPT_index, &CGPT_rindex);
1420  else
1421  strong_couplings_agg2(A, S, &Sh, vertices, &CGPT_index, &CGPT_rindex);
1422 
1423  fasp_icsr_trans(&Sh, &ShT);
1424 
1425  CGPT_index.row = num_c;
1426  CGPT_rindex.row = row;
1427  cp_index = CGPT_index.val;
1428 
1429  // 1. Initialize lambda
1430 #ifdef _OPENMP
1431 #pragma omp parallel for if(num_c>OPENMP_HOLDS)
1432 #endif
1433  for ( ci = 0; ci < num_c; ++ci ) lambda[ci] = ShT.IA[ci+1]-ShT.IA[ci];
1434 
1435  // 2. Form linked list for lambda (max to min)
1436  for ( ci = 0; ci < num_c; ++ci ) {
1437 
1438  i = cp_index[ci];
1439  measure = lambda[ci];
1440 
1441  if ( vec[i] == ISPT ) continue; // skip isolated points
1442 
1443  if ( measure > 0 ) {
1444  enter_list(&LoL_head, &LoL_tail, lambda[ci], ci, lists, where);
1445  num_left++;
1446  }
1447  else {
1448  if ( measure < 0) printf("### WARNING: Negative lambda[%d]!\n", i);
1449 
1450  vec[i] = FGPT; // set i as fine node
1451 
1452  // update the lambda value in the CGPT neighbor of i
1453  for ( ck = Sh.IA[ci]; ck < Sh.IA[ci+1]; ++ck ) {
1454 
1455  cj = Sh.JA[ck];
1456  j = cp_index[cj];
1457 
1458  if ( vec[j] == ISPT ) continue;
1459 
1460  if ( cj < ci ) {
1461  newmeas = lambda[cj];
1462  if ( newmeas > 0 ) {
1463  remove_node(&LoL_head, &LoL_tail, newmeas, cj, lists, where);
1464  num_left--;
1465  }
1466  newmeas = ++(lambda[cj]);
1467  enter_list(&LoL_head, &LoL_tail, newmeas, cj, lists, where);
1468  num_left++;
1469  }
1470  else {
1471  newmeas = ++(lambda[cj]);
1472  } // end if cj<ci
1473 
1474  } // end for ck
1475 
1476  } // end if
1477 
1478  } // end for ci
1479 
1480  // 3. Main loop
1481  while ( num_left > 0 ) {
1482 
1483  // pick $i\in U$ with $\max\lambda_i: C:=C\cup\{i\}, U:=U\\{i\}$
1484  maxnode = LoL_head->head;
1485  maxmeas = lambda[maxnode];
1486  if ( maxmeas == 0 ) printf("### WARNING: Head of the list has measure 0!\n");
1487 
1488  // mark maxnode as real coarse node, labelled as number 3
1489  vec[cp_index[maxnode]] = 3;
1490  --num_left;
1491  remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
1492  lambda[maxnode] = 0;
1493  col++; // count for the real coarse node after aggressive coarsening
1494 
1495  // for all $j\in S_i^T\cap U: F:=F\cup\{j\}, U:=U\backslash\{j\}$
1496  for ( ci = ShT.IA[maxnode]; ci < ShT.IA[maxnode+1]; ++ci ) {
1497 
1498  cj = ShT.JA[ci];
1499  j = cp_index[cj];
1500 
1501  if ( vec[j] != CGPT ) continue; // skip if j is not C-point
1502 
1503  vec[j] = 4; // set j as 4--fake CGPT
1504  remove_node(&LoL_head, &LoL_tail, lambda[cj], cj, lists, where);
1505  --num_left;
1506 
1507  // update the measure for neighboring points
1508  for ( cl = Sh.IA[cj]; cl < Sh.IA[cj+1]; cl++ ) {
1509  ck = Sh.JA[cl];
1510  k = cp_index[ck];
1511  if ( vec[k] == CGPT ) { // k is temporary CGPT
1512  remove_node(&LoL_head, &LoL_tail, lambda[ck], ck, lists, where);
1513  newmeas = ++(lambda[ck]);
1514  enter_list(&LoL_head, &LoL_tail, newmeas, ck, lists, where);
1515  }
1516  }
1517 
1518  } // end for ci
1519 
1520  // Update lambda and linked list after maxnode->C
1521  for ( ci = Sh.IA[maxnode]; ci < Sh.IA[maxnode+1]; ++ci ) {
1522 
1523  cj = Sh.JA[ci];
1524  j = cp_index[cj];
1525 
1526  if ( vec[j] != CGPT ) continue; // skip if j is not C-point
1527 
1528  measure = lambda[cj];
1529  remove_node(&LoL_head, &LoL_tail, measure, cj, lists, where);
1530  lambda[cj] = --measure;
1531 
1532  if ( measure > 0 ) {
1533  enter_list(&LoL_head, &LoL_tail,measure, cj, lists, where);
1534  }
1535  else {
1536  vec[j] = 4; // set j as fake CGPT variable
1537  --num_left;
1538  for ( cl = Sh.IA[cj]; cl < Sh.IA[cj+1]; cl++ ) {
1539  ck = Sh.JA[cl];
1540  k = cp_index[ck];
1541  if ( vec[k] == CGPT ) { // k is temporary CGPT
1542  remove_node(&LoL_head, &LoL_tail, lambda[ck], ck, lists, where);
1543  newmeas = ++(lambda[ck]);
1544  enter_list(&LoL_head, &LoL_tail, newmeas, ck, lists, where);
1545  }
1546  } // end for l
1547  } // end if
1548 
1549  } // end for
1550 
1551  } // while
1552 
1553  // 4. reorganize the variable type: mark temporary CGPT--1 and fake CGPT--4 as
1554  // FGPT; mark real CGPT--3 to be CGPT
1555 #ifdef _OPENMP
1556 #pragma omp parallel for if(row>OPENMP_HOLDS)
1557 #endif
1558  for ( i = 0; i < row; i++ ) {
1559  if ( vec[i] == CGPT || vec[i] == 4 ) vec[i] = FGPT;
1560  }
1561 
1562 #ifdef _OPENMP
1563 #pragma omp parallel for if(row>OPENMP_HOLDS)
1564 #endif
1565  for ( i = 0; i < row; i++ ) {
1566  if ( vec[i] == 3 ) vec[i] = CGPT;
1567  }
1568 
1569  /************************************************************/
1570  /* Coarsening Phase THREE: all the FGPTs which have no CGPT */
1571  /* neighbors within distance 2. Change them into CGPT such */
1572  /* that the standard interpolation works! */
1573  /************************************************************/
1574 
1575  for ( i = 0; i < row; i++ ) {
1576 
1577  if ( vec[i] != FGPT ) continue;
1578 
1579  IS_CNEIGH = FALSE; // whether there exist CGPT neighbors within distance of 2
1580 
1581  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1582 
1583  if ( IS_CNEIGH ) break;
1584 
1585  k = S->JA[j];
1586 
1587  if ( vec[k] == CGPT ) {
1588  IS_CNEIGH = TRUE;
1589  }
1590  else if ( vec[k] == FGPT ) {
1591  for ( l = S->IA[k]; l < S->IA[k+1]; l++ ) {
1592  m = S->JA[l];
1593  if ( vec[m] == CGPT ) {
1594  IS_CNEIGH = TRUE; break;
1595  }
1596  } // end for l
1597  }
1598 
1599  } // end for j
1600 
1601  // no CGPT neighbors in distance <= 2, mark i as CGPT
1602  if ( !IS_CNEIGH ) {
1603  vec[i] = CGPT; col++;
1604  }
1605 
1606  } // end for i
1607 
1608  if ( LoL_head ) {
1609  list_ptr = LoL_head;
1610  LoL_head->prev_node = NULL;
1611  LoL_head->next_node = NULL;
1612  LoL_head = list_ptr->next_node;
1613  fasp_mem_free(list_ptr);
1614  }
1615 
1616  fasp_ivec_free(&CGPT_index);
1617  fasp_ivec_free(&CGPT_rindex);
1618  fasp_icsr_free(&Sh);
1619  fasp_icsr_free(&ST);
1620  fasp_icsr_free(&ShT);
1621  fasp_mem_free(work);
1622 
1623 #if DEBUG_MODE > 0
1624  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1625 #endif
1626 
1627  return col;
1628 }
1629 
1653 static INT clean_ff_couplings (iCSRmat *S,
1654  ivector *vertices,
1655  INT row,
1656  INT col)
1657 {
1658  // local variables
1659  INT *vec = vertices->val;
1660  INT *cindex = (INT *)fasp_mem_calloc(row, sizeof(INT));
1661  INT set_empty = TRUE, C_i_nonempty = FALSE;
1662  INT ci_tilde = -1, ci_tilde_mark = -1;
1663  INT ji, jj, i, j, index;
1664 
1665  fasp_iarray_set(row, cindex, -1);
1666 
1667  for ( i = 0; i < row; ++i ) {
1668 
1669  if ( vec[i] != FGPT ) continue; // skip non F-variables
1670 
1671  for ( ji = S->IA[i]; ji < S->IA[i+1]; ++ji ) {
1672  j = S->JA[ji];
1673  if ( vec[j] == CGPT ) cindex[j] = i; // mark C-neighbors
1674  else cindex[j] = -1; // reset cindex --Chensong 06/02/2013
1675  }
1676 
1677  if ( ci_tilde_mark != i ) ci_tilde = -1;//???
1678 
1679  for ( ji = S->IA[i]; ji < S->IA[i+1]; ++ji ) {
1680 
1681  j = S->JA[ji];
1682 
1683  if ( vec[j] != FGPT ) continue; // skip non F-variables
1684 
1685  // check whether there is a C-connection
1686  set_empty = TRUE;
1687  for ( jj = S->IA[j]; jj < S->IA[j+1]; ++jj ) {
1688  index = S->JA[jj];
1689  if ( cindex[index] == i ) {
1690  set_empty = FALSE; break;
1691  }
1692  } // end for jj
1693 
1694  // change the point i (if only F-F exists) to C
1695  if ( set_empty ) {
1696  if ( C_i_nonempty ) {
1697  vec[i] = CGPT; col++;
1698  if ( ci_tilde > -1 ) {
1699  vec[ci_tilde] = FGPT; col--;
1700  ci_tilde = -1;
1701  }
1702  C_i_nonempty = FALSE;
1703  break;
1704  }
1705  else { // temporary set j->C and roll back
1706  vec[j] = CGPT; col++;
1707  ci_tilde = j;
1708  ci_tilde_mark = i;
1709  C_i_nonempty = TRUE;
1710  i--; // roll back to check i-point again
1711  break;
1712  } // end if C_i_nonempty
1713  } // end if set_empty
1714 
1715  } // end for ji
1716 
1717  } // end for i
1718 
1719  fasp_mem_free(cindex);
1720 
1721  return col;
1722 }
1723 
1742 static void form_P_pattern_dir (dCSRmat *P,
1743  iCSRmat *S,
1744  ivector *vertices,
1745  INT row,
1746  INT col )
1747 {
1748  // local variables
1749  INT i, j, k, index;
1750  INT *vec = vertices->val;
1751 
1752  INT nthreads = 1, use_openmp = FALSE;
1753 
1754 #ifdef _OPENMP
1755  if ( row > OPENMP_HOLDS ) {
1756  use_openmp = TRUE;
1757  nthreads = FASP_GET_NUM_THREADS();
1758  }
1759 #endif
1760 
1761  // Initialize P matrix
1762  P->row = row; P->col = col;
1763  P->IA = (INT*)fasp_mem_calloc(row+1, sizeof(INT));
1764 
1765  // step 1: Find the structure IA of P first: using P as a counter
1766  if ( use_openmp ) {
1767 
1768  INT mybegin,myend,myid;
1769 #ifdef _OPENMP
1770 #pragma omp parallel for private(myid, mybegin,myend,i,j,k)
1771 #endif
1772  for ( myid = 0; myid < nthreads; myid++ ) {
1773  FASP_GET_START_END(myid, nthreads, row, &mybegin, &myend);
1774  for ( i = mybegin; i < myend; ++i ) {
1775  switch ( vec[i] ) {
1776  case FGPT: // fine grid points
1777  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1778  k = S->JA[j];
1779  if ( vec[k] == CGPT ) P->IA[i+1]++;
1780  }
1781  break;
1782 
1783  case CGPT: // coarse grid points
1784  P->IA[i+1] = 1; break;
1785 
1786  default: // treat everything else as isolated
1787  P->IA[i+1] = 0; break;
1788  }
1789  }
1790  }
1791 
1792  }
1793 
1794  else {
1795 
1796  for ( i = 0; i < row; ++i ) {
1797  switch ( vec[i] ) {
1798  case FGPT: // fine grid points
1799  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1800  k = S->JA[j];
1801  if ( vec[k] == CGPT ) P->IA[i+1]++;
1802  }
1803  break;
1804 
1805  case CGPT: // coarse grid points
1806  P->IA[i+1] = 1; break;
1807 
1808  default: // treat everything else as isolated
1809  P->IA[i+1] = 0; break;
1810  }
1811  } // end for i
1812 
1813  } // end if
1814 
1815  // Form P->IA from the counter P
1816  for ( i = 0; i < P->row; ++i ) P->IA[i+1] += P->IA[i];
1817  P->nnz = P->IA[P->row]-P->IA[0];
1818 
1819  // step 2: Find the structure JA of P
1820  P->JA = (INT*)fasp_mem_calloc(P->nnz,sizeof(INT));
1821  P->val = (REAL*)fasp_mem_calloc(P->nnz,sizeof(REAL));
1822 
1823  for ( index = i = 0; i < row; ++i ) {
1824  if ( vec[i] == FGPT ) { // fine grid point
1825  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1826  k = S->JA[j];
1827  if ( vec[k] == CGPT ) P->JA[index++] = k;
1828  } // end for j
1829  } // end if
1830  else if ( vec[i] == CGPT ) { // coarse grid point -- one entry only
1831  P->JA[index++] = i;
1832  }
1833  }
1834 
1835 }
1836 
1855 static void form_P_pattern_std (dCSRmat *P,
1856 iCSRmat *S,
1857 ivector *vertices,
1858 INT row,
1859 INT col)
1860 {
1861  // local variables
1862  INT i, j, k, l, h, index;
1863  INT *vec = vertices->val;
1864 
1865  // number of times a C-point is visited
1866  INT *visited = (INT*)fasp_mem_calloc(row,sizeof(INT));
1867 
1868  P->row = row; P->col = col;
1869  P->IA = (INT*)fasp_mem_calloc(row+1, sizeof(INT));
1870 
1871  fasp_iarray_set(row, visited, -1);
1872 
1873  // Step 1: Find the structure IA of P first: use P as a counter
1874  for ( i = 0; i < row; ++i ) {
1875 
1876  if ( vec[i] == FGPT ) { // if node i is a F point
1877  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1878 
1879  k = S->JA[j];
1880 
1881  // if neighbor of i is a C point, good
1882  if ( (vec[k] == CGPT) && (visited[k] != i) ) {
1883  visited[k] = i;
1884  P->IA[i+1]++;
1885  }
1886 
1887  // if k is a F point and k is not i, look for indirect C neighbors
1888  else if ( (vec[k] == FGPT) && (k != i) ) {
1889  for ( l = S->IA[k]; l < S->IA[k+1]; l++ ) { // neighbors of k
1890  h = S->JA[l];
1891  if ( (vec[h] == CGPT) && (visited[h] != i) ) {
1892  visited[h] = i;
1893  P->IA[i+1]++;
1894  }
1895  } // end for(l=S->IA[k];l<S->IA[k+1];l++)
1896  } // end if (vec[k]==CGPT)
1897 
1898  } // end for (j=S->IA[i];j<S->IA[i+1];j++)
1899  }
1900 
1901  else if ( vec[i] == CGPT ) { // if node i is a C point
1902  P->IA[i+1] = 1;
1903  }
1904 
1905  else { // treat everything else as isolated points
1906  P->IA[i+1] = 0;
1907  } // end if (vec[i]==FGPT)
1908 
1909  } // end for (i=0;i<row;++i)
1910 
1911  // Form P->IA from the counter P
1912  for ( i = 0; i < P->row; ++i ) P->IA[i+1] += P->IA[i];
1913  P->nnz = P->IA[P->row]-P->IA[0];
1914 
1915  // Step 2: Find the structure JA of P
1916  P->JA = (INT*)fasp_mem_calloc(P->nnz,sizeof(INT));
1917  P->val = (REAL*)fasp_mem_calloc(P->nnz,sizeof(REAL));
1918 
1919  fasp_iarray_set(row, visited, -1); // re-init visited array
1920 
1921  for ( i = 0; i < row; ++i ) {
1922 
1923  if ( vec[i] == FGPT ) { // if node i is a F point
1924 
1925  index = 0;
1926 
1927  for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1928 
1929  k = S->JA[j];
1930 
1931  // if neighbor k of i is a C point
1932  if ( (vec[k] == CGPT) && (visited[k] != i) ) {
1933  visited[k] = i;
1934  P->JA[P->IA[i]+index] = k;
1935  index++;
1936  }
1937 
1938  // if neighbor k of i is a F point and k is not i
1939  else if ( (vec[k] == FGPT) && (k != i) ) {
1940  for ( l = S->IA[k]; l < S->IA[k+1]; l++ ) { // neighbors of k
1941  h = S->JA[l];
1942  if ( (vec[h] == CGPT) && (visited[h] != i) ) {
1943  visited[h] = i;
1944  P->JA[P->IA[i]+index] = h;
1945  index++;
1946  }
1947 
1948  } // end for (l=S->IA[k];l<S->IA[k+1];l++)
1949 
1950  } // end if (vec[k]==CGPT)
1951 
1952  } // end for (j=S->IA[i];j<S->IA[i+1];j++)
1953  }
1954 
1955  else if ( vec[i] == CGPT ) {
1956  P->JA[P->IA[i]] = i;
1957  }
1958  }
1959 
1960  // clean up
1961  fasp_mem_free(visited);
1962 }
1963 
1978 static INT cfsplitting_mis (iCSRmat *S,
1979  ivector *vertices,
1980  ivector *order)
1981 {
1982  const INT n = S->row;
1983 
1984  INT col = 0;
1985  INT *ord = order->val;
1986  INT *vec = vertices->val;
1987  INT *IS = S->IA;
1988  INT *JS = S->JA;
1989 
1990  INT i, j, ind;
1991  INT row_begin, row_end;
1992 
1993  fasp_ivec_set (UNPT, vertices);
1994 
1995  for (i=0; i<n ; i++)
1996  {
1997  ind = ord[i];
1998  if (vec[ind] == UNPT) {
1999  vec[ind] = CGPT;
2000  row_begin = IS[ind]; row_end = IS[ind+1];
2001  for (j = row_begin; j<row_end; j++)
2002  {
2003  if (vec[JS[j]] == CGPT ) {
2004  vec[ind] = FGPT;
2005  break;
2006  }
2007  }
2008  if (vec[ind] == CGPT) {
2009  col++;
2010  for (j = row_begin; j<row_end; j++)
2011  {
2012  vec[JS[j]] = FGPT;
2013  }
2014  }
2015  }
2016  }
2017  return col;
2018 }
2019 
2034 static void ordering1 (iCSRmat *S,
2035  ivector *order)
2036 {
2037  const INT n = order->row;
2038  INT * IS = S->IA;
2039  INT * ord = order->val;
2040  INT maxind, maxdeg, degree;
2041 
2042  INT i;
2043  for (i = 0; i < n; i++)
2044  {
2045  ord[i] = i;
2046  }
2047  maxdeg=0;
2048  for (i = 0; i < n; i++)
2049  {
2050  degree = IS[i+1] - IS[i];
2051  if (degree > maxdeg)
2052  {
2053  maxind = i;
2054  maxdeg = degree;
2055  }
2056  }
2057  ord[0] = maxind;
2058  ord[maxind] = 0;
2059  return;
2060 }
2061 
2062 /*---------------------------------*/
2063 /*-- End of File --*/
2064 /*---------------------------------*/
#define TRUE
Definition of logic type.
Definition: fasp_const.h:67
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: message.c:199
Parameters for AMG solver.
Definition: fasp.h:583
#define ERROR_AMG_COARSEING
Definition: fasp_const.h:45
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
SHORT coarsening_type
coarsening type
Definition: fasp.h:643
struct linked_list * prev_node
previous node
Definition: fasp.h:1187
#define REAL
Definition: fasp.h:67
Vector with n entries of INT type.
Definition: fasp.h:356
INT head
starting of the list
Definition: fasp.h:1178
#define COARSE_AC
Definition: fasp_const.h:200
#define ERROR_UNKNOWN
Definition: fasp_const.h:62
#define ERROR_AMG_INTERP_TYPE
Definition: fasp_const.h:42
void fasp_dvec_free(dvector *u)
Free vector data space of REAL type.
Definition: vec.c:139
REAL * val
nonzero entries of A
Definition: fasp.h:166
#define INTERP_DIR
Definition of interpolation types.
Definition: fasp_const.h:206
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
INT aggressive_path
number of paths use to determine strongly coupled C points
Definition: fasp.h:664
Vector with n entries of REAL type.
Definition: fasp.h:342
#define INTERP_STD
Definition: fasp_const.h:207
#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
REAL strong_threshold
strong connection threshold for coarsening
Definition: fasp.h:652
INT row
row number of matrix A, m
Definition: fasp.h:181
void FASP_GET_START_END(INT procid, INT nprocs, INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: threads.c:83
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void fasp_iarray_cp(const INT n, INT *x, INT *y)
Copy an array to the other y=x.
Definition: array.c:185
void fasp_ivec_free(ivector *u)
Free vector data space of INT type.
Definition: vec.c:159
#define OPENMP_HOLDS
Definition: fasp_const.h:248
INT nnz
number of nonzero entries
Definition: fasp.h:157
void fasp_icsr_trans(iCSRmat *A, iCSRmat *AT)
Find transpose of iCSRmat matrix A.
Definition: sparse_csr.c:750
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
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
INT * val
nonzero entries of A
Definition: fasp.h:196
INT * val
actual vector entries
Definition: fasp.h:362
#define COARSE_RSP
Definition: fasp_const.h:198
INT row
row number of matrix A, m
Definition: fasp.h:151
struct linked_list * next_node
next node
Definition: fasp.h:1184
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
#define ERROR_AMG_COARSE_TYPE
Definition: fasp_const.h:44
void fasp_dcsr_getdiag(INT n, dCSRmat *A, dvector *diag)
Get first n diagonal entries of a CSR matrix A.
Definition: sparse_csr.c:410
void fasp_ivec_set(const INT m, ivector *u)
Set ivector value to be m.
Definition: vec.c:304
#define CGPT
Definition: fasp_const.h:216
#define COARSE_MIS
Definition: fasp_const.h:201
void fasp_ivec_alloc(const INT m, ivector *u)
Create vector data space of INT type.
Definition: vec.c:119
void fasp_icsr_free(iCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: sparse_csr.c:185
REAL max_row_sum
maximal row sum parameter
Definition: fasp.h:655
#define ABS(a)
Definition: fasp.h:74
#define UNPT
Definition: fasp_const.h:214
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define MIN(a, b)
Definition: fasp.h:73
SHORT interpolation_type
interpolation type
Definition: fasp.h:649
#define COARSE_RS
Definition of coarsening types.
Definition: fasp_const.h:197
INT count
INT row
number of rows
Definition: fasp.h:359
#define ISPT
Definition: fasp_const.h:217
#define INTERP_ENG
Definition: fasp_const.h:208
INT nnz
number of nonzero entries
Definition: fasp.h:187
#define FALSE
Definition: fasp_const.h:68
INT col
column of matrix A, n
Definition: fasp.h:184
#define COARSE_CR
Definition: fasp_const.h:199
ivector fasp_ivec_create(const INT m)
Create vector data space of INT type.
Definition: vec.c:78
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
Definition: array.c:107
INT fasp_amg_coarsening_cr(const INT i_0, const INT i_n, dCSRmat *A, ivector *vertices, AMG_param *param)
CR coarsening.
Definition: coarsening_cr.c:42
Sparse matrix of INT type in CSR format.
Definition: fasp.h:178
A linked list node.
Definition: fasp.h:1171
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72
#define FGPT
Definition: fasp_const.h:215