Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
PreAMGCoarsenRS.c
Go to the documentation of this file.
1
20#ifdef _OPENMP
21#include <omp.h>
22#endif
23
24#include "fasp.h"
25#include "fasp_functs.h"
26
27/*---------------------------------*/
28/*-- Declare Private Functions --*/
29/*---------------------------------*/
30
31#include "PreAMGUtil.inl"
32
33static INT cfsplitting_cls (dCSRmat *, iCSRmat *, ivector *);
34static INT cfsplitting_clsp (dCSRmat *, iCSRmat *, ivector *);
35static INT cfsplitting_agg (dCSRmat *, iCSRmat *, ivector *, INT);
36static INT cfsplitting_mis (iCSRmat *, ivector *, ivector *);
37static INT clean_ff_couplings (iCSRmat *, ivector *, INT, INT);
38static INT compress_S (iCSRmat *);
39
40static void strong_couplings (dCSRmat *, iCSRmat *, AMG_param *);
41static void form_P_pattern_dir (dCSRmat *, iCSRmat *, ivector *, INT, INT);
42static void form_P_pattern_std (dCSRmat *, iCSRmat *, ivector *, INT, INT);
43static void ordering1 (iCSRmat *, ivector *);
44
45/*---------------------------------*/
46/*-- Public Functions --*/
47/*---------------------------------*/
48
74 ivector *vertices,
75 dCSRmat *P,
76 iCSRmat *S,
77 AMG_param *param)
78{
79 const SHORT coarse_type = param->coarsening_type;
80 const INT agg_path = param->aggressive_path;
81 const INT row = A->row;
82
83 // local variables
84 SHORT interp_type = param->interpolation_type;
85 INT col = 0;
86
87#if DEBUG_MODE > 0
88 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
89#endif
90
91#if DEBUG_MODE > 1
92 printf("### DEBUG: Step 1. Find strong connections ......\n");
93#endif
94
95 // make sure standard interp is used for aggressive coarsening
96 if ( coarse_type == COARSE_AC ) interp_type = INTERP_STD;
97
98 // find strong couplings and return them in S
99 strong_couplings(A, S, param);
100
101#if DEBUG_MODE > 1
102 printf("### DEBUG: Step 2. C/F splitting ......\n");
103#endif
104
105 switch ( coarse_type ) {
106
107 case COARSE_RSP: // Classical coarsening with positive connections
108 col = cfsplitting_clsp(A, S, vertices); break;
109
110 case COARSE_AC: // Aggressive coarsening
111 col = cfsplitting_agg(A, S, vertices, agg_path); break;
112
113 case COARSE_CR: // Compatible relaxation
114 col = fasp_amg_coarsening_cr(0, row-1, A, vertices, param); break;
115
116 case COARSE_MIS: // Maximal independent set
117 {
118 ivector order = fasp_ivec_create(row);
119 compress_S(S);
120 ordering1(S, &order);
121 col = cfsplitting_mis(S, vertices, &order);
122 fasp_ivec_free(&order);
123 break;
124 }
125
126 default: // Classical coarsening
127 col = cfsplitting_cls(A, S, vertices);
128
129 }
130
131#if DEBUG_MODE > 1
132 printf("### DEBUG: col = %d\n", col);
133#endif
134 if ( col <= 0 ) return ERROR_UNKNOWN;
135
136#if DEBUG_MODE > 1
137 printf("### DEBUG: Step 3. Find support of C points ......\n");
138#endif
139
140 switch ( interp_type ) {
141
142 case INTERP_DIR: // Direct interpolation or ...
143 case INTERP_ENG: // Energy-min interpolation
144 col = clean_ff_couplings(S, vertices, row, col);
145 form_P_pattern_dir(P, S, vertices, row, col);
146 break;
147
148 case INTERP_STD: // Standard interpolation
149 case INTERP_EXT: // Extended interpolation
150 form_P_pattern_std(P, S, vertices, row, col); break;
151
152 default:
153 fasp_chkerr(ERROR_AMG_INTERP_TYPE, __FUNCTION__);
154
155 }
156
157#if DEBUG_MODE > 0
158 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
159#endif
160
161 return FASP_SUCCESS;
162}
163
164
165/*---------------------------------*/
166/*-- Private Functions --*/
167/*---------------------------------*/
168
186static void strong_couplings (dCSRmat *A,
187 iCSRmat *S,
188 AMG_param *param )
189{
190 const SHORT coarse_type = param->coarsening_type;
191 const REAL max_row_sum = param->max_row_sum;
192 const REAL epsilon_str = param->strong_threshold;
193 const INT row = A->row, col = A->col, row1 = row+1;
194 const INT nnz = A->nnz;
195
196 INT *ia = A->IA, *ja = A->JA;
197 REAL *aj = A->val;
198
199 // local variables
200 INT i, j, begin_row, end_row;
201 REAL row_scl, row_sum;
202
203 SHORT nthreads = 1, use_openmp = FALSE;
204
205#ifdef _OPENMP
206 if ( row > OPENMP_HOLDS ) {
207 use_openmp = TRUE;
208 nthreads = fasp_get_num_threads();
209 }
210#endif
211
212 // get the diagonal entry of A: assume all connections are strong
213 dvector diag; fasp_dcsr_getdiag(0, A, &diag);
214
215 // copy the structure of A to S
216 S->row = row; S->col = col; S->nnz = nnz; S->val = NULL;
217 S->IA = (INT *)fasp_mem_calloc(row1, sizeof(INT));
218 S->JA = (INT *)fasp_mem_calloc(nnz, sizeof(INT));
219 fasp_iarray_cp(row1, ia, S->IA);
220 fasp_iarray_cp(nnz, ja, S->JA);
221
222 if ( use_openmp ) {
223
224 // This part is still old! Need to be updated. --Chensong 09/18/2016
225
226 INT mybegin, myend, myid;
227#ifdef _OPENMP
228#pragma omp parallel for private(myid, mybegin,myend,i,row_scl,row_sum,begin_row,end_row,j)
229#endif
230 for ( myid = 0; myid < nthreads; myid++ ) {
231 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
232 for ( i = mybegin; i < myend; i++) {
233
234 // Compute most negative entry in each row and row sum
235 row_scl = row_sum = 0.0;
236 begin_row = ia[i]; end_row = ia[i+1];
237 for ( j = begin_row; j < end_row; j++ ) {
238 row_scl = MIN(row_scl, aj[j]);
239 row_sum += aj[j];
240 }
241
242 // Find diagonal entries of S and remove them later
243 for ( j = begin_row; j < end_row; j++ ) {
244 if ( ja[j] == i ) { S->JA[j] = -1; break; }
245 }
246
247 // Mark entire row as weak couplings if strongly diagonal-dominant
248 if ( ABS(row_sum) > max_row_sum * ABS(diag.val[i]) ) {
249 for ( j = begin_row; j < end_row; j++ ) S->JA[j] = -1;
250 }
251 else {
252 for ( j = begin_row; j < end_row; j++) {
253 // If a_{ij} >= \epsilon_{str} * \min a_{ij}, the connection
254 // j->i is set to be weak; positive entries result in weak
255 // connections
256 if ( A->val[j] >= epsilon_str*row_scl ) S->JA[j] = -1;
257 }
258 }
259
260 } // end for i
261 } // end for myid
262
263 }
264
265 else {
266
267 for ( i = 0; i < row; ++i ) {
268
269 // Compute row scale and row sum
270 row_scl = row_sum = 0.0;
271 begin_row = ia[i]; end_row = ia[i+1];
272
273 for ( j = begin_row; j < end_row; j++ ) {
274
275 // Originally: Not consider positive entries
276 // row_sum += aj[j];
277 // Now changed to --Chensong 05/17/2013
278 row_sum += ABS(aj[j]);
279
280 // Originally: Not consider positive entries
281 // row_scl = MAX(row_scl, -aj[j]); // smallest negative
282 // Now changed to --Chensong 06/01/2013
283 if ( ja[j] != i ) row_scl = MAX(row_scl, ABS(aj[j])); // largest abs
284
285 }
286
287 // Multiply by the strength threshold
288 row_scl *= epsilon_str;
289
290 // Find diagonal entries of S and remove them later
291 for ( j = begin_row; j < end_row; j++ ) {
292 if ( ja[j] == i ) { S->JA[j] = -1; break; }
293 }
294
295 // Mark entire row as weak couplings if strongly diagonal-dominant
296 // Originally: Not consider positive entries
297 // if ( ABS(row_sum) > max_row_sum * ABS(diag.val[i]) ) {
298 // Now changed to --Chensong 05/17/2013
299 if ( row_sum < (2 - max_row_sum) * ABS(diag.val[i]) ) {
300
301 for ( j = begin_row; j < end_row; j++ ) S->JA[j] = -1;
302
303 }
304 else {
305
306 switch ( coarse_type ) {
307
308 case COARSE_RSP: // consider positive off-diag as well
309 for ( j = begin_row; j < end_row; j++ ) {
310 if ( ABS(A->val[j]) <= row_scl ) S->JA[j] = -1;
311 }
312 break;
313
314 default: // only consider n-couplings
315 for ( j = begin_row; j < end_row; j++ ) {
316 if ( -A->val[j] <= row_scl ) S->JA[j] = -1;
317 }
318 break;
319
320 }
321 }
322 } // end for i
323
324 } // end if openmp
325
326 fasp_dvec_free(&diag);
327}
328
343static INT compress_S (iCSRmat *S)
344{
345 const INT row = S->row;
346 INT * ia = S->IA;
347
348 // local variables
349 INT index, i, j, begin_row, end_row;
350
351 // compress S: remove weak connections and form strong coupling matrix
352 for ( index = i = 0; i < row; ++i ) {
353
354 begin_row = ia[i]; end_row = ia[i+1];
355
356 ia[i] = index;
357 for ( j = begin_row; j < end_row; j++ ) {
358 if ( S->JA[j] > -1 ) S->JA[index++] = S->JA[j]; // strong couplings
359 }
360
361 }
362
363 S->nnz = S->IA[row] = index;
364
365 if ( S->nnz <= 0 ) {
366 return ERROR_UNKNOWN;
367 }
368 else {
369 return FASP_SUCCESS;
370 }
371}
372
385static void rem_positive_ff (dCSRmat *A,
386 iCSRmat *Stemp,
387 ivector *vertices)
388{
389 const INT row = A->row;
390 INT *ia = A->IA, *vec = vertices->val;
391
392 REAL row_scl, max_entry;
393 INT i, j, ji, max_index;
394
395 for ( i = 0; i < row; ++i ) {
396
397 if ( vec[i] != FGPT ) continue; // skip non F-variables
398
399 row_scl = 0.0;
400 for ( ji = ia[i]; ji < ia[i+1]; ++ji ) {
401 j = A->JA[ji];
402 if ( j == i ) continue; // skip diagonal
403 row_scl = MAX(row_scl, ABS(A->val[ji])); // max abs entry
404 } // end for ji
405 row_scl *= 0.75;
406
407 // looking for strong F-F connections
408 max_index = -1; max_entry = 0.0;
409 for ( ji = ia[i]; ji < ia[i+1]; ++ji ) {
410 j = A->JA[ji];
411 if ( j == i ) continue; // skip diagonal
412 if ( vec[j] != FGPT ) continue; // skip F-C connections
413 if ( A->val[ji] > row_scl ) {
414 Stemp->JA[ji] = j;
415 if ( A->val[ji] > max_entry ) {
416 max_entry = A->val[ji];
417 max_index = j; // max positive entry
418 }
419 }
420 } // end for ji
421
422 // mark max positive entry as C-point
423 if ( max_index != -1 ) vec[max_index] = CGPT;
424
425 } // end for i
426
427}
428
450static INT cfsplitting_cls (dCSRmat *A,
451 iCSRmat *S,
452 ivector *vertices)
453{
454 const INT row = A->row;
455
456 // local variables
457 INT col = 0;
458 INT maxmeas, maxnode, num_left = 0;
459 INT measure, newmeas;
460 INT i, j, k, l;
461 INT myid, mybegin, myend;
462 INT *vec = vertices->val;
463 INT *work = (INT*)fasp_mem_calloc(3*row,sizeof(INT));
464 INT *lists = work, *where = lists+row, *lambda = where+row;
465
466#if RS_C1
467 INT set_empty = 1;
468 INT jkeep = 0, cnt, index;
469 INT row_end_S, ji, row_end_S_nabor, jj;
470 INT *graph_array = lambda;
471#else
472 INT *ia = A->IA;
473#endif
474
475 LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
476
477 SHORT nthreads = 1, use_openmp = FALSE;
478
479#if DEBUG_MODE > 0
480 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
481#endif
482
483#ifdef _OPENMP
484 if ( row > OPENMP_HOLDS ) {
485 use_openmp = TRUE;
486 nthreads = fasp_get_num_threads();
487 }
488#endif
489
490 // 0. Compress S and form S_transpose
491 col = compress_S(S);
492 if ( col < 0 ) goto FINISHED; // compression failed!!!
493
494 iCSRmat ST; fasp_icsr_trans(S, &ST);
495
496 // 1. Initialize lambda
497 if ( use_openmp ) {
498#ifdef _OPENMP
499#pragma omp parallel for private(myid, mybegin, myend, i)
500#endif
501 for ( myid = 0; myid < nthreads; myid++ ) {
502 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
503 for ( i = mybegin; i < myend; i++ ) lambda[i] = ST.IA[i+1] - ST.IA[i];
504 }
505 }
506 else {
507 for ( i = 0; i < row; ++i ) lambda[i] = ST.IA[i+1] - ST.IA[i];
508 }
509
510 // 2. Before C/F splitting algorithm starts, filter out the variables which
511 // have no connections at all and mark them as special F-variables.
512 if ( use_openmp ) {
513
514#ifdef _OPENMP
515#pragma omp parallel for reduction(+:num_left) private(myid, mybegin, myend, i)
516#endif
517 for ( myid = 0; myid < nthreads; myid++ ) {
518 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
519 for ( i = mybegin; i < myend; i++ ) {
520#if RS_C1 // Check C1 criteria or not
521 if ( S->IA[i+1] == S->IA[i] )
522#else
523 if ( (ia[i+1]-ia[i]) <= 1 )
524#endif
525 {
526 vec[i] = ISPT; // set i as an ISOLATED fine node
527 lambda[i] = 0;
528 }
529 else {
530 vec[i] = UNPT; // set i as a undecided node
531 num_left++;
532 }
533 }
534 } // end for myid
535
536 }
537
538 else {
539
540 for ( i = 0; i < row; ++i ) {
541
542#if RS_C1
543 if ( S->IA[i+1] == S->IA[i] )
544#else
545 if ( (ia[i+1]-ia[i]) <= 1 )
546#endif
547 {
548 vec[i] = ISPT; // set i as an ISOLATED fine node
549 lambda[i] = 0;
550 }
551 else {
552 vec[i] = UNPT; // set i as a undecided node
553 num_left++;
554 }
555 } // end for i
556
557 }
558
559 // 3. Form linked list for lambda (max to min)
560 for ( i = 0; i < row; ++i ) {
561
562 if ( vec[i] == ISPT ) continue; // skip isolated variables
563
564 measure = lambda[i];
565
566 if ( measure > 0 ) {
567 enter_list(&LoL_head, &LoL_tail, lambda[i], i, lists, where);
568 }
569 else {
570
571 if ( measure < 0 ) printf("### WARNING: Negative lambda[%d]!\n", i);
572
573 // Set variables with non-positive measure as F-variables
574 vec[i] = FGPT; // no strong connections, set i as fine node
575 --num_left;
576
577 // Update lambda and linked list after i->F
578 for ( k = S->IA[i]; k < S->IA[i+1]; ++k ) {
579 j = S->JA[k];
580 if ( vec[j] == ISPT ) continue; // skip isolate variables
581 if ( j < i ) {
582 newmeas = lambda[j];
583 if ( newmeas > 0 ) {
584 remove_node(&LoL_head, &LoL_tail, newmeas, j, lists, where);
585 }
586 newmeas = ++(lambda[j]);
587 enter_list(&LoL_head, &LoL_tail, newmeas, j, lists, where);
588 }
589 else {
590 newmeas = ++(lambda[j]);
591 }
592 }
593
594 } // end if measure
595
596 } // end for i
597
598 // 4. Main loop
599 while ( num_left > 0 ) {
600
601 // pick $i\in U$ with $\max\lambda_i: C:=C\cup\{i\}, U:=U\\{i\}$
602 maxnode = LoL_head->head;
603 maxmeas = lambda[maxnode];
604 if ( maxmeas == 0 )
605 printf("### WARNING: Head of the list has measure 0!\n");
606
607 vec[maxnode] = CGPT; // set maxnode as coarse node
608 lambda[maxnode] = 0;
609 --num_left;
610 remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
611 col++;
612
613 // for all $j\in S_i^T\cap U: F:=F\cup\{j\}, U:=U\backslash\{j\}$
614 for ( i = ST.IA[maxnode]; i < ST.IA[maxnode+1]; ++i ) {
615
616 j = ST.JA[i];
617
618 if ( vec[j] != UNPT ) continue; // skip decided variables
619
620 vec[j] = FGPT; // set j as fine node
621 remove_node(&LoL_head, &LoL_tail, lambda[j], j, lists, where);
622 --num_left;
623
624 // Update lambda and linked list after j->F
625 for ( l = S->IA[j]; l < S->IA[j+1]; l++ ) {
626 k = S->JA[l];
627 if ( vec[k] == UNPT ) { // k is unknown
628 remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
629 newmeas = ++(lambda[k]);
630 enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
631 }
632 }
633
634 } // end for i
635
636 // Update lambda and linked list after maxnode->C
637 for ( i = S->IA[maxnode]; i < S->IA[maxnode+1]; ++i ) {
638
639 j = S->JA[i];
640
641 if ( vec[j] != UNPT ) continue; // skip decided variables
642
643 measure = lambda[j];
644 remove_node(&LoL_head, &LoL_tail, measure, j, lists, where);
645 lambda[j] = --measure;
646
647 if ( measure > 0 ) {
648 enter_list(&LoL_head, &LoL_tail, measure, j, lists, where);
649 }
650 else { // j is the only point left, set as fine variable
651 vec[j] = FGPT;
652 --num_left;
653
654 // Update lambda and linked list after j->F
655 for ( l = S->IA[j]; l < S->IA[j+1]; l++ ) {
656 k = S->JA[l];
657 if ( vec[k] == UNPT ) { // k is unknown
658 remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
659 newmeas = ++(lambda[k]);
660 enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
661 }
662 } // end for l
663 } // end if
664
665 } // end for
666
667 } // end while
668
669#if RS_C1
670
671 // C/F splitting of RS coarsening check C1 Criterion
672 fasp_iarray_set(row, graph_array, -1);
673 for (i = 0; i < row; i ++)
674 {
675 if (vec[i] == FGPT)
676 {
677 row_end_S = S->IA[i+1];
678 for (ji = S->IA[i]; ji < row_end_S; ji ++)
679 {
680 j = S->JA[ji];
681 if (vec[j] == CGPT)
682 {
683 graph_array[j] = i;
684 }
685 }
686 cnt = 0;
687 for (ji = S->IA[i]; ji < row_end_S; ji ++)
688 {
689 j = S->JA[ji];
690 if (vec[j] == FGPT)
691 {
692 set_empty = 1;
693 row_end_S_nabor = S->IA[j+1];
694 for (jj = S->IA[j]; jj < row_end_S_nabor; jj ++)
695 {
696 index = S->JA[jj];
697 if (graph_array[index] == i)
698 {
699 set_empty = 0;
700 break;
701 }
702 }
703 if (set_empty)
704 {
705 if (cnt == 0)
706 {
707 vec[j] = CGPT;
708 col++;
709 graph_array[j] = i;
710 jkeep = j;
711 cnt = 1;
712 }
713 else
714 {
715 vec[i] = CGPT;
716 vec[jkeep] = FGPT;
717 break;
718 }
719 }
720 }
721 }
722 }
723 }
724
725#endif
726
727 fasp_icsr_free(&ST);
728
729 if ( LoL_head ) {
730 list_ptr = LoL_head;
731 LoL_head->prev_node = NULL;
732 LoL_head->next_node = NULL;
733 LoL_head = list_ptr->next_node;
734 fasp_mem_free(list_ptr); list_ptr = NULL;
735 }
736
737FINISHED:
738 fasp_mem_free(work); work = NULL;
739
740#if DEBUG_MODE > 0
741 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
742#endif
743
744 return col;
745}
746
766static INT cfsplitting_clsp (dCSRmat *A,
767 iCSRmat *S,
768 ivector *vertices)
769{
770 const INT row = A->row;
771
772 // local variables
773 INT col = 0;
774 INT maxmeas, maxnode, num_left = 0;
775 INT measure, newmeas;
776 INT i, j, k, l;
777 INT myid, mybegin, myend;
778
779 INT *ia = A->IA, *vec = vertices->val;
780 INT *work = (INT*)fasp_mem_calloc(3*row,sizeof(INT));
781 INT *lists = work, *where = lists+row, *lambda = where+row;
782
783 LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
784
785 SHORT nthreads = 1, use_openmp = FALSE;
786
787#if DEBUG_MODE > 0
788 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
789#endif
790
791#ifdef _OPENMP
792 if ( row > OPENMP_HOLDS ) {
793 use_openmp = TRUE;
794 nthreads = fasp_get_num_threads();
795 }
796#endif
797
798 // 0. Compress S and form S_transpose (not complete, just IA and JA)
799 iCSRmat Stemp;
800 Stemp.row = S->row; Stemp.col = S->col; Stemp.nnz = S->nnz;
801 Stemp.IA = (INT *)fasp_mem_calloc(S->row+1, sizeof(INT));
802 fasp_iarray_cp (S->row+1, S->IA, Stemp.IA);
803 Stemp.JA = (INT *)fasp_mem_calloc(S->nnz, sizeof(INT));
804 fasp_iarray_cp (S->nnz, S->JA, Stemp.JA);
805
806 if ( compress_S(S) < 0 ) goto FINISHED; // compression failed!!!
807
808 iCSRmat ST; fasp_icsr_trans(S, &ST);
809
810 // 1. Initialize lambda
811 if ( use_openmp ) {
812#ifdef _OPENMP
813#pragma omp parallel for private(myid, mybegin,myend,i)
814#endif
815 for ( myid = 0; myid < nthreads; myid++ ) {
816 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
817 for ( i = mybegin; i < myend; i++ ) lambda[i] = ST.IA[i+1] - ST.IA[i];
818 }
819 }
820 else {
821 for ( i = 0; i < row; ++i ) lambda[i] = ST.IA[i+1] - ST.IA[i];
822 }
823
824 // 2. Before C/F splitting algorithm starts, filter out the variables which
825 // have no connections at all and mark them as special F-variables.
826 if ( use_openmp ) {
827
828#ifdef _OPENMP
829#pragma omp parallel for reduction(+:num_left) private(myid, mybegin, myend, i)
830#endif
831 for ( myid = 0; myid < nthreads; myid++ ) {
832 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
833 for ( i = mybegin; i < myend; i++ ) {
834 if ( (ia[i+1]-ia[i]) <= 1 ) {
835 vec[i] = ISPT; // set i as an ISOLATED fine node
836 lambda[i] = 0;
837 }
838 else {
839 vec[i] = UNPT; // set i as a undecided node
840 num_left++;
841 }
842 }
843 } // end for myid
844
845 }
846 else {
847
848 for ( i = 0; i < row; ++i ) {
849 if ( (ia[i+1]-ia[i]) <= 1 ) {
850 vec[i] = ISPT; // set i as an ISOLATED fine node
851 lambda[i] = 0;
852 }
853 else {
854 vec[i] = UNPT; // set i as a undecided node
855 num_left++;
856 }
857 } // end for i
858
859 }
860
861 // 3. Form linked list for lambda (max to min)
862 for ( i = 0; i < row; ++i ) {
863
864 if ( vec[i] == ISPT ) continue; // skip isolated variables
865
866 measure = lambda[i];
867
868 if ( measure > 0 ) {
869 enter_list(&LoL_head, &LoL_tail, lambda[i], i, lists, where);
870 }
871 else {
872
873 if ( measure < 0 ) printf("### WARNING: Negative lambda[%d]!\n", i);
874
875 // Set variables with non-positive measure as F-variables
876 vec[i] = FGPT; // no strong connections, set i as fine node
877 --num_left;
878
879 // Update lambda and linked list after i->F
880 for ( k = S->IA[i]; k < S->IA[i+1]; ++k ) {
881
882 j = S->JA[k];
883 if ( vec[j] == ISPT ) continue; // skip isolate variables
884
885 if ( j < i ) { // only look at the previous points!!
886 newmeas = lambda[j];
887 if ( newmeas > 0 ) {
888 remove_node(&LoL_head, &LoL_tail, newmeas, j, lists, where);
889 }
890 newmeas = ++(lambda[j]);
891 enter_list(&LoL_head, &LoL_tail, newmeas, j, lists, where);
892 }
893 else { // will be checked later on
894 newmeas = ++(lambda[j]);
895 } // end if
896
897 } // end for k
898
899 } // end if measure
900
901 } // end for i
902
903 // 4. Main loop
904 while ( num_left > 0 ) {
905
906 // pick $i\in U$ with $\max\lambda_i: C:=C\cup\{i\}, U:=U\\{i\}$
907 maxnode = LoL_head->head;
908 maxmeas = lambda[maxnode];
909 if ( maxmeas == 0 )
910 printf("### WARNING: Head of the list has measure 0!\n");
911
912 vec[maxnode] = CGPT; // set maxnode as coarse node
913 lambda[maxnode] = 0;
914 --num_left;
915 remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
916 col++;
917
918 // for all $j\in S_i^T\cap U: F:=F\cup\{j\}, U:=U\backslash\{j\}$
919 for ( i = ST.IA[maxnode]; i < ST.IA[maxnode+1]; ++i ) {
920
921 j = ST.JA[i];
922
923 if ( vec[j] != UNPT ) continue; // skip decided variables
924
925 vec[j] = FGPT; // set j as fine node
926 remove_node(&LoL_head, &LoL_tail, lambda[j], j, lists, where);
927 --num_left;
928
929 // Update lambda and linked list after j->F
930 for ( l = S->IA[j]; l < S->IA[j+1]; l++ ) {
931 k = S->JA[l];
932 if ( vec[k] == UNPT ) { // k is unknown
933 remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
934 newmeas = ++(lambda[k]);
935 enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
936 }
937 }
938
939 } // end for i
940
941 // Update lambda and linked list after maxnode->C
942 for ( i = S->IA[maxnode]; i < S->IA[maxnode+1]; ++i ) {
943
944 j = S->JA[i];
945
946 if ( vec[j] != UNPT ) continue; // skip decided variables
947
948 measure = lambda[j];
949 remove_node(&LoL_head, &LoL_tail, measure, j, lists, where);
950 lambda[j] = --measure;
951
952 if ( measure > 0 ) {
953 enter_list(&LoL_head, &LoL_tail, measure, j, lists, where);
954 }
955 else { // j is the only point left, set as fine variable
956 vec[j] = FGPT;
957 --num_left;
958
959 // Update lambda and linked list after j->F
960 for ( l = S->IA[j]; l < S->IA[j+1]; l++ ) {
961 k = S->JA[l];
962 if ( vec[k] == UNPT ) { // k is unknown
963 remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
964 newmeas = ++(lambda[k]);
965 enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
966 }
967 } // end for l
968 } // end if
969
970 } // end for
971
972 } // end while
973
974 fasp_icsr_free(&ST);
975
976 if ( LoL_head ) {
977 list_ptr = LoL_head;
978 LoL_head->prev_node = NULL;
979 LoL_head->next_node = NULL;
980 LoL_head = list_ptr->next_node;
981 fasp_mem_free(list_ptr); list_ptr = NULL;
982 }
983
984 // Enforce F-C connections. Adding this step helps for the ExxonMobil test
985 // problems! Need more tests though --Chensong 06/08/2013
986 // col = clean_ff_couplings(S, vertices, row, col);
987
988 rem_positive_ff(A, &Stemp, vertices);
989
990 if ( compress_S(&Stemp) < 0 ) goto FINISHED; // compression failed!!!
991
992 S->row = Stemp.row;
993 S->col = Stemp.col;
994 S->nnz = Stemp.nnz;
995
996 fasp_mem_free(S->IA); S->IA = Stemp.IA;
997 fasp_mem_free(S->JA); S->JA = Stemp.JA;
998
999FINISHED:
1000 fasp_mem_free(work); work = NULL;
1001
1002#if DEBUG_MODE > 0
1003 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1004#endif
1005
1006 return col;
1007}
1008
1029static void strong_couplings_agg1 (dCSRmat *A,
1030 iCSRmat *S,
1031 iCSRmat *Sh,
1032 ivector *vertices,
1033 ivector *CGPT_index,
1034 ivector *CGPT_rindex)
1035{
1036 const INT row = A->row;
1037
1038 // local variables
1039 INT i, j, k;
1040 INT num_c, count, ci, cj, ck, fj, cck;
1041 INT *cp_index, *cp_rindex, *visited;
1042 INT *vec = vertices->val;
1043
1044 // count the number of coarse grid points
1045 for ( num_c = i = 0; i < row; i++ ) {
1046 if ( vec[i] == CGPT ) num_c++;
1047 }
1048
1049 // for the reverse indexing of coarse grid points
1050 fasp_ivec_alloc(row, CGPT_rindex);
1051 cp_rindex = CGPT_rindex->val;
1052
1053 // generate coarse grid point index
1054 fasp_ivec_alloc(num_c, CGPT_index);
1055 cp_index = CGPT_index->val;
1056 for ( j = i = 0; i < row; i++ ) {
1057 if ( vec[i] == CGPT ) {
1058 cp_index[j] = i;
1059 cp_rindex[i] = j;
1060 j++;
1061 }
1062 }
1063
1064 // allocate space for Sh
1065 Sh->row = Sh->col = num_c;
1066 Sh->val = Sh->JA = NULL;
1067 Sh->IA = (INT*)fasp_mem_calloc(Sh->row+1, sizeof(INT));
1068
1069 // record the number of times some coarse point is visited
1070 visited = (INT*)fasp_mem_calloc(num_c, sizeof(INT));
1071 fasp_iarray_set(num_c, visited, -1);
1072
1073 /**********************************************/
1074 /* step 1: Find first the structure IA of Sh */
1075 /**********************************************/
1076
1077 Sh->IA[0] = 0;
1078
1079 for ( ci = 0; ci < Sh->row; ci++ ) {
1080
1081 i = cp_index[ci]; // find the index of the ci-th coarse grid point
1082
1083 // number of coarse point that i is strongly connected to w.r.t. S(p,2)
1084 count = 0;
1085
1086 // visit all the fine neighbors that ci is strongly connected to
1087 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1088
1089 fj = S->JA[j];
1090
1091 if ( vec[fj] == CGPT && fj != i ) {
1092 cj = cp_rindex[fj];
1093 if ( visited[cj] != ci ) {
1094 visited[cj] = ci; // mark as strongly connected from ci
1095 count++;
1096 }
1097
1098 }
1099
1100 else if ( vec[fj] == FGPT ) { // fine grid point,
1101
1102 // find all the coarse neighbors that fj is strongly connected to
1103 for ( k = S->IA[fj]; k < S->IA[fj+1]; k++ ) {
1104 ck = S->JA[k];
1105 if ( vec[ck] == CGPT && ck != i ) { // it is a coarse grid point
1106 if ( cp_rindex[ck] >= num_c ) {
1107 printf("### ERROR: ck=%d, num_c=%d, out of bound!\n",
1108 ck, num_c);
1109 fasp_chkerr(ERROR_AMG_COARSEING, __FUNCTION__);
1110 }
1111 cck = cp_rindex[ck];
1112
1113 if ( visited[cck] != ci ) {
1114 visited[cck] = ci; // mark as strongly connected from ci
1115 count++;
1116 }
1117 } //end if
1118 } //end for k
1119
1120 } //end if
1121
1122 } //end for j
1123
1124 Sh->IA[ci+1] = Sh->IA[ci] + count;
1125
1126 } //end for i
1127
1128 /*************************/
1129 /* step 2: Find JA of Sh */
1130 /*************************/
1131
1132 fasp_iarray_set(num_c, visited, -1); // reset visited
1133
1134 Sh->nnz = Sh->IA[Sh->row];
1135 Sh->JA = (INT*)fasp_mem_calloc(Sh->nnz, sizeof(INT));
1136
1137 for ( ci = 0; ci < Sh->row; ci++ ) {
1138
1139 i = cp_index[ci]; // find the index of the i-th coarse grid point
1140 count = Sh->IA[ci]; // count for coarse points
1141
1142 // visit all the fine neighbors that ci is strongly connected to
1143 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1144
1145 fj = S->JA[j];
1146
1147 if ( vec[fj] == CGPT && fj != i ) {
1148 cj = cp_rindex[fj];
1149 if ( visited[cj] != ci ) { // not visited yet
1150 visited[cj] = ci;
1151 Sh->JA[count] = cj;
1152 count++;
1153 }
1154 }
1155 else if ( vec[fj] == FGPT ) { // fine grid point,
1156 //find all the coarse neighbors that fj is strongly connected to
1157 for ( k = S->IA[fj]; k < S->IA[fj+1]; k++ ) {
1158 ck = S->JA[k];
1159 if ( vec[ck] == CGPT && ck != i ) { // coarse grid point
1160 cck = cp_rindex[ck];
1161 if ( visited[cck] != ci ) { // not visited yet
1162 visited[cck] = ci;
1163 Sh->JA[count] = cck;
1164 count++;
1165 }
1166 } // end if
1167 } // end for k
1168 } // end if
1169
1170 } // end for j
1171
1172 if ( count != Sh->IA[ci+1] ) {
1173 printf("### WARNING: Inconsistent numbers of nonzeros!\n ");
1174 }
1175
1176 } // end for ci
1177
1178 fasp_mem_free(visited); visited = NULL;
1179}
1180
1207static void strong_couplings_agg2 (dCSRmat *A,
1208 iCSRmat *S,
1209 iCSRmat *Sh,
1210 ivector *vertices,
1211 ivector *CGPT_index,
1212 ivector *CGPT_rindex)
1213{
1214 const INT row = A->row;
1215
1216 // local variables
1217 INT i, j, k;
1218 INT num_c, count, ci, cj, ck, fj, cck;
1219 INT *cp_index, *cp_rindex, *visited;
1220 INT *vec = vertices->val;
1221
1222 // count the number of coarse grid points
1223 for ( num_c = i = 0; i < row; i++ ) {
1224 if ( vec[i] == CGPT ) num_c++;
1225 }
1226
1227 // for the reverse indexing of coarse grid points
1228 fasp_ivec_alloc(row, CGPT_rindex);
1229 cp_rindex = CGPT_rindex->val;
1230
1231 // generate coarse grid point index
1232 fasp_ivec_alloc(num_c, CGPT_index);
1233 cp_index = CGPT_index->val;
1234 for ( j = i = 0; i < row; i++ ) {
1235 if ( vec[i] == CGPT ) {
1236 cp_index[j] = i;
1237 cp_rindex[i] = j;
1238 j++;
1239 }
1240 }
1241
1242 // allocate space for Sh
1243 Sh->row = Sh->col = num_c;
1244 Sh->val = Sh->JA = NULL;
1245 Sh->IA = (INT*)fasp_mem_calloc(Sh->row+1, sizeof(INT));
1246
1247 // record the number of times some coarse point is visited
1248 visited = (INT*)fasp_mem_calloc(num_c, sizeof(INT));
1249 memset(visited, 0, sizeof(INT)*num_c);
1250
1251 /**********************************************/
1252 /* step 1: Find first the structure IA of Sh */
1253 /**********************************************/
1254
1255 Sh->IA[0] = 0;
1256
1257 for ( ci = 0; ci < Sh->row; ci++ ) {
1258
1259 i = cp_index[ci]; // find the index of the ci-th coarse grid point
1260
1261 // number of coarse point that i is strongly connected to w.r.t. S(p,2)
1262 count = 0;
1263
1264 // visit all the fine neighbors that ci is strongly connected to
1265 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1266
1267 fj = S->JA[j];
1268
1269 if ( vec[fj] == CGPT && fj != i ) {
1270 cj = cp_rindex[fj];
1271 if ( visited[cj] != ci+1 ) { // not visited yet
1272 visited[cj] = ci+1; // mark as strongly connected from ci
1273 count++;
1274 }
1275 }
1276
1277 else if ( vec[fj] == FGPT ) { // fine grid point
1278
1279 // find all the coarse neighbors that fj is strongly connected to
1280 for ( k = S->IA[fj]; k < S->IA[fj+1]; k++ ) {
1281
1282 ck = S->JA[k];
1283
1284 if ( vec[ck] == CGPT && ck != i ) { // coarse grid point
1285 if ( cp_rindex[ck] >= num_c ) {
1286 printf("### ERROR: ck=%d, num_c=%d, out of bound!\n",
1287 ck, num_c);
1288 fasp_chkerr(ERROR_AMG_COARSEING, __FUNCTION__);
1289 }
1290 cck = cp_rindex[ck];
1291
1292 if ( visited[cck] == ci+1 ) {
1293 // visited already!
1294 }
1295 else if ( visited[cck] == -ci-1 ) {
1296 visited[cck] = ci+1; // mark as strongly connected from ci
1297 count++;
1298 }
1299 else {
1300 visited[cck] = -ci-1; // mark as visited
1301 }
1302
1303 } //end if vec[ck]
1304
1305 } // end for k
1306
1307 } // end if vec[fj]
1308
1309 } // end for j
1310
1311 Sh->IA[ci+1] = Sh->IA[ci] + count;
1312
1313 } //end for i
1314
1315 /*************************/
1316 /* step 2: Find JA of Sh */
1317 /*************************/
1318
1319 memset(visited, 0, sizeof(INT)*num_c); // reset visited
1320
1321 Sh->nnz = Sh->IA[Sh->row];
1322 Sh->JA = (INT*)fasp_mem_calloc(Sh->nnz,sizeof(INT));
1323
1324 for ( ci = 0; ci < Sh->row; ci++ ) {
1325
1326 i = cp_index[ci]; // find the index of the i-th coarse grid point
1327 count = Sh->IA[ci]; // count for coarse points
1328
1329 // visit all the fine neighbors that ci is strongly connected to
1330 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1331
1332 fj = S->JA[j];
1333
1334 if ( vec[fj] == CGPT && fj != i ) {
1335 cj = cp_rindex[fj];
1336 if ( visited[cj] != ci+1 ) { // not visited yet
1337 visited[cj] = ci+1;
1338 Sh->JA[count] = cj;
1339 count++;
1340 }
1341 }
1342
1343 else if ( vec[fj] == FGPT ) { // fine grid point
1344
1345 // find all the coarse neighbors that fj is strongly connected to
1346 for ( k = S->IA[fj]; k < S->IA[fj+1]; k++ ) {
1347
1348 ck = S->JA[k];
1349
1350 if ( vec[ck] == CGPT && ck != i ) { // coarse grid point
1351 cck = cp_rindex[ck];
1352 if ( visited[cck] == ci+1 ) {
1353 // visited before
1354 }
1355 else if ( visited[cck] == -ci-1 ) {
1356 visited[cck] = ci+1;
1357 Sh->JA[count] = cck;
1358 count++;
1359 }
1360 else {
1361 visited[cck] = -ci-1;
1362 }
1363 } // end if vec[ck]
1364
1365 } // end for k
1366
1367 } // end if vec[fj]
1368
1369 } // end for j
1370
1371 if ( count != Sh->IA[ci+1] ) {
1372 printf("### WARNING: Inconsistent numbers of nonzeros!\n ");
1373 }
1374
1375 } // end for ci
1376
1377 fasp_mem_free(visited); visited = NULL;
1378}
1379
1401static INT cfsplitting_agg (dCSRmat *A,
1402 iCSRmat *S,
1403 ivector *vertices,
1404 INT aggressive_path)
1405{
1406 const INT row = A->row;
1407 INT col = 0; // initialize col(P): returning output
1408
1409 // local variables
1410 INT *vec = vertices->val, *cp_index;
1411 INT maxmeas, maxnode, num_left = 0;
1412 INT measure, newmeas;
1413 INT i, j, k, l, m, ci, cj, ck, cl, num_c;
1414 SHORT IS_CNEIGH;
1415
1416 INT *work = (INT*)fasp_mem_calloc(3*row,sizeof(INT));
1417 INT *lists = work, *where = lists+row, *lambda = where+row;
1418
1419 ivector CGPT_index, CGPT_rindex;
1420 LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
1421
1422 // Sh is for the strong coupling matrix between temporary CGPTs
1423 // ShT is the transpose of Sh
1424 // Snew is for combining the information from S and Sh
1425 iCSRmat ST, Sh, ShT;
1426
1427#if DEBUG_MODE > 0
1428 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1429#endif
1430
1431 /************************************************************/
1432 /* Coarsening Phase ONE: find temporary coarse level points */
1433 /************************************************************/
1434
1435 num_c = cfsplitting_cls(A, S, vertices);
1436 fasp_icsr_trans(S, &ST);
1437
1438 /************************************************************/
1439 /* Coarsening Phase TWO: find real coarse level points */
1440 /************************************************************/
1441
1442 // find Sh, the strong coupling between coarse grid points S(path,2)
1443 if ( aggressive_path < 2 )
1444 strong_couplings_agg1(A, S, &Sh, vertices, &CGPT_index, &CGPT_rindex);
1445 else
1446 strong_couplings_agg2(A, S, &Sh, vertices, &CGPT_index, &CGPT_rindex);
1447
1448 fasp_icsr_trans(&Sh, &ShT);
1449
1450 CGPT_index.row = num_c;
1451 CGPT_rindex.row = row;
1452 cp_index = CGPT_index.val;
1453
1454 // 1. Initialize lambda
1455#ifdef _OPENMP
1456#pragma omp parallel for if(num_c>OPENMP_HOLDS)
1457#endif
1458 for ( ci = 0; ci < num_c; ++ci ) lambda[ci] = ShT.IA[ci+1]-ShT.IA[ci];
1459
1460 // 2. Form linked list for lambda (max to min)
1461 for ( ci = 0; ci < num_c; ++ci ) {
1462
1463 i = cp_index[ci];
1464 measure = lambda[ci];
1465
1466 if ( vec[i] == ISPT ) continue; // skip isolated points
1467
1468 if ( measure > 0 ) {
1469 enter_list(&LoL_head, &LoL_tail, lambda[ci], ci, lists, where);
1470 num_left++;
1471 }
1472 else {
1473 if ( measure < 0) printf("### WARNING: Negative lambda[%d]!\n", i);
1474
1475 vec[i] = FGPT; // set i as fine node
1476
1477 // update the lambda value in the CGPT neighbor of i
1478 for ( ck = Sh.IA[ci]; ck < Sh.IA[ci+1]; ++ck ) {
1479
1480 cj = Sh.JA[ck];
1481 j = cp_index[cj];
1482
1483 if ( vec[j] == ISPT ) continue;
1484
1485 if ( cj < ci ) {
1486 newmeas = lambda[cj];
1487 if ( newmeas > 0 ) {
1488 remove_node(&LoL_head, &LoL_tail, newmeas, cj, lists, where);
1489 num_left--;
1490 }
1491 newmeas = ++(lambda[cj]);
1492 enter_list(&LoL_head, &LoL_tail, newmeas, cj, lists, where);
1493 num_left++;
1494 }
1495 else {
1496 newmeas = ++(lambda[cj]);
1497 } // end if cj<ci
1498
1499 } // end for ck
1500
1501 } // end if
1502
1503 } // end for ci
1504
1505 // 3. Main loop
1506 while ( num_left > 0 ) {
1507
1508 // pick $i\in U$ with $\max\lambda_i: C:=C\cup\{i\}, U:=U\\{i\}$
1509 maxnode = LoL_head->head;
1510 maxmeas = lambda[maxnode];
1511 if ( maxmeas == 0 ) printf("### WARNING: Head of the list has measure 0!\n");
1512
1513 // mark maxnode as real coarse node, labelled as number 3
1514 vec[cp_index[maxnode]] = 3;
1515 --num_left;
1516 remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
1517 lambda[maxnode] = 0;
1518 col++; // count for the real coarse node after aggressive coarsening
1519
1520 // for all $j\in S_i^T\cap U: F:=F\cup\{j\}, U:=U\backslash\{j\}$
1521 for ( ci = ShT.IA[maxnode]; ci < ShT.IA[maxnode+1]; ++ci ) {
1522
1523 cj = ShT.JA[ci];
1524 j = cp_index[cj];
1525
1526 if ( vec[j] != CGPT ) continue; // skip if j is not C-point
1527
1528 vec[j] = 4; // set j as 4--fake CGPT
1529 remove_node(&LoL_head, &LoL_tail, lambda[cj], cj, lists, where);
1530 --num_left;
1531
1532 // update the measure for neighboring points
1533 for ( cl = Sh.IA[cj]; cl < Sh.IA[cj+1]; cl++ ) {
1534 ck = Sh.JA[cl];
1535 k = cp_index[ck];
1536 if ( vec[k] == CGPT ) { // k is temporary CGPT
1537 remove_node(&LoL_head, &LoL_tail, lambda[ck], ck, lists, where);
1538 newmeas = ++(lambda[ck]);
1539 enter_list(&LoL_head, &LoL_tail, newmeas, ck, lists, where);
1540 }
1541 }
1542
1543 } // end for ci
1544
1545 // Update lambda and linked list after maxnode->C
1546 for ( ci = Sh.IA[maxnode]; ci < Sh.IA[maxnode+1]; ++ci ) {
1547
1548 cj = Sh.JA[ci];
1549 j = cp_index[cj];
1550
1551 if ( vec[j] != CGPT ) continue; // skip if j is not C-point
1552
1553 measure = lambda[cj];
1554 remove_node(&LoL_head, &LoL_tail, measure, cj, lists, where);
1555 lambda[cj] = --measure;
1556
1557 if ( measure > 0 ) {
1558 enter_list(&LoL_head, &LoL_tail,measure, cj, lists, where);
1559 }
1560 else {
1561 vec[j] = 4; // set j as fake CGPT variable
1562 --num_left;
1563 for ( cl = Sh.IA[cj]; cl < Sh.IA[cj+1]; cl++ ) {
1564 ck = Sh.JA[cl];
1565 k = cp_index[ck];
1566 if ( vec[k] == CGPT ) { // k is temporary CGPT
1567 remove_node(&LoL_head, &LoL_tail, lambda[ck], ck, lists, where);
1568 newmeas = ++(lambda[ck]);
1569 enter_list(&LoL_head, &LoL_tail, newmeas, ck, lists, where);
1570 }
1571 } // end for l
1572 } // end if
1573
1574 } // end for
1575
1576 } // while
1577
1578 // 4. reorganize the variable type: mark temporary CGPT--1 and fake CGPT--4 as
1579 // FGPT; mark real CGPT--3 to be CGPT
1580#ifdef _OPENMP
1581#pragma omp parallel for if(row>OPENMP_HOLDS)
1582#endif
1583 for ( i = 0; i < row; i++ ) {
1584 if ( vec[i] == CGPT || vec[i] == 4 ) vec[i] = FGPT;
1585 }
1586
1587#ifdef _OPENMP
1588#pragma omp parallel for if(row>OPENMP_HOLDS)
1589#endif
1590 for ( i = 0; i < row; i++ ) {
1591 if ( vec[i] == 3 ) vec[i] = CGPT;
1592 }
1593
1594 /************************************************************/
1595 /* Coarsening Phase THREE: all the FGPTs which have no CGPT */
1596 /* neighbors within distance 2. Change them into CGPT such */
1597 /* that the standard interpolation works! */
1598 /************************************************************/
1599
1600 for ( i = 0; i < row; i++ ) {
1601
1602 if ( vec[i] != FGPT ) continue;
1603
1604 IS_CNEIGH = FALSE; // whether there exist CGPT neighbors within distance of 2
1605
1606 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1607
1608 if ( IS_CNEIGH ) break;
1609
1610 k = S->JA[j];
1611
1612 if ( vec[k] == CGPT ) {
1613 IS_CNEIGH = TRUE;
1614 }
1615 else if ( vec[k] == FGPT ) {
1616 for ( l = S->IA[k]; l < S->IA[k+1]; l++ ) {
1617 m = S->JA[l];
1618 if ( vec[m] == CGPT ) {
1619 IS_CNEIGH = TRUE; break;
1620 }
1621 } // end for l
1622 }
1623
1624 } // end for j
1625
1626 // no CGPT neighbors in distance <= 2, mark i as CGPT
1627 if ( !IS_CNEIGH ) {
1628 vec[i] = CGPT; col++;
1629 }
1630
1631 } // end for i
1632
1633 if ( LoL_head ) {
1634 list_ptr = LoL_head;
1635 LoL_head->prev_node = NULL;
1636 LoL_head->next_node = NULL;
1637 LoL_head = list_ptr->next_node;
1638 fasp_mem_free(list_ptr); list_ptr = NULL;
1639 }
1640
1641 fasp_ivec_free(&CGPT_index);
1642 fasp_ivec_free(&CGPT_rindex);
1643 fasp_icsr_free(&Sh);
1644 fasp_icsr_free(&ST);
1645 fasp_icsr_free(&ShT);
1646 fasp_mem_free(work); work = NULL;
1647
1648#if DEBUG_MODE > 0
1649 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1650#endif
1651
1652 return col;
1653}
1654
1678static INT clean_ff_couplings (iCSRmat *S,
1679 ivector *vertices,
1680 INT row,
1681 INT col)
1682{
1683 // local variables
1684 INT *vec = vertices->val;
1685 INT *cindex = (INT *)fasp_mem_calloc(row, sizeof(INT));
1686 INT set_empty = TRUE, C_i_nonempty = FALSE;
1687 INT ci_tilde = -1, ci_tilde_mark = -1;
1688 INT ji, jj, i, j, index;
1689
1690 fasp_iarray_set(row, cindex, -1);
1691
1692 for ( i = 0; i < row; ++i ) {
1693
1694 if ( vec[i] != FGPT ) continue; // skip non F-variables
1695
1696 for ( ji = S->IA[i]; ji < S->IA[i+1]; ++ji ) {
1697 j = S->JA[ji];
1698 if ( vec[j] == CGPT ) cindex[j] = i; // mark C-neighbors
1699 else cindex[j] = -1; // reset cindex --Chensong 06/02/2013
1700 }
1701
1702 if ( ci_tilde_mark != i ) ci_tilde = -1;//???
1703
1704 for ( ji = S->IA[i]; ji < S->IA[i+1]; ++ji ) {
1705
1706 j = S->JA[ji];
1707
1708 if ( vec[j] != FGPT ) continue; // skip non F-variables
1709
1710 // check whether there is a C-connection
1711 set_empty = TRUE;
1712 for ( jj = S->IA[j]; jj < S->IA[j+1]; ++jj ) {
1713 index = S->JA[jj];
1714 if ( cindex[index] == i ) {
1715 set_empty = FALSE; break;
1716 }
1717 } // end for jj
1718
1719 // change the point i (if only F-F exists) to C
1720 if ( set_empty ) {
1721 if ( C_i_nonempty ) {
1722 vec[i] = CGPT; col++;
1723 if ( ci_tilde > -1 ) {
1724 vec[ci_tilde] = FGPT; col--;
1725 ci_tilde = -1;
1726 }
1727 C_i_nonempty = FALSE;
1728 break;
1729 }
1730 else { // temporary set j->C and roll back
1731 vec[j] = CGPT; col++;
1732 ci_tilde = j;
1733 ci_tilde_mark = i;
1734 C_i_nonempty = TRUE;
1735 i--; // roll back to check i-point again
1736 break;
1737 } // end if C_i_nonempty
1738 } // end if set_empty
1739
1740 } // end for ji
1741
1742 } // end for i
1743
1744 fasp_mem_free(cindex); cindex = NULL;
1745
1746 return col;
1747}
1748
1767static void form_P_pattern_dir (dCSRmat *P,
1768 iCSRmat *S,
1769 ivector *vertices,
1770 INT row,
1771 INT col)
1772{
1773 // local variables
1774 INT i, j, k, index;
1775 INT *vec = vertices->val;
1776
1777 SHORT nthreads = 1, use_openmp = FALSE;
1778
1779#ifdef _OPENMP
1780 if ( row > OPENMP_HOLDS ) {
1781 use_openmp = TRUE;
1782 nthreads = fasp_get_num_threads();
1783 }
1784#endif
1785
1786 // Initialize P matrix
1787 P->row = row; P->col = col;
1788 P->IA = (INT*)fasp_mem_calloc(row+1, sizeof(INT));
1789
1790 // step 1: Find the structure IA of P first: using P as a counter
1791 if ( use_openmp ) {
1792
1793 INT mybegin, myend, myid;
1794#ifdef _OPENMP
1795#pragma omp parallel for private(myid, mybegin,myend,i,j,k)
1796#endif
1797 for ( myid = 0; myid < nthreads; myid++ ) {
1798 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
1799 for ( i = mybegin; i < myend; ++i ) {
1800 switch ( vec[i] ) {
1801 case FGPT: // fine grid points
1802 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1803 k = S->JA[j];
1804 if ( vec[k] == CGPT ) P->IA[i+1]++;
1805 }
1806 break;
1807
1808 case CGPT: // coarse grid points
1809 P->IA[i+1] = 1; break;
1810
1811 default: // treat everything else as isolated
1812 P->IA[i+1] = 0; break;
1813 }
1814 }
1815 }
1816
1817 }
1818
1819 else {
1820
1821 for ( i = 0; i < row; ++i ) {
1822 switch ( vec[i] ) {
1823 case FGPT: // fine grid points
1824 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1825 k = S->JA[j];
1826 if ( vec[k] == CGPT ) P->IA[i+1]++;
1827 }
1828 break;
1829
1830 case CGPT: // coarse grid points
1831 P->IA[i+1] = 1; break;
1832
1833 default: // treat everything else as isolated
1834 P->IA[i+1] = 0; break;
1835 }
1836 } // end for i
1837
1838 } // end if
1839
1840 // Form P->IA from the counter P
1841 for ( i = 0; i < P->row; ++i ) P->IA[i+1] += P->IA[i];
1842 P->nnz = P->IA[P->row]-P->IA[0];
1843
1844 // step 2: Find the structure JA of P
1845 P->JA = (INT*)fasp_mem_calloc(P->nnz,sizeof(INT));
1846 P->val = (REAL*)fasp_mem_calloc(P->nnz,sizeof(REAL));
1847
1848 for ( index = i = 0; i < row; ++i ) {
1849 if ( vec[i] == FGPT ) { // fine grid point
1850 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1851 k = S->JA[j];
1852 if ( vec[k] == CGPT ) P->JA[index++] = k;
1853 } // end for j
1854 } // end if
1855 else if ( vec[i] == CGPT ) { // coarse grid point -- one entry only
1856 P->JA[index++] = i;
1857 }
1858 }
1859
1860}
1861
1880static void form_P_pattern_std (dCSRmat *P,
1881 iCSRmat *S,
1882 ivector *vertices,
1883 INT row,
1884 INT col)
1885{
1886 // local variables
1887 INT i, j, k, l, h, index;
1888 INT *vec = vertices->val;
1889
1890 // number of times a C-point is visited
1891 INT *visited = (INT*)fasp_mem_calloc(row,sizeof(INT));
1892
1893 P->row = row; P->col = col;
1894 P->IA = (INT*)fasp_mem_calloc(row+1, sizeof(INT));
1895
1896 fasp_iarray_set(row, visited, -1);
1897
1898 // Step 1: Find the structure IA of P first: use P as a counter
1899 for ( i = 0; i < row; ++i ) {
1900
1901 if ( vec[i] == FGPT ) { // if node i is a F point
1902 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1903
1904 k = S->JA[j];
1905
1906 // if neighbor of i is a C point, good
1907 if ( (vec[k] == CGPT) && (visited[k] != i) ) {
1908 visited[k] = i;
1909 P->IA[i+1]++;
1910 }
1911
1912 // if k is a F point and k is not i, look for indirect C neighbors
1913 else if ( (vec[k] == FGPT) && (k != i) ) {
1914 for ( l = S->IA[k]; l < S->IA[k+1]; l++ ) { // neighbors of k
1915 h = S->JA[l];
1916 if ( (vec[h] == CGPT) && (visited[h] != i) ) {
1917 visited[h] = i;
1918 P->IA[i+1]++;
1919 }
1920 } // end for(l=S->IA[k];l<S->IA[k+1];l++)
1921 } // end if (vec[k]==CGPT)
1922
1923 } // end for (j=S->IA[i];j<S->IA[i+1];j++)
1924 }
1925
1926 else if ( vec[i] == CGPT ) { // if node i is a C point
1927 P->IA[i+1] = 1;
1928 }
1929
1930 else { // treat everything else as isolated points
1931 P->IA[i+1] = 0;
1932 } // end if (vec[i]==FGPT)
1933
1934 } // end for (i=0;i<row;++i)
1935
1936 // Form P->IA from the counter P
1937 for ( i = 0; i < P->row; ++i ) P->IA[i+1] += P->IA[i];
1938 P->nnz = P->IA[P->row]-P->IA[0];
1939
1940 // Step 2: Find the structure JA of P
1941 P->JA = (INT*)fasp_mem_calloc(P->nnz,sizeof(INT));
1942 P->val = (REAL*)fasp_mem_calloc(P->nnz,sizeof(REAL));
1943
1944 fasp_iarray_set(row, visited, -1); // re-init visited array
1945
1946 for ( i = 0; i < row; ++i ) {
1947
1948 if ( vec[i] == FGPT ) { // if node i is a F point
1949
1950 index = 0;
1951
1952 for ( j = S->IA[i]; j < S->IA[i+1]; j++ ) {
1953
1954 k = S->JA[j];
1955
1956 // if neighbor k of i is a C point
1957 if ( (vec[k] == CGPT) && (visited[k] != i) ) {
1958 visited[k] = i;
1959 P->JA[P->IA[i]+index] = k;
1960 index++;
1961 }
1962
1963 // if neighbor k of i is a F point and k is not i
1964 else if ( (vec[k] == FGPT) && (k != i) ) {
1965 for ( l = S->IA[k]; l < S->IA[k+1]; l++ ) { // neighbors of k
1966 h = S->JA[l];
1967 if ( (vec[h] == CGPT) && (visited[h] != i) ) {
1968 visited[h] = i;
1969 P->JA[P->IA[i]+index] = h;
1970 index++;
1971 }
1972
1973 } // end for (l=S->IA[k];l<S->IA[k+1];l++)
1974
1975 } // end if (vec[k]==CGPT)
1976
1977 } // end for (j=S->IA[i];j<S->IA[i+1];j++)
1978 }
1979
1980 else if ( vec[i] == CGPT ) {
1981 P->JA[P->IA[i]] = i;
1982 }
1983 }
1984
1985 // clean up
1986 fasp_mem_free(visited); visited = NULL;
1987}
1988
2003static INT cfsplitting_mis (iCSRmat *S,
2004 ivector *vertices,
2005 ivector *order)
2006{
2007 const INT n = S->row;
2008
2009 INT col = 0;
2010 INT *ord = order->val;
2011 INT *vec = vertices->val;
2012 INT *IS = S->IA;
2013 INT *JS = S->JA;
2014
2015 INT i, j, ind;
2016 INT row_begin, row_end;
2017
2018 fasp_ivec_set (n, vertices, UNPT);
2019
2020 for (i=0; i<n ; i++)
2021 {
2022 ind = ord[i];
2023 if (vec[ind] == UNPT) {
2024 vec[ind] = CGPT;
2025 row_begin = IS[ind]; row_end = IS[ind+1];
2026 for (j = row_begin; j<row_end; j++)
2027 {
2028 if (vec[JS[j]] == CGPT ) {
2029 vec[ind] = FGPT;
2030 break;
2031 }
2032 }
2033 if (vec[ind] == CGPT) {
2034 col++;
2035 for (j = row_begin; j<row_end; j++)
2036 {
2037 vec[JS[j]] = FGPT;
2038 }
2039 }
2040 }
2041 }
2042 return col;
2043}
2044
2059static void ordering1 (iCSRmat *S,
2060 ivector *order)
2061{
2062 const INT n = order->row;
2063 INT * IS = S->IA;
2064 INT * ord = order->val;
2065 INT maxind, maxdeg, degree;
2066 INT i;
2067
2068 for (i = 0; i < n; i++) ord[i] = i;
2069
2070 for (maxind = maxdeg = i = 0; i < n; i++)
2071 {
2072 degree = IS[i+1] - IS[i];
2073 if (degree > maxdeg)
2074 {
2075 maxind = i;
2076 maxdeg = degree;
2077 }
2078 }
2079
2080 ord[0] = maxind;
2081 ord[maxind] = 0;
2082
2083 return;
2084}
2085
2086/*---------------------------------*/
2087/*-- End of File --*/
2088/*---------------------------------*/
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
Definition: AuxArray.c:103
void fasp_iarray_cp(const INT n, const INT *x, INT *y)
Copy an array to the other y=x.
Definition: AuxArray.c:184
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:155
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
void fasp_get_start_end(const INT procid, const INT nprocs, const INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: AuxThreads.c:92
void fasp_dvec_free(dvector *u)
Free vector data space of REAL type.
Definition: AuxVector.c:145
void fasp_ivec_free(ivector *u)
Free vector data space of INT type.
Definition: AuxVector.c:164
void fasp_ivec_alloc(const INT m, ivector *u)
Create vector data space of INT type.
Definition: AuxVector.c:125
ivector fasp_ivec_create(const INT m)
Create vector data space of INT type.
Definition: AuxVector.c:84
void fasp_ivec_set(INT n, ivector *u, const INT m)
Set ivector value to be m.
Definition: AuxVector.c:291
void fasp_icsr_trans(const iCSRmat *A, iCSRmat *AT)
Find transpose of iCSRmat matrix A.
Definition: BlaSparseCSR.c:875
void fasp_icsr_free(iCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:219
void fasp_dcsr_getdiag(INT n, const dCSRmat *A, dvector *diag)
Get first n diagonal entries of a CSR matrix A.
Definition: BlaSparseCSR.c:537
INT fasp_amg_coarsening_cr(const INT i_0, const INT i_n, dCSRmat *A, ivector *vertices, AMG_param *param)
CR coarsening.
SHORT fasp_amg_coarsening_rs(dCSRmat *A, ivector *vertices, dCSRmat *P, iCSRmat *S, AMG_param *param)
Standard and aggressive coarsening schemes.
Main header file for the FASP project.
#define MIN(a, b)
Definition: fasp.h:75
#define REAL
Definition: fasp.h:67
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define ABS(a)
Definition: fasp.h:76
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:74
#define INT
Definition: fasp.h:64
#define ISPT
Definition: fasp_const.h:228
#define COARSE_MIS
Definition: fasp_const.h:211
#define ERROR_AMG_INTERP_TYPE
Definition: fasp_const.h:35
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
#define FGPT
Definition: fasp_const.h:226
#define COARSE_RSP
Definition: fasp_const.h:208
#define INTERP_EXT
Definition: fasp_const.h:219
#define CGPT
Definition: fasp_const.h:227
#define COARSE_AC
Definition: fasp_const.h:210
#define UNPT
Definition: fasp_const.h:225
#define INTERP_STD
Definition: fasp_const.h:217
#define INTERP_DIR
Definition of interpolation types.
Definition: fasp_const.h:216
#define OPENMP_HOLDS
Definition: fasp_const.h:260
#define TRUE
Definition of logic type.
Definition: fasp_const.h:61
#define FALSE
Definition: fasp_const.h:62
#define ERROR_AMG_COARSEING
Definition: fasp_const.h:38
#define COARSE_CR
Definition: fasp_const.h:209
#define ERROR_UNKNOWN
Definition: fasp_const.h:56
#define INTERP_ENG
Definition: fasp_const.h:218
Parameters for AMG methods.
Definition: fasp.h:447
REAL strong_threshold
strong connection threshold for coarsening
Definition: fasp.h:516
INT aggressive_path
number of paths use to determine strongly coupled C points
Definition: fasp.h:528
SHORT interpolation_type
interpolation type
Definition: fasp.h:513
SHORT coarsening_type
coarsening type
Definition: fasp.h:507
REAL max_row_sum
maximal row sum parameter
Definition: fasp.h:519
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:143
INT col
column of matrix A, n
Definition: fasp.h:149
REAL * val
nonzero entries of A
Definition: fasp.h:161
INT row
row number of matrix A, m
Definition: fasp.h:146
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:155
INT nnz
number of nonzero entries
Definition: fasp.h:152
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:158
Vector with n entries of REAL type.
Definition: fasp.h:346
REAL * val
actual vector entries
Definition: fasp.h:352
Sparse matrix of INT type in CSR format.
Definition: fasp.h:182
INT col
column of matrix A, n
Definition: fasp.h:188
INT row
row number of matrix A, m
Definition: fasp.h:185
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:194
INT nnz
number of nonzero entries
Definition: fasp.h:191
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:197
INT * val
nonzero entries of A
Definition: fasp.h:200
Vector with n entries of INT type.
Definition: fasp.h:361
INT row
number of rows
Definition: fasp.h:364
INT * val
actual vector entries
Definition: fasp.h:367