Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
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#include <math.h>
33
34static INT cfsplitting_cls(dCSRmat*, iCSRmat*, ivector*);
35static INT cfsplitting_clsp(dCSRmat*, iCSRmat*, ivector*);
36static INT cfsplitting_agg(dCSRmat*, iCSRmat*, ivector*, INT);
37static INT cfsplitting_mis(iCSRmat*, ivector*, ivector*);
38static INT clean_ff_couplings(iCSRmat*, ivector*, INT, INT);
39static INT compress_S(iCSRmat*);
40
41static void strong_couplings(dCSRmat*, iCSRmat*, AMG_param*);
42static void form_P_pattern_dir(dCSRmat*, iCSRmat*, ivector*, INT, INT);
43static void form_P_pattern_std(dCSRmat*, iCSRmat*, ivector*, INT, INT);
44static void ordering1(iCSRmat*, ivector*);
45
46static void form_P_pattern_rdc(dCSRmat*, dCSRmat*, double*, ivector*, INT, INT);
47
48/*---------------------------------*/
49/*-- Public Functions --*/
50/*---------------------------------*/
51
76SHORT fasp_amg_coarsening_rs(
77 dCSRmat* A, ivector* vertices, dCSRmat* P, iCSRmat* S, 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 // printf("### DEBUG: Step 2. C/F splitting ......\n");
105 // printf("coarse_type:%d \n", coarse_type);
106 switch (coarse_type) {
107
108 case COARSE_RSP: // Classical coarsening with positive connections
109 col = cfsplitting_clsp(A, S, vertices);
110 break;
111
112 case COARSE_AC: // Aggressive coarsening
113 col = cfsplitting_agg(A, S, vertices, agg_path);
114 break;
115
116 case COARSE_CR: // Compatible relaxation
117 col = fasp_amg_coarsening_cr(0, row - 1, A, vertices, param);
118 break;
119
120 case COARSE_MIS: // Maximal independent set
121 {
122 ivector order = fasp_ivec_create(row);
123 compress_S(S);
124 ordering1(S, &order);
125 col = cfsplitting_mis(S, vertices, &order);
126 fasp_ivec_free(&order);
127 break;
128 }
129
130 default: // Classical coarsening
131 col = cfsplitting_cls(A, S, vertices);
132 }
133
134#if DEBUG_MODE > 1
135 printf("### DEBUG: col = %d\n", col);
136#endif
137 // printf("### DEBUG: col = %d\n", col);
138 if (col <= 0) return ERROR_UNKNOWN;
139
140#if DEBUG_MODE > 1
141 printf("### DEBUG: Step 3. Find support of C points ......\n");
142#endif
143
144 switch (interp_type) {
145
146 case INTERP_DIR: // Direct interpolation or ...
147 case INTERP_ENG: // Energy-min interpolation
148 col = clean_ff_couplings(S, vertices, row, col);
149 form_P_pattern_dir(P, S, vertices, row, col);
150 break;
151
152 case INTERP_STD: // Standard interpolation
153 case INTERP_EXT: // Extended interpolation
154 form_P_pattern_std(P, S, vertices, row, col);
155 break;
156
157 case INTERP_RDC: // Reduction-based amg interpolation
158 {
159 // printf("### DEBUG: Reduction-based interpolation\n");
160 INT i;
161 double* theta = (double*)fasp_mem_calloc(row, sizeof(double));
162 form_P_pattern_rdc(P, A, theta, vertices, row, col);
163 // theta will be used to
164 // 1. compute entries of interpolation matrix
165 // 2. calculate relaxation parameter
166 // 3. approximate convergence factor
167 param->theta = 1.0;
168 for (i = 0; i < row; ++i)
169 if (theta[i] < param->theta) param->theta = theta[i];
170 printf("### DEBUG: theta = %e\n", param->theta);
171 fasp_mem_free(theta);
172 theta = NULL;
173 double eps = (2 - 2 * param->theta) / (2 * param->theta - 1);
174 // assume: two-grid
175 double nu = (param->presmooth_iter + param->postsmooth_iter) / 2.0;
176 double conv_factor = (eps / (1 + eps)) *
177 (1 + pow(eps, 2 * nu - 1) / pow(2 + eps, 2 * nu));
178 conv_factor = sqrt(conv_factor);
179 printf("### DEBUG: Upper bound for conv_factor = %e\n", conv_factor);
180 if (param->theta <= 0.5) {
181 REAL reset_value = 0.5 + 1e-5;
182 printf("### WARNING: theta = %e <= 0.5, use %e instead \n",
183 param->theta, reset_value);
184 param->theta = reset_value;
185 }
186 if (param->theta >= 0.0) {
187 // weighted Jacobi smoother only on F-points
188 REAL theta = param->theta;
189 REAL eps = (2 - 2 * theta) / (2 * theta - 1);
190 REAL sigma = 2 / (2 + eps);
191 REAL weight = sigma / (2 - 1 / theta);
192 printf("### DEBUG: theta = %e, eps = %e, sigma = %e\n", theta, eps,
193 sigma);
194
195 // set reduction-based AMG smoother parameters
196 if (param->smoother == SMOOTHER_JACOBIF) {
197 param->relaxation = weight;
198 printf("### DEBUG: Weight for JACOBI_F = %e\n", weight);
199 }
200 }
201 break;
202 }
203
204 default:
205 fasp_chkerr(ERROR_AMG_INTERP_TYPE, __FUNCTION__);
206 }
207
208#if DEBUG_MODE > 0
209 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
210#endif
211
212 return FASP_SUCCESS;
213}
214
215/*---------------------------------*/
216/*-- Private Functions --*/
217/*---------------------------------*/
218
236static void strong_couplings(dCSRmat* A, iCSRmat* S, AMG_param* param)
237{
238 const SHORT coarse_type = param->coarsening_type;
239 const REAL max_row_sum = param->max_row_sum;
240 const REAL epsilon_str = param->strong_threshold;
241 const INT row = A->row, col = A->col, row1 = row + 1;
242 const INT nnz = A->nnz;
243
244 INT * ia = A->IA, *ja = A->JA;
245 REAL* aj = A->val;
246
247 // local variables
248 INT i, j, begin_row, end_row;
249 REAL row_scl, row_sum;
250
251 SHORT nthreads = 1, use_openmp = FALSE;
252
253#ifdef _OPENMP
254 if (row > OPENMP_HOLDS) {
255 use_openmp = TRUE;
256 nthreads = fasp_get_num_threads();
257 }
258#endif
259
260 // get the diagonal entry of A: assume all connections are strong
261 dvector diag;
262 fasp_dcsr_getdiag(0, A, &diag);
263
264 // copy the structure of A to S
265 S->row = row;
266 S->col = col;
267 S->nnz = nnz;
268 S->val = NULL;
269 S->IA = (INT*)fasp_mem_calloc(row1, sizeof(INT));
270 S->JA = (INT*)fasp_mem_calloc(nnz, sizeof(INT));
271 fasp_iarray_cp(row1, ia, S->IA);
272 fasp_iarray_cp(nnz, ja, S->JA);
273
274 if (use_openmp) {
275
276 // This part is still old! Need to be updated. --Chensong 09/18/2016
277
278 INT mybegin, myend, myid;
279#ifdef _OPENMP
280#pragma omp parallel for private(myid, mybegin, myend, i, row_scl, row_sum, begin_row, \
281 end_row, j)
282#endif
283 for (myid = 0; myid < nthreads; myid++) {
284 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
285 for (i = mybegin; i < myend; i++) {
286
287 // Compute most negative entry in each row and row sum
288 row_scl = row_sum = 0.0;
289 begin_row = ia[i];
290 end_row = ia[i + 1];
291 for (j = begin_row; j < end_row; j++) {
292 row_scl = MIN(row_scl, aj[j]);
293 row_sum += aj[j];
294 }
295
296 // Find diagonal entries of S and remove them later
297 for (j = begin_row; j < end_row; j++) {
298 if (ja[j] == i) {
299 S->JA[j] = -1;
300 break;
301 }
302 }
303
304 // Mark entire row as weak couplings if strongly diagonal-dominant
305 if (ABS(row_sum) > max_row_sum * ABS(diag.val[i])) {
306 for (j = begin_row; j < end_row; j++) S->JA[j] = -1;
307 } else {
308 for (j = begin_row; j < end_row; j++) {
309 // If a_{ij} >= \epsilon_{str} * \min a_{ij}, the connection
310 // j->i is set to be weak; positive entries result in weak
311 // connections
312 if (A->val[j] >= epsilon_str * row_scl) S->JA[j] = -1;
313 }
314 }
315
316 } // end for i
317 } // end for myid
318
319 }
320
321 else {
322
323 for (i = 0; i < row; ++i) {
324
325 // Compute row scale and row sum
326 row_scl = row_sum = 0.0;
327 begin_row = ia[i];
328 end_row = ia[i + 1];
329
330 for (j = begin_row; j < end_row; j++) {
331
332 // Originally: Not consider positive entries
333 // row_sum += aj[j];
334 // Now changed to --Chensong 05/17/2013
335 row_sum += ABS(aj[j]);
336
337 // Originally: Not consider positive entries
338 // row_scl = MAX(row_scl, -aj[j]); // smallest negative
339 // Now changed to --Chensong 06/01/2013
340 if (ja[j] != i) row_scl = MAX(row_scl, ABS(aj[j])); // largest abs
341 }
342
343 // Multiply by the strength threshold
344 row_scl *= epsilon_str;
345 // printf("row_sum:%e, row_scl:%e \n", row_sum, row_scl);
346 // Find diagonal entries of S and remove them later
347 for (j = begin_row; j < end_row; j++) {
348 if (ja[j] == i) {
349 S->JA[j] = -1;
350 break;
351 }
352 }
353
354 // Mark entire row as weak couplings if strongly diagonal-dominant
355 // Originally: Not consider positive entries
356 // if ( ABS(row_sum) > max_row_sum * ABS(diag.val[i]) ) {
357 // Now changed to --Chensong 05/17/2013
358 /* printf(" Mark entire row:%e, \n", (2 - max_row_sum) *
359 * ABS(diag.val[i]));*/
360 if (row_sum < (2 - max_row_sum) * ABS(diag.val[i])) {
361
362 for (j = begin_row; j < end_row; j++) S->JA[j] = -1;
363
364 } else {
365
366 switch (coarse_type) {
367
368 case COARSE_RSP: // consider positive off-diag as well
369 for (j = begin_row; j < end_row; j++) {
370 if (ABS(A->val[j]) <= row_scl) S->JA[j] = -1;
371 }
372 break;
373
374 default: // only consider n-couplings
375 for (j = begin_row; j < end_row; j++) {
376 if (-A->val[j] <= row_scl) S->JA[j] = -1;
377 // printf("j : %d , val[j] : %e", j, A->val[j]);
378 }
379 break;
380 }
381 }
382 } // end for i
383
384 } // end if openmp
385
386 fasp_dvec_free(&diag);
387}
388
403static INT compress_S(iCSRmat* S)
404{
405 const INT row = S->row;
406 INT* ia = S->IA;
407
408 // local variables
409 INT index, i, j, begin_row, end_row;
410
411 // compress S: remove weak connections and form strong coupling matrix
412 for (index = i = 0; i < row; ++i) {
413
414 begin_row = ia[i];
415 end_row = ia[i + 1];
416
417 ia[i] = index;
418 for (j = begin_row; j < end_row; j++) {
419 if (S->JA[j] > -1) S->JA[index++] = S->JA[j]; // strong couplings
420 }
421 }
422
423 S->nnz = S->IA[row] = index;
424
425 if (S->nnz <= 0) {
426 return ERROR_UNKNOWN;
427 } else {
428 return FASP_SUCCESS;
429 }
430}
431
444static void rem_positive_ff(dCSRmat* A, iCSRmat* Stemp, ivector* vertices)
445{
446 const INT row = A->row;
447 INT * ia = A->IA, *vec = vertices->val;
448
449 REAL row_scl, max_entry;
450 INT i, j, ji, max_index;
451
452 for (i = 0; i < row; ++i) {
453
454 if (vec[i] != FGPT) continue; // skip non F-variables
455
456 row_scl = 0.0;
457 for (ji = ia[i]; ji < ia[i + 1]; ++ji) {
458 j = A->JA[ji];
459 if (j == i) continue; // skip diagonal
460 row_scl = MAX(row_scl, ABS(A->val[ji])); // max abs entry
461 } // end for ji
462 row_scl *= 0.75;
463
464 // looking for strong F-F connections
465 max_index = -1;
466 max_entry = 0.0;
467 for (ji = ia[i]; ji < ia[i + 1]; ++ji) {
468 j = A->JA[ji];
469 if (j == i) continue; // skip diagonal
470 if (vec[j] != FGPT) continue; // skip F-C connections
471 if (A->val[ji] > row_scl) {
472 Stemp->JA[ji] = j;
473 if (A->val[ji] > max_entry) {
474 max_entry = A->val[ji];
475 max_index = j; // max positive entry
476 }
477 }
478 } // end for ji
479
480 // mark max positive entry as C-point
481 if (max_index != -1) vec[max_index] = CGPT;
482
483 } // end for i
484}
485
507static INT cfsplitting_cls(dCSRmat* A, iCSRmat* S, ivector* vertices)
508{
509 const INT row = A->row;
510
511 // local variables
512 INT col = 0;
513 INT maxmeas, maxnode, num_left = 0;
514 INT measure, newmeas;
515 INT i, j, k, l;
516 INT myid, mybegin, myend;
517 INT* vec = vertices->val;
518 INT* work = (INT*)fasp_mem_calloc(3 * row, sizeof(INT));
519 INT *lists = work, *where = lists + row, *lambda = where + row;
520
521#if RS_C1
522 INT set_empty = 1;
523 INT jkeep = 0, cnt, index;
524 INT row_end_S, ji, row_end_S_nabor, jj;
525 INT* graph_array = lambda;
526#else
527 INT* ia = A->IA;
528#endif
529
530 LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
531
532 SHORT nthreads = 1, use_openmp = FALSE;
533
534#if DEBUG_MODE > 0
535 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
536#endif
537
538#ifdef _OPENMP
539 if (row > OPENMP_HOLDS) {
540 use_openmp = TRUE;
541 nthreads = fasp_get_num_threads();
542 }
543#endif
544
545 // 0. Compress S and form S_transpose
546
547 col = compress_S(S);
548
549 if (col < 0) goto FINISHED; // compression failed!!!
550
551 iCSRmat ST;
552 fasp_icsr_trans(S, &ST);
553
554 // 1. Initialize lambda
555 if (use_openmp) {
556#ifdef _OPENMP
557#pragma omp parallel for private(myid, mybegin, myend, i)
558#endif
559 for (myid = 0; myid < nthreads; myid++) {
560 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
561 for (i = mybegin; i < myend; i++) lambda[i] = ST.IA[i + 1] - ST.IA[i];
562 }
563 } else {
564 for (i = 0; i < row; ++i) lambda[i] = ST.IA[i + 1] - ST.IA[i];
565 }
566
567 // 2. Before C/F splitting algorithm starts, filter out the variables which
568 // have no connections at all and mark them as special F-variables.
569 if (use_openmp) {
570
571#ifdef _OPENMP
572#pragma omp parallel for reduction(+ : num_left) private(myid, mybegin, myend, i)
573#endif
574 for (myid = 0; myid < nthreads; myid++) {
575 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
576 for (i = mybegin; i < myend; i++) {
577#if RS_C1 // Check C1 criteria or not
578 if (S->IA[i + 1] == S->IA[i])
579#else
580 if ((ia[i + 1] - ia[i]) <= 1)
581#endif
582 {
583 vec[i] = ISPT; // set i as an ISOLATED fine node
584 lambda[i] = 0;
585 } else {
586 vec[i] = UNPT; // set i as a undecided node
587 num_left++;
588 }
589 }
590 } // end for myid
591
592 }
593
594 else {
595
596 for (i = 0; i < row; ++i) {
597
598#if RS_C1
599 if (S->IA[i + 1] == S->IA[i])
600#else
601 if ((ia[i + 1] - ia[i]) <= 1)
602#endif
603 {
604 vec[i] = ISPT; // set i as an ISOLATED fine node
605 lambda[i] = 0;
606 } else {
607 vec[i] = UNPT; // set i as a undecided node
608 num_left++;
609 }
610 } // end for i
611 }
612
613 // 3. Form linked list for lambda (max to min)
614 for (i = 0; i < row; ++i) {
615
616 if (vec[i] == ISPT) continue; // skip isolated variables
617
618 measure = lambda[i];
619
620 if (measure > 0) {
621 enter_list(&LoL_head, &LoL_tail, lambda[i], i, lists, where);
622 } else {
623
624 if (measure < 0) printf("### WARNING: Negative lambda[%d]!\n", i);
625
626 // Set variables with non-positive measure as F-variables
627 vec[i] = FGPT; // no strong connections, set i as fine node
628 --num_left;
629
630 // Update lambda and linked list after i->F
631 for (k = S->IA[i]; k < S->IA[i + 1]; ++k) {
632 j = S->JA[k];
633 if (vec[j] == ISPT) continue; // skip isolate variables
634 if (j < i) {
635 newmeas = lambda[j];
636 if (newmeas > 0) {
637 remove_node(&LoL_head, &LoL_tail, newmeas, j, lists, where);
638 }
639 newmeas = ++(lambda[j]);
640 enter_list(&LoL_head, &LoL_tail, newmeas, j, lists, where);
641 } else {
642 newmeas = ++(lambda[j]);
643 }
644 }
645
646 } // end if measure
647
648 } // end for i
649
650 // 4. Main loop
651 while (num_left > 0) {
652
653 // pick $i\in U$ with $\max\lambda_i: C:=C\cup\{i\}, U:=U\\{i\}$
654 maxnode = LoL_head->head;
655 maxmeas = lambda[maxnode];
656 if (maxmeas == 0) printf("### WARNING: Head of the list has measure 0!\n");
657
658 vec[maxnode] = CGPT; // set maxnode as coarse node
659 lambda[maxnode] = 0;
660 --num_left;
661 remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
662 col++;
663
664 // for all $j\in S_i^T\cap U: F:=F\cup\{j\}, U:=U\backslash\{j\}$
665 for (i = ST.IA[maxnode]; i < ST.IA[maxnode + 1]; ++i) {
666
667 j = ST.JA[i];
668
669 if (vec[j] != UNPT) continue; // skip decided variables
670
671 vec[j] = FGPT; // set j as fine node
672 remove_node(&LoL_head, &LoL_tail, lambda[j], j, lists, where);
673 --num_left;
674
675 // Update lambda and linked list after j->F
676 for (l = S->IA[j]; l < S->IA[j + 1]; l++) {
677 k = S->JA[l];
678 if (vec[k] == UNPT) { // k is unknown
679 remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
680 newmeas = ++(lambda[k]);
681 enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
682 }
683 }
684
685 } // end for i
686
687 // Update lambda and linked list after maxnode->C
688 for (i = S->IA[maxnode]; i < S->IA[maxnode + 1]; ++i) {
689
690 j = S->JA[i];
691
692 if (vec[j] != UNPT) continue; // skip decided variables
693
694 measure = lambda[j];
695 remove_node(&LoL_head, &LoL_tail, measure, j, lists, where);
696 lambda[j] = --measure;
697
698 if (measure > 0) {
699 enter_list(&LoL_head, &LoL_tail, measure, j, lists, where);
700 } else { // j is the only point left, set as fine variable
701 vec[j] = FGPT;
702 --num_left;
703
704 // Update lambda and linked list after j->F
705 for (l = S->IA[j]; l < S->IA[j + 1]; l++) {
706 k = S->JA[l];
707 if (vec[k] == UNPT) { // k is unknown
708 remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
709 newmeas = ++(lambda[k]);
710 enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
711 }
712 } // end for l
713 } // end if
714
715 } // end for
716
717 } // end while
718
719#if RS_C1
720
721 // C/F splitting of RS coarsening check C1 Criterion
722 fasp_iarray_set(row, graph_array, -1);
723 for (i = 0; i < row; i++) {
724 if (vec[i] == FGPT) {
725 row_end_S = S->IA[i + 1];
726 for (ji = S->IA[i]; ji < row_end_S; ji++) {
727 j = S->JA[ji];
728 if (vec[j] == CGPT) {
729 graph_array[j] = i;
730 }
731 }
732 cnt = 0;
733 for (ji = S->IA[i]; ji < row_end_S; ji++) {
734 j = S->JA[ji];
735 if (vec[j] == FGPT) {
736 set_empty = 1;
737 row_end_S_nabor = S->IA[j + 1];
738 for (jj = S->IA[j]; jj < row_end_S_nabor; jj++) {
739 index = S->JA[jj];
740 if (graph_array[index] == i) {
741 set_empty = 0;
742 break;
743 }
744 }
745 if (set_empty) {
746 if (cnt == 0) {
747 vec[j] = CGPT;
748 col++;
749 graph_array[j] = i;
750 jkeep = j;
751 cnt = 1;
752 } else {
753 vec[i] = CGPT;
754 vec[jkeep] = FGPT;
755 break;
756 }
757 }
758 }
759 }
760 }
761 }
762
763#endif
764
765 fasp_icsr_free(&ST);
766
767 if (LoL_head) {
768 list_ptr = LoL_head;
769 LoL_head->prev_node = NULL;
770 LoL_head->next_node = NULL;
771 LoL_head = list_ptr->next_node;
772 fasp_mem_free(list_ptr);
773 list_ptr = NULL;
774 }
775
776FINISHED:
777 fasp_mem_free(work);
778 work = NULL;
779
780#if DEBUG_MODE > 0
781 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
782#endif
783
784 return col;
785}
786
806static INT cfsplitting_clsp(dCSRmat* A, iCSRmat* S, ivector* vertices)
807{
808 const INT row = A->row;
809
810 // local variables
811 INT col = 0;
812 INT maxmeas, maxnode, num_left = 0;
813 INT measure, newmeas;
814 INT i, j, k, l;
815 INT myid, mybegin, myend;
816
817 INT *ia = A->IA, *vec = vertices->val;
818 INT* work = (INT*)fasp_mem_calloc(3 * row, sizeof(INT));
819 INT *lists = work, *where = lists + row, *lambda = where + row;
820
821 LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
822
823 SHORT nthreads = 1, use_openmp = FALSE;
824
825#if DEBUG_MODE > 0
826 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
827#endif
828
829#ifdef _OPENMP
830 if (row > OPENMP_HOLDS) {
831 use_openmp = TRUE;
832 nthreads = fasp_get_num_threads();
833 }
834#endif
835
836 // 0. Compress S and form S_transpose (not complete, just IA and JA)
837 iCSRmat Stemp;
838 Stemp.row = S->row;
839 Stemp.col = S->col;
840 Stemp.nnz = S->nnz;
841 Stemp.IA = (INT*)fasp_mem_calloc(S->row + 1, sizeof(INT));
842 fasp_iarray_cp(S->row + 1, S->IA, Stemp.IA);
843 Stemp.JA = (INT*)fasp_mem_calloc(S->nnz, sizeof(INT));
844 fasp_iarray_cp(S->nnz, S->JA, Stemp.JA);
845
846 if (compress_S(S) < 0) goto FINISHED; // compression failed!!!
847
848 iCSRmat ST;
849 fasp_icsr_trans(S, &ST);
850
851 // 1. Initialize lambda
852 if (use_openmp) {
853#ifdef _OPENMP
854#pragma omp parallel for private(myid, mybegin, myend, i)
855#endif
856 for (myid = 0; myid < nthreads; myid++) {
857 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
858 for (i = mybegin; i < myend; i++) lambda[i] = ST.IA[i + 1] - ST.IA[i];
859 }
860 } else {
861 for (i = 0; i < row; ++i) lambda[i] = ST.IA[i + 1] - ST.IA[i];
862 }
863
864 // 2. Before C/F splitting algorithm starts, filter out the variables which
865 // have no connections at all and mark them as special F-variables.
866 if (use_openmp) {
867
868#ifdef _OPENMP
869#pragma omp parallel for reduction(+ : num_left) private(myid, mybegin, myend, i)
870#endif
871 for (myid = 0; myid < nthreads; myid++) {
872 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
873 for (i = mybegin; i < myend; i++) {
874 if ((ia[i + 1] - ia[i]) <= 1) {
875 vec[i] = ISPT; // set i as an ISOLATED fine node
876 lambda[i] = 0;
877 } else {
878 vec[i] = UNPT; // set i as a undecided node
879 num_left++;
880 }
881 }
882 } // end for myid
883
884 } else {
885
886 for (i = 0; i < row; ++i) {
887 if ((ia[i + 1] - ia[i]) <= 1) {
888 vec[i] = ISPT; // set i as an ISOLATED fine node
889 lambda[i] = 0;
890 } else {
891 vec[i] = UNPT; // set i as a undecided node
892 num_left++;
893 }
894 } // end for i
895 }
896
897 // 3. Form linked list for lambda (max to min)
898 for (i = 0; i < row; ++i) {
899
900 if (vec[i] == ISPT) continue; // skip isolated variables
901
902 measure = lambda[i];
903
904 if (measure > 0) {
905 enter_list(&LoL_head, &LoL_tail, lambda[i], i, lists, where);
906 } else {
907
908 if (measure < 0) printf("### WARNING: Negative lambda[%d]!\n", i);
909
910 // Set variables with non-positive measure as F-variables
911 vec[i] = FGPT; // no strong connections, set i as fine node
912 --num_left;
913
914 // Update lambda and linked list after i->F
915 for (k = S->IA[i]; k < S->IA[i + 1]; ++k) {
916
917 j = S->JA[k];
918 if (vec[j] == ISPT) continue; // skip isolate variables
919
920 if (j < i) { // only look at the previous points!!
921 newmeas = lambda[j];
922 if (newmeas > 0) {
923 remove_node(&LoL_head, &LoL_tail, newmeas, j, lists, where);
924 }
925 newmeas = ++(lambda[j]);
926 enter_list(&LoL_head, &LoL_tail, newmeas, j, lists, where);
927 } else { // will be checked later on
928 newmeas = ++(lambda[j]);
929 } // end if
930
931 } // end for k
932
933 } // end if measure
934
935 } // end for i
936
937 // 4. Main loop
938 while (num_left > 0) {
939
940 // pick $i\in U$ with $\max\lambda_i: C:=C\cup\{i\}, U:=U\\{i\}$
941 maxnode = LoL_head->head;
942 maxmeas = lambda[maxnode];
943 if (maxmeas == 0) printf("### WARNING: Head of the list has measure 0!\n");
944
945 vec[maxnode] = CGPT; // set maxnode as coarse node
946 lambda[maxnode] = 0;
947 --num_left;
948 remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
949 col++;
950
951 // for all $j\in S_i^T\cap U: F:=F\cup\{j\}, U:=U\backslash\{j\}$
952 for (i = ST.IA[maxnode]; i < ST.IA[maxnode + 1]; ++i) {
953
954 j = ST.JA[i];
955
956 if (vec[j] != UNPT) continue; // skip decided variables
957
958 vec[j] = FGPT; // set j as fine node
959 remove_node(&LoL_head, &LoL_tail, lambda[j], j, lists, where);
960 --num_left;
961
962 // Update lambda and linked list after j->F
963 for (l = S->IA[j]; l < S->IA[j + 1]; l++) {
964 k = S->JA[l];
965 if (vec[k] == UNPT) { // k is unknown
966 remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
967 newmeas = ++(lambda[k]);
968 enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
969 }
970 }
971
972 } // end for i
973
974 // Update lambda and linked list after maxnode->C
975 for (i = S->IA[maxnode]; i < S->IA[maxnode + 1]; ++i) {
976
977 j = S->JA[i];
978
979 if (vec[j] != UNPT) continue; // skip decided variables
980
981 measure = lambda[j];
982 remove_node(&LoL_head, &LoL_tail, measure, j, lists, where);
983 lambda[j] = --measure;
984
985 if (measure > 0) {
986 enter_list(&LoL_head, &LoL_tail, measure, j, lists, where);
987 } else { // j is the only point left, set as fine variable
988 vec[j] = FGPT;
989 --num_left;
990
991 // Update lambda and linked list after j->F
992 for (l = S->IA[j]; l < S->IA[j + 1]; l++) {
993 k = S->JA[l];
994 if (vec[k] == UNPT) { // k is unknown
995 remove_node(&LoL_head, &LoL_tail, lambda[k], k, lists, where);
996 newmeas = ++(lambda[k]);
997 enter_list(&LoL_head, &LoL_tail, newmeas, k, lists, where);
998 }
999 } // end for l
1000 } // end if
1001
1002 } // end for
1003
1004 } // end while
1005
1006 fasp_icsr_free(&ST);
1007
1008 if (LoL_head) {
1009 list_ptr = LoL_head;
1010 LoL_head->prev_node = NULL;
1011 LoL_head->next_node = NULL;
1012 LoL_head = list_ptr->next_node;
1013 fasp_mem_free(list_ptr);
1014 list_ptr = NULL;
1015 }
1016
1017 // Enforce F-C connections. Adding this step helps for the ExxonMobil test
1018 // problems! Need more tests though --Chensong 06/08/2013
1019 // col = clean_ff_couplings(S, vertices, row, col);
1020
1021 rem_positive_ff(A, &Stemp, vertices);
1022
1023 if (compress_S(&Stemp) < 0) goto FINISHED; // compression failed!!!
1024
1025 S->row = Stemp.row;
1026 S->col = Stemp.col;
1027 S->nnz = Stemp.nnz;
1028
1029 fasp_mem_free(S->IA);
1030 S->IA = Stemp.IA;
1031 fasp_mem_free(S->JA);
1032 S->JA = Stemp.JA;
1033
1034FINISHED:
1035 fasp_mem_free(work);
1036 work = NULL;
1037
1038#if DEBUG_MODE > 0
1039 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1040#endif
1041
1042 return col;
1043}
1044
1065static void strong_couplings_agg1(dCSRmat* A,
1066 iCSRmat* S,
1067 iCSRmat* Sh,
1068 ivector* vertices,
1069 ivector* CGPT_index,
1070 ivector* CGPT_rindex)
1071{
1072 const INT row = A->row;
1073
1074 // local variables
1075 INT i, j, k;
1076 INT num_c, count, ci, cj, ck, fj, cck;
1077 INT *cp_index, *cp_rindex, *visited;
1078 INT* vec = vertices->val;
1079
1080 // count the number of coarse grid points
1081 for (num_c = i = 0; i < row; i++) {
1082 if (vec[i] == CGPT) num_c++;
1083 }
1084
1085 // for the reverse indexing of coarse grid points
1086 fasp_ivec_alloc(row, CGPT_rindex);
1087 cp_rindex = CGPT_rindex->val;
1088
1089 // generate coarse grid point index
1090 fasp_ivec_alloc(num_c, CGPT_index);
1091 cp_index = CGPT_index->val;
1092 for (j = i = 0; i < row; i++) {
1093 if (vec[i] == CGPT) {
1094 cp_index[j] = i;
1095 cp_rindex[i] = j;
1096 j++;
1097 }
1098 }
1099
1100 // allocate space for Sh
1101 Sh->row = Sh->col = num_c;
1102 Sh->val = Sh->JA = NULL;
1103 Sh->IA = (INT*)fasp_mem_calloc(Sh->row + 1, sizeof(INT));
1104
1105 // record the number of times some coarse point is visited
1106 visited = (INT*)fasp_mem_calloc(num_c, sizeof(INT));
1107 fasp_iarray_set(num_c, visited, -1);
1108
1109 /**********************************************/
1110 /* step 1: Find first the structure IA of Sh */
1111 /**********************************************/
1112
1113 Sh->IA[0] = 0;
1114
1115 for (ci = 0; ci < Sh->row; ci++) {
1116
1117 i = cp_index[ci]; // find the index of the ci-th coarse grid point
1118
1119 // number of coarse point that i is strongly connected to w.r.t. S(p,2)
1120 count = 0;
1121
1122 // visit all the fine neighbors that ci is strongly connected to
1123 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
1124
1125 fj = S->JA[j];
1126
1127 if (vec[fj] == CGPT && fj != i) {
1128 cj = cp_rindex[fj];
1129 if (visited[cj] != ci) {
1130 visited[cj] = ci; // mark as strongly connected from ci
1131 count++;
1132 }
1133
1134 }
1135
1136 else if (vec[fj] == FGPT) { // fine grid point,
1137
1138 // find all the coarse neighbors that fj is strongly connected to
1139 for (k = S->IA[fj]; k < S->IA[fj + 1]; k++) {
1140 ck = S->JA[k];
1141 if (vec[ck] == CGPT && ck != i) { // it is a coarse grid point
1142 if (cp_rindex[ck] >= num_c) {
1143 printf("### ERROR: ck=%d, num_c=%d, out of bound!\n", ck,
1144 num_c);
1145 fasp_chkerr(ERROR_AMG_COARSEING, __FUNCTION__);
1146 }
1147 cck = cp_rindex[ck];
1148
1149 if (visited[cck] != ci) {
1150 visited[cck] = ci; // mark as strongly connected from ci
1151 count++;
1152 }
1153 } // end if
1154 } // end for k
1155
1156 } // end if
1157
1158 } // end for j
1159
1160 Sh->IA[ci + 1] = Sh->IA[ci] + count;
1161
1162 } // end for i
1163
1164 /*************************/
1165 /* step 2: Find JA of Sh */
1166 /*************************/
1167
1168 fasp_iarray_set(num_c, visited, -1); // reset visited
1169
1170 Sh->nnz = Sh->IA[Sh->row];
1171 Sh->JA = (INT*)fasp_mem_calloc(Sh->nnz, sizeof(INT));
1172
1173 for (ci = 0; ci < Sh->row; ci++) {
1174
1175 i = cp_index[ci]; // find the index of the i-th coarse grid point
1176 count = Sh->IA[ci]; // count for coarse points
1177
1178 // visit all the fine neighbors that ci is strongly connected to
1179 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
1180
1181 fj = S->JA[j];
1182
1183 if (vec[fj] == CGPT && fj != i) {
1184 cj = cp_rindex[fj];
1185 if (visited[cj] != ci) { // not visited yet
1186 visited[cj] = ci;
1187 Sh->JA[count] = cj;
1188 count++;
1189 }
1190 } else if (vec[fj] == FGPT) { // fine grid point,
1191 // find all the coarse neighbors that fj is strongly connected to
1192 for (k = S->IA[fj]; k < S->IA[fj + 1]; k++) {
1193 ck = S->JA[k];
1194 if (vec[ck] == CGPT && ck != i) { // coarse grid point
1195 cck = cp_rindex[ck];
1196 if (visited[cck] != ci) { // not visited yet
1197 visited[cck] = ci;
1198 Sh->JA[count] = cck;
1199 count++;
1200 }
1201 } // end if
1202 } // end for k
1203 } // end if
1204
1205 } // end for j
1206
1207 if (count != Sh->IA[ci + 1]) {
1208 printf("### WARNING: Inconsistent numbers of nonzeros!\n ");
1209 }
1210
1211 } // end for ci
1212
1213 fasp_mem_free(visited);
1214 visited = NULL;
1215}
1216
1243static void strong_couplings_agg2(dCSRmat* A,
1244 iCSRmat* S,
1245 iCSRmat* Sh,
1246 ivector* vertices,
1247 ivector* CGPT_index,
1248 ivector* CGPT_rindex)
1249{
1250 const INT row = A->row;
1251
1252 // local variables
1253 INT i, j, k;
1254 INT num_c, count, ci, cj, ck, fj, cck;
1255 INT *cp_index, *cp_rindex, *visited;
1256 INT* vec = vertices->val;
1257
1258 // count the number of coarse grid points
1259 for (num_c = i = 0; i < row; i++) {
1260 if (vec[i] == CGPT) num_c++;
1261 }
1262
1263 // for the reverse indexing of coarse grid points
1264 fasp_ivec_alloc(row, CGPT_rindex);
1265 cp_rindex = CGPT_rindex->val;
1266
1267 // generate coarse grid point index
1268 fasp_ivec_alloc(num_c, CGPT_index);
1269 cp_index = CGPT_index->val;
1270 for (j = i = 0; i < row; i++) {
1271 if (vec[i] == CGPT) {
1272 cp_index[j] = i;
1273 cp_rindex[i] = j;
1274 j++;
1275 }
1276 }
1277
1278 // allocate space for Sh
1279 Sh->row = Sh->col = num_c;
1280 Sh->val = Sh->JA = NULL;
1281 Sh->IA = (INT*)fasp_mem_calloc(Sh->row + 1, sizeof(INT));
1282
1283 // record the number of times some coarse point is visited
1284 visited = (INT*)fasp_mem_calloc(num_c, sizeof(INT));
1285 memset(visited, 0, sizeof(INT) * num_c);
1286
1287 /**********************************************/
1288 /* step 1: Find first the structure IA of Sh */
1289 /**********************************************/
1290
1291 Sh->IA[0] = 0;
1292
1293 for (ci = 0; ci < Sh->row; ci++) {
1294
1295 i = cp_index[ci]; // find the index of the ci-th coarse grid point
1296
1297 // number of coarse point that i is strongly connected to w.r.t. S(p,2)
1298 count = 0;
1299
1300 // visit all the fine neighbors that ci is strongly connected to
1301 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
1302
1303 fj = S->JA[j];
1304
1305 if (vec[fj] == CGPT && fj != i) {
1306 cj = cp_rindex[fj];
1307 if (visited[cj] != ci + 1) { // not visited yet
1308 visited[cj] = ci + 1; // mark as strongly connected from ci
1309 count++;
1310 }
1311 }
1312
1313 else if (vec[fj] == FGPT) { // fine grid point
1314
1315 // find all the coarse neighbors that fj is strongly connected to
1316 for (k = S->IA[fj]; k < S->IA[fj + 1]; k++) {
1317
1318 ck = S->JA[k];
1319
1320 if (vec[ck] == CGPT && ck != i) { // coarse grid point
1321 if (cp_rindex[ck] >= num_c) {
1322 printf("### ERROR: ck=%d, num_c=%d, out of bound!\n", ck,
1323 num_c);
1324 fasp_chkerr(ERROR_AMG_COARSEING, __FUNCTION__);
1325 }
1326 cck = cp_rindex[ck];
1327
1328 if (visited[cck] == ci + 1) {
1329 // visited already!
1330 } else if (visited[cck] == -ci - 1) {
1331 visited[cck] = ci + 1; // mark as strongly connected from ci
1332 count++;
1333 } else {
1334 visited[cck] = -ci - 1; // mark as visited
1335 }
1336
1337 } // end if vec[ck]
1338
1339 } // end for k
1340
1341 } // end if vec[fj]
1342
1343 } // end for j
1344
1345 Sh->IA[ci + 1] = Sh->IA[ci] + count;
1346
1347 } // end for i
1348
1349 /*************************/
1350 /* step 2: Find JA of Sh */
1351 /*************************/
1352
1353 memset(visited, 0, sizeof(INT) * num_c); // reset visited
1354
1355 Sh->nnz = Sh->IA[Sh->row];
1356 Sh->JA = (INT*)fasp_mem_calloc(Sh->nnz, sizeof(INT));
1357
1358 for (ci = 0; ci < Sh->row; ci++) {
1359
1360 i = cp_index[ci]; // find the index of the i-th coarse grid point
1361 count = Sh->IA[ci]; // count for coarse points
1362
1363 // visit all the fine neighbors that ci is strongly connected to
1364 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
1365
1366 fj = S->JA[j];
1367
1368 if (vec[fj] == CGPT && fj != i) {
1369 cj = cp_rindex[fj];
1370 if (visited[cj] != ci + 1) { // not visited yet
1371 visited[cj] = ci + 1;
1372 Sh->JA[count] = cj;
1373 count++;
1374 }
1375 }
1376
1377 else if (vec[fj] == FGPT) { // fine grid point
1378
1379 // find all the coarse neighbors that fj is strongly connected to
1380 for (k = S->IA[fj]; k < S->IA[fj + 1]; k++) {
1381
1382 ck = S->JA[k];
1383
1384 if (vec[ck] == CGPT && ck != i) { // coarse grid point
1385 cck = cp_rindex[ck];
1386 if (visited[cck] == ci + 1) {
1387 // visited before
1388 } else if (visited[cck] == -ci - 1) {
1389 visited[cck] = ci + 1;
1390 Sh->JA[count] = cck;
1391 count++;
1392 } else {
1393 visited[cck] = -ci - 1;
1394 }
1395 } // end if vec[ck]
1396
1397 } // end for k
1398
1399 } // end if vec[fj]
1400
1401 } // end for j
1402
1403 if (count != Sh->IA[ci + 1]) {
1404 printf("### WARNING: Inconsistent numbers of nonzeros!\n ");
1405 }
1406
1407 } // end for ci
1408
1409 fasp_mem_free(visited);
1410 visited = NULL;
1411}
1412
1434static INT
1435cfsplitting_agg(dCSRmat* A, iCSRmat* S, ivector* vertices, INT aggressive_path)
1436{
1437 const INT row = A->row;
1438 INT col = 0; // initialize col(P): returning output
1439
1440 // local variables
1441 INT * vec = vertices->val, *cp_index;
1442 INT maxmeas, maxnode, num_left = 0;
1443 INT measure, newmeas;
1444 INT i, j, k, l, m, ci, cj, ck, cl, num_c;
1445 SHORT IS_CNEIGH;
1446
1447 INT* work = (INT*)fasp_mem_calloc(3 * row, sizeof(INT));
1448 INT *lists = work, *where = lists + row, *lambda = where + row;
1449
1450 ivector CGPT_index, CGPT_rindex;
1451 LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
1452
1453 // Sh is for the strong coupling matrix between temporary CGPTs
1454 // ShT is the transpose of Sh
1455 // Snew is for combining the information from S and Sh
1456 iCSRmat ST, Sh, ShT;
1457
1458#if DEBUG_MODE > 0
1459 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1460#endif
1461
1462 /************************************************************/
1463 /* Coarsening Phase ONE: find temporary coarse level points */
1464 /************************************************************/
1465
1466 num_c = cfsplitting_cls(A, S, vertices);
1467 fasp_icsr_trans(S, &ST);
1468
1469 /************************************************************/
1470 /* Coarsening Phase TWO: find real coarse level points */
1471 /************************************************************/
1472
1473 // find Sh, the strong coupling between coarse grid points S(path,2)
1474 if (aggressive_path < 2)
1475 strong_couplings_agg1(A, S, &Sh, vertices, &CGPT_index, &CGPT_rindex);
1476 else
1477 strong_couplings_agg2(A, S, &Sh, vertices, &CGPT_index, &CGPT_rindex);
1478
1479 fasp_icsr_trans(&Sh, &ShT);
1480
1481 CGPT_index.row = num_c;
1482 CGPT_rindex.row = row;
1483 cp_index = CGPT_index.val;
1484
1485 // 1. Initialize lambda
1486#ifdef _OPENMP
1487#pragma omp parallel for if (num_c > OPENMP_HOLDS)
1488#endif
1489 for (ci = 0; ci < num_c; ++ci) lambda[ci] = ShT.IA[ci + 1] - ShT.IA[ci];
1490
1491 // 2. Form linked list for lambda (max to min)
1492 for (ci = 0; ci < num_c; ++ci) {
1493
1494 i = cp_index[ci];
1495 measure = lambda[ci];
1496
1497 if (vec[i] == ISPT) continue; // skip isolated points
1498
1499 if (measure > 0) {
1500 enter_list(&LoL_head, &LoL_tail, lambda[ci], ci, lists, where);
1501 num_left++;
1502 } else {
1503 if (measure < 0) printf("### WARNING: Negative lambda[%d]!\n", i);
1504
1505 vec[i] = FGPT; // set i as fine node
1506
1507 // update the lambda value in the CGPT neighbor of i
1508 for (ck = Sh.IA[ci]; ck < Sh.IA[ci + 1]; ++ck) {
1509
1510 cj = Sh.JA[ck];
1511 j = cp_index[cj];
1512
1513 if (vec[j] == ISPT) continue;
1514
1515 if (cj < ci) {
1516 newmeas = lambda[cj];
1517 if (newmeas > 0) {
1518 remove_node(&LoL_head, &LoL_tail, newmeas, cj, lists, where);
1519 num_left--;
1520 }
1521 newmeas = ++(lambda[cj]);
1522 enter_list(&LoL_head, &LoL_tail, newmeas, cj, lists, where);
1523 num_left++;
1524 } else {
1525 newmeas = ++(lambda[cj]);
1526 } // end if cj<ci
1527
1528 } // end for ck
1529
1530 } // end if
1531
1532 } // end for ci
1533
1534 // 3. Main loop
1535 while (num_left > 0) {
1536
1537 // pick $i\in U$ with $\max\lambda_i: C:=C\cup\{i\}, U:=U\\{i\}$
1538 maxnode = LoL_head->head;
1539 maxmeas = lambda[maxnode];
1540 if (maxmeas == 0) printf("### WARNING: Head of the list has measure 0!\n");
1541
1542 // mark maxnode as real coarse node, labelled as number 3
1543 vec[cp_index[maxnode]] = 3;
1544 --num_left;
1545 remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
1546 lambda[maxnode] = 0;
1547 col++; // count for the real coarse node after aggressive coarsening
1548
1549 // for all $j\in S_i^T\cap U: F:=F\cup\{j\}, U:=U\backslash\{j\}$
1550 for (ci = ShT.IA[maxnode]; ci < ShT.IA[maxnode + 1]; ++ci) {
1551
1552 cj = ShT.JA[ci];
1553 j = cp_index[cj];
1554
1555 if (vec[j] != CGPT) continue; // skip if j is not C-point
1556
1557 vec[j] = 4; // set j as 4--fake CGPT
1558 remove_node(&LoL_head, &LoL_tail, lambda[cj], cj, lists, where);
1559 --num_left;
1560
1561 // update the measure for neighboring points
1562 for (cl = Sh.IA[cj]; cl < Sh.IA[cj + 1]; cl++) {
1563 ck = Sh.JA[cl];
1564 k = cp_index[ck];
1565 if (vec[k] == CGPT) { // k is temporary CGPT
1566 remove_node(&LoL_head, &LoL_tail, lambda[ck], ck, lists, where);
1567 newmeas = ++(lambda[ck]);
1568 enter_list(&LoL_head, &LoL_tail, newmeas, ck, lists, where);
1569 }
1570 }
1571
1572 } // end for ci
1573
1574 // Update lambda and linked list after maxnode->C
1575 for (ci = Sh.IA[maxnode]; ci < Sh.IA[maxnode + 1]; ++ci) {
1576
1577 cj = Sh.JA[ci];
1578 j = cp_index[cj];
1579
1580 if (vec[j] != CGPT) continue; // skip if j is not C-point
1581
1582 measure = lambda[cj];
1583 remove_node(&LoL_head, &LoL_tail, measure, cj, lists, where);
1584 lambda[cj] = --measure;
1585
1586 if (measure > 0) {
1587 enter_list(&LoL_head, &LoL_tail, measure, cj, lists, where);
1588 } else {
1589 vec[j] = 4; // set j as fake CGPT variable
1590 --num_left;
1591 for (cl = Sh.IA[cj]; cl < Sh.IA[cj + 1]; cl++) {
1592 ck = Sh.JA[cl];
1593 k = cp_index[ck];
1594 if (vec[k] == CGPT) { // k is temporary CGPT
1595 remove_node(&LoL_head, &LoL_tail, lambda[ck], ck, lists, where);
1596 newmeas = ++(lambda[ck]);
1597 enter_list(&LoL_head, &LoL_tail, newmeas, ck, lists, where);
1598 }
1599 } // end for l
1600 } // end if
1601
1602 } // end for
1603
1604 } // while
1605
1606 // 4. reorganize the variable type: mark temporary CGPT--1 and fake CGPT--4 as
1607 // FGPT; mark real CGPT--3 to be CGPT
1608#ifdef _OPENMP
1609#pragma omp parallel for if (row > OPENMP_HOLDS)
1610#endif
1611 for (i = 0; i < row; i++) {
1612 if (vec[i] == CGPT || vec[i] == 4) vec[i] = FGPT;
1613 }
1614
1615#ifdef _OPENMP
1616#pragma omp parallel for if (row > OPENMP_HOLDS)
1617#endif
1618 for (i = 0; i < row; i++) {
1619 if (vec[i] == 3) vec[i] = CGPT;
1620 }
1621
1622 /************************************************************/
1623 /* Coarsening Phase THREE: all the FGPTs which have no CGPT */
1624 /* neighbors within distance 2. Change them into CGPT such */
1625 /* that the standard interpolation works! */
1626 /************************************************************/
1627
1628 for (i = 0; i < row; i++) {
1629
1630 if (vec[i] != FGPT) continue;
1631
1632 IS_CNEIGH = FALSE; // whether there exist CGPT neighbors within distance of 2
1633
1634 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
1635
1636 if (IS_CNEIGH) break;
1637
1638 k = S->JA[j];
1639
1640 if (vec[k] == CGPT) {
1641 IS_CNEIGH = TRUE;
1642 } else if (vec[k] == FGPT) {
1643 for (l = S->IA[k]; l < S->IA[k + 1]; l++) {
1644 m = S->JA[l];
1645 if (vec[m] == CGPT) {
1646 IS_CNEIGH = TRUE;
1647 break;
1648 }
1649 } // end for l
1650 }
1651
1652 } // end for j
1653
1654 // no CGPT neighbors in distance <= 2, mark i as CGPT
1655 if (!IS_CNEIGH) {
1656 vec[i] = CGPT;
1657 col++;
1658 }
1659
1660 } // end for i
1661
1662 if (LoL_head) {
1663 list_ptr = LoL_head;
1664 LoL_head->prev_node = NULL;
1665 LoL_head->next_node = NULL;
1666 LoL_head = list_ptr->next_node;
1667 fasp_mem_free(list_ptr);
1668 list_ptr = NULL;
1669 }
1670
1671 fasp_ivec_free(&CGPT_index);
1672 fasp_ivec_free(&CGPT_rindex);
1673 fasp_icsr_free(&Sh);
1674 fasp_icsr_free(&ST);
1675 fasp_icsr_free(&ShT);
1676 fasp_mem_free(work);
1677 work = NULL;
1678
1679#if DEBUG_MODE > 0
1680 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1681#endif
1682
1683 return col;
1684}
1685
1709static INT clean_ff_couplings(iCSRmat* S, ivector* vertices, INT row, INT col)
1710{
1711 // local variables
1712 INT* vec = vertices->val;
1713 INT* cindex = (INT*)fasp_mem_calloc(row, sizeof(INT));
1714 INT set_empty = TRUE, C_i_nonempty = FALSE;
1715 INT ci_tilde = -1, ci_tilde_mark = -1;
1716 INT ji, jj, i, j, index;
1717
1718 fasp_iarray_set(row, cindex, -1);
1719
1720 for (i = 0; i < row; ++i) {
1721
1722 if (vec[i] != FGPT) continue; // skip non F-variables
1723
1724 for (ji = S->IA[i]; ji < S->IA[i + 1]; ++ji) {
1725 j = S->JA[ji];
1726 if (vec[j] == CGPT)
1727 cindex[j] = i; // mark C-neighbors
1728 else
1729 cindex[j] = -1; // reset cindex --Chensong 06/02/2013
1730 }
1731
1732 if (ci_tilde_mark != i) ci_tilde = -1; //???
1733
1734 for (ji = S->IA[i]; ji < S->IA[i + 1]; ++ji) {
1735
1736 j = S->JA[ji];
1737
1738 if (vec[j] != FGPT) continue; // skip non F-variables
1739
1740 // check whether there is a C-connection
1741 set_empty = TRUE;
1742 for (jj = S->IA[j]; jj < S->IA[j + 1]; ++jj) {
1743 index = S->JA[jj];
1744 if (cindex[index] == i) {
1745 set_empty = FALSE;
1746 break;
1747 }
1748 } // end for jj
1749
1750 // change the point i (if only F-F exists) to C
1751 if (set_empty) {
1752 if (C_i_nonempty) {
1753 vec[i] = CGPT;
1754 col++;
1755 if (ci_tilde > -1) {
1756 vec[ci_tilde] = FGPT;
1757 col--;
1758 ci_tilde = -1;
1759 }
1760 C_i_nonempty = FALSE;
1761 break;
1762 } else { // temporary set j->C and roll back
1763 vec[j] = CGPT;
1764 col++;
1765 ci_tilde = j;
1766 ci_tilde_mark = i;
1767 C_i_nonempty = TRUE;
1768 i--; // roll back to check i-point again
1769 break;
1770 } // end if C_i_nonempty
1771 } // end if set_empty
1772
1773 } // end for ji
1774
1775 } // end for i
1776
1777 fasp_mem_free(cindex);
1778 cindex = NULL;
1779
1780 return col;
1781}
1782
1783REAL rabs(REAL x) { return (x > 0) ? x : -x; }
1796static void form_P_pattern_rdc(
1797 dCSRmat* P, dCSRmat* A, double* theta, ivector* vertices, INT row, INT col)
1798{
1799 // local variables
1800 INT* vec = vertices->val;
1801 INT i, j, k, index = 0;
1802
1803 // Initialize P matrix
1804 P->row = row;
1805 P->col = col;
1806 P->IA = (INT*)fasp_mem_calloc(row + 1, sizeof(INT));
1807 P->IA[0] = 0;
1808
1809 /* Generate sparsity pattern of P & calculate theta */
1810 // firt pass: P->IA & theta
1811 for (i = 0; i < row; ++i) {
1812 if (vec[i] == CGPT) { // identity interpolation for C-points
1813 P->IA[i + 1] = P->IA[i] + 1;
1814 theta[i] = 1000000.0;
1815 } else { // D_FF^-1 * A_FC for F-points
1816 P->IA[i + 1] = P->IA[i];
1817 double sum = 0.0;
1818 INT diagptr = -1;
1819 for (k = A->IA[i]; k < A->IA[i + 1]; ++k) {
1820 j = A->JA[k];
1821 if (vec[j] == CGPT)
1822 P->IA[i + 1]++;
1823 else {
1824 sum += rabs(A->val[k]);
1825 }
1826 if (j == i) diagptr = k;
1827 }
1828 if (diagptr > -1)
1829 theta[i] = rabs(A->val[diagptr]) / sum;
1830 else {
1831 printf("### ERROR: Diagonal element is zero! [%s:%d]\n", __FUNCTION__,
1832 __LINE__);
1833 theta[i] = 1000000.0;
1834 }
1835 if (theta[i] < 0.5) {
1836 printf("WARNING: theta[%d] = %f < 0.5! not diagonal dominant!\n", i,
1837 theta[i]);
1838 // view the matrix row entries
1839 // printf("### DEBUG: row %d: ", i);
1840 int ii;
1841 int jj;
1842 double sum_row = 0.0;
1843 for (ii = A->IA[i]; ii < A->IA[i + 1]; ++ii) {
1844 jj = A->JA[ii];
1845 sum_row += A->val[ii];
1846 // printf("A[%d,%d] = %f, ", i, jj, A->val[ii]);
1847 }
1848 printf("(no abs op)sum_row = %f\n", sum_row);
1849 }
1850 }
1851 }
1852
1853 // second pass: P->JA
1854 // TODO: given nnz, we can combine the two passes
1855 P->nnz = P->IA[row];
1856 P->JA = (INT*)fasp_mem_calloc(P->nnz, sizeof(INT));
1857 for (i = 0; i < row; ++i) {
1858 if (vec[i] == CGPT) { // identity interpolation for C-points
1859 P->JA[index++] = i; // use index on fine grid (need to be replaced with
1860 // coarse grid index later)
1861 } else { // D_FF^-1 * A_FC for F-points
1862 for (k = A->IA[i]; k < A->IA[i + 1]; ++k) {
1863 j = A->JA[k];
1864 if (vec[j] == CGPT) P->JA[index++] = j;
1865 }
1866 }
1867 }
1868
1869 // val allocated here
1870 P->val = (REAL*)fasp_mem_calloc(P->nnz, sizeof(REAL));
1871}
1872
1891static void
1892form_P_pattern_dir(dCSRmat* P, iCSRmat* S, ivector* vertices, INT row, INT col)
1893{
1894 // local variables
1895 INT i, j, k, index;
1896 INT* vec = vertices->val;
1897
1898 SHORT nthreads = 1, use_openmp = FALSE;
1899
1900#ifdef _OPENMP
1901 if (row > OPENMP_HOLDS) {
1902 use_openmp = TRUE;
1903 nthreads = fasp_get_num_threads();
1904 }
1905#endif
1906
1907 // Initialize P matrix
1908 P->row = row;
1909 P->col = col;
1910 P->IA = (INT*)fasp_mem_calloc(row + 1, sizeof(INT));
1911
1912 // step 1: Find the structure IA of P first: using P as a counter
1913 if (use_openmp) {
1914
1915 INT mybegin, myend, myid;
1916#ifdef _OPENMP
1917#pragma omp parallel for private(myid, mybegin, myend, i, j, k)
1918#endif
1919 for (myid = 0; myid < nthreads; myid++) {
1920 fasp_get_start_end(myid, nthreads, row, &mybegin, &myend);
1921 for (i = mybegin; i < myend; ++i) {
1922 switch (vec[i]) {
1923 case FGPT: // fine grid points
1924 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
1925 k = S->JA[j];
1926 if (vec[k] == CGPT) P->IA[i + 1]++;
1927 }
1928 break;
1929
1930 case CGPT: // coarse grid points
1931 P->IA[i + 1] = 1;
1932 break;
1933
1934 default: // treat everything else as isolated
1935 P->IA[i + 1] = 0;
1936 break;
1937 }
1938 }
1939 }
1940
1941 }
1942
1943 else {
1944
1945 for (i = 0; i < row; ++i) {
1946 switch (vec[i]) {
1947 case FGPT: // fine grid points
1948 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
1949 k = S->JA[j];
1950 if (vec[k] == CGPT) P->IA[i + 1]++;
1951 }
1952 break;
1953
1954 case CGPT: // coarse grid points
1955 P->IA[i + 1] = 1;
1956 break;
1957
1958 default: // treat everything else as isolated
1959 P->IA[i + 1] = 0;
1960 break;
1961 }
1962 } // end for i
1963
1964 } // end if
1965
1966 // Form P->IA from the counter P
1967 for (i = 0; i < P->row; ++i) P->IA[i + 1] += P->IA[i];
1968 P->nnz = P->IA[P->row] - P->IA[0];
1969
1970 // step 2: Find the structure JA of P
1971 P->JA = (INT*)fasp_mem_calloc(P->nnz, sizeof(INT));
1972 P->val = (REAL*)fasp_mem_calloc(P->nnz, sizeof(REAL));
1973
1974 for (index = i = 0; i < row; ++i) {
1975 if (vec[i] == FGPT) { // fine grid point
1976 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
1977 k = S->JA[j];
1978 if (vec[k] == CGPT) P->JA[index++] = k;
1979 } // end for j
1980 } // end if
1981 else if (vec[i] == CGPT) { // coarse grid point -- one entry only
1982 P->JA[index++] = i;
1983 }
1984 }
1985}
1986
2005static void
2006form_P_pattern_std(dCSRmat* P, iCSRmat* S, ivector* vertices, INT row, INT col)
2007{
2008 // local variables
2009 INT i, j, k, l, h, index;
2010 INT* vec = vertices->val;
2011
2012 // number of times a C-point is visited
2013 INT* visited = (INT*)fasp_mem_calloc(row, sizeof(INT));
2014
2015 P->row = row;
2016 P->col = col;
2017 P->IA = (INT*)fasp_mem_calloc(row + 1, sizeof(INT));
2018
2019 fasp_iarray_set(row, visited, -1);
2020
2021 // Step 1: Find the structure IA of P first: use P as a counter
2022 for (i = 0; i < row; ++i) {
2023
2024 if (vec[i] == FGPT) { // if node i is a F point
2025 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
2026
2027 k = S->JA[j];
2028
2029 // if neighbor of i is a C point, good
2030 if ((vec[k] == CGPT) && (visited[k] != i)) {
2031 visited[k] = i;
2032 P->IA[i + 1]++;
2033 }
2034
2035 // if k is a F point and k is not i, look for indirect C neighbors
2036 else if ((vec[k] == FGPT) && (k != i)) {
2037 for (l = S->IA[k]; l < S->IA[k + 1]; l++) { // neighbors of k
2038 h = S->JA[l];
2039 if ((vec[h] == CGPT) && (visited[h] != i)) {
2040 visited[h] = i;
2041 P->IA[i + 1]++;
2042 }
2043 } // end for(l=S->IA[k];l<S->IA[k+1];l++)
2044 } // end if (vec[k]==CGPT)
2045
2046 } // end for (j=S->IA[i];j<S->IA[i+1];j++)
2047 }
2048
2049 else if (vec[i] == CGPT) { // if node i is a C point
2050 P->IA[i + 1] = 1;
2051 }
2052
2053 else { // treat everything else as isolated points
2054 P->IA[i + 1] = 0;
2055 } // end if (vec[i]==FGPT)
2056
2057 } // end for (i=0;i<row;++i)
2058
2059 // Form P->IA from the counter P
2060 for (i = 0; i < P->row; ++i) P->IA[i + 1] += P->IA[i];
2061 P->nnz = P->IA[P->row] - P->IA[0];
2062
2063 // Step 2: Find the structure JA of P
2064 P->JA = (INT*)fasp_mem_calloc(P->nnz, sizeof(INT));
2065 P->val = (REAL*)fasp_mem_calloc(P->nnz, sizeof(REAL));
2066
2067 fasp_iarray_set(row, visited, -1); // re-init visited array
2068
2069 for (i = 0; i < row; ++i) {
2070
2071 if (vec[i] == FGPT) { // if node i is a F point
2072
2073 index = 0;
2074
2075 for (j = S->IA[i]; j < S->IA[i + 1]; j++) {
2076
2077 k = S->JA[j];
2078
2079 // if neighbor k of i is a C point
2080 if ((vec[k] == CGPT) && (visited[k] != i)) {
2081 visited[k] = i;
2082 P->JA[P->IA[i] + index] = k;
2083 index++;
2084 }
2085
2086 // if neighbor k of i is a F point and k is not i
2087 else if ((vec[k] == FGPT) && (k != i)) {
2088 for (l = S->IA[k]; l < S->IA[k + 1]; l++) { // neighbors of k
2089 h = S->JA[l];
2090 if ((vec[h] == CGPT) && (visited[h] != i)) {
2091 visited[h] = i;
2092 P->JA[P->IA[i] + index] = h;
2093 index++;
2094 }
2095
2096 } // end for (l=S->IA[k];l<S->IA[k+1];l++)
2097
2098 } // end if (vec[k]==CGPT)
2099
2100 } // end for (j=S->IA[i];j<S->IA[i+1];j++)
2101 }
2102
2103 else if (vec[i] == CGPT) {
2104 P->JA[P->IA[i]] = i;
2105 }
2106 }
2107
2108 // clean up
2109 fasp_mem_free(visited);
2110 visited = NULL;
2111}
2112
2127static INT cfsplitting_mis(iCSRmat* S, ivector* vertices, ivector* order)
2128{
2129 const INT n = S->row;
2130
2131 INT col = 0;
2132 INT* ord = order->val;
2133 INT* vec = vertices->val;
2134 INT* IS = S->IA;
2135 INT* JS = S->JA;
2136
2137 INT i, j, ind;
2138 INT row_begin, row_end;
2139
2140 fasp_ivec_set(n, vertices, UNPT);
2141
2142 for (i = 0; i < n; i++) {
2143 ind = ord[i];
2144 if (vec[ind] == UNPT) {
2145 vec[ind] = CGPT;
2146 row_begin = IS[ind];
2147 row_end = IS[ind + 1];
2148 for (j = row_begin; j < row_end; j++) {
2149 if (vec[JS[j]] == CGPT) {
2150 vec[ind] = FGPT;
2151 break;
2152 }
2153 }
2154 if (vec[ind] == CGPT) {
2155 col++;
2156 for (j = row_begin; j < row_end; j++) {
2157 vec[JS[j]] = FGPT;
2158 }
2159 }
2160 }
2161 }
2162 return col;
2163}
2164
2179static void ordering1(iCSRmat* S, ivector* order)
2180{
2181 const INT n = order->row;
2182 INT* IS = S->IA;
2183 INT* ord = order->val;
2184 INT maxind, maxdeg, degree;
2185 INT i;
2186
2187 for (i = 0; i < n; i++) ord[i] = i;
2188
2189 for (maxind = maxdeg = i = 0; i < n; i++) {
2190 degree = IS[i + 1] - IS[i];
2191 if (degree > maxdeg) {
2192 maxind = i;
2193 maxdeg = degree;
2194 }
2195 }
2196
2197 ord[0] = maxind;
2198 ord[maxind] = 0;
2199
2200 return;
2201}
2202
2203/*---------------------------------*/
2204/*-- End of File --*/
2205/*---------------------------------*/
fasp_iarray_set
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:98
fasp_iarray_cp
void fasp_iarray_cp(const INT n, const INT *x, INT *y)
Copy an array to the other y=x.
Definition: AuxArray.c:227
fasp_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_chkerr
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
fasp_get_start_end
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:93
fasp_dvec_free
void fasp_dvec_free(dvector *u)
Free vector data space of REAL type.
Definition: AuxVector.c:145
fasp_ivec_free
void fasp_ivec_free(ivector *u)
Free vector data space of INT type.
Definition: AuxVector.c:164
fasp_ivec_alloc
void fasp_ivec_alloc(const INT m, ivector *u)
Create vector data space of INT type.
Definition: AuxVector.c:125
fasp_ivec_create
ivector fasp_ivec_create(const INT m)
Create vector data space of INT type.
Definition: AuxVector.c:84
fasp_ivec_set
void fasp_ivec_set(INT n, ivector *u, const INT m)
Set ivector value to be m.
Definition: AuxVector.c:291
fasp_icsr_trans
void fasp_icsr_trans(const iCSRmat *A, iCSRmat *AT)
Find transpose of iCSRmat matrix A.
Definition: BlaSparseCSR.c:875
fasp_icsr_free
void fasp_icsr_free(iCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:219
fasp_dcsr_getdiag
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
fasp_amg_coarsening_cr
INT fasp_amg_coarsening_cr(const INT i_0, const INT i_n, dCSRmat *A, ivector *vertices, AMG_param *param)
CR coarsening.
Definition: PreAMGCoarsenCR.c:62
fasp_amg_coarsening_rs
SHORT fasp_amg_coarsening_rs(dCSRmat *A, ivector *vertices, dCSRmat *P, iCSRmat *S, AMG_param *param)
Standard and aggressive coarsening schemes.
Definition: PreAMGCoarsenRS.c:76
fasp.h
Main header file for the FASP project.
MIN
#define MIN(a, b)
Definition: fasp.h:83
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
ABS
#define ABS(a)
Definition: fasp.h:84
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
ISPT
#define ISPT
Definition: fasp_const.h:235
COARSE_MIS
#define COARSE_MIS
Definition: fasp_const.h:217
ERROR_AMG_INTERP_TYPE
#define ERROR_AMG_INTERP_TYPE
Definition: fasp_const.h:35
FASP_SUCCESS
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
FGPT
#define FGPT
Definition: fasp_const.h:233
COARSE_RSP
#define COARSE_RSP
Definition: fasp_const.h:214
INTERP_EXT
#define INTERP_EXT
Definition: fasp_const.h:225
CGPT
#define CGPT
Definition: fasp_const.h:234
INTERP_RDC
#define INTERP_RDC
Definition: fasp_const.h:226
COARSE_AC
#define COARSE_AC
Definition: fasp_const.h:216
UNPT
#define UNPT
Definition: fasp_const.h:232
INTERP_STD
#define INTERP_STD
Definition: fasp_const.h:223
INTERP_DIR
#define INTERP_DIR
Definition of interpolation types.
Definition: fasp_const.h:222
OPENMP_HOLDS
#define OPENMP_HOLDS
Definition: fasp_const.h:269
TRUE
#define TRUE
Definition of logic type.
Definition: fasp_const.h:61
FALSE
#define FALSE
Definition: fasp_const.h:62
ERROR_AMG_COARSEING
#define ERROR_AMG_COARSEING
Definition: fasp_const.h:38
COARSE_CR
#define COARSE_CR
Definition: fasp_const.h:215
SMOOTHER_JACOBIF
#define SMOOTHER_JACOBIF
Definition of standard smoother types.
Definition: fasp_const.h:189
ERROR_UNKNOWN
#define ERROR_UNKNOWN
Definition: fasp_const.h:56
INTERP_ENG
#define INTERP_ENG
Definition: fasp_const.h:224
AMG_param
Parameters for AMG methods.
Definition: fasp.h:455
AMG_param::strong_threshold
REAL strong_threshold
strong connection threshold for coarsening
Definition: fasp.h:527
AMG_param::aggressive_path
INT aggressive_path
number of paths use to determine strongly coupled C points
Definition: fasp.h:539
AMG_param::interpolation_type
SHORT interpolation_type
interpolation type
Definition: fasp.h:524
AMG_param::coarsening_type
SHORT coarsening_type
coarsening type
Definition: fasp.h:515
AMG_param::theta
REAL theta
theta for reduction-based amg
Definition: fasp.h:593
AMG_param::relaxation
REAL relaxation
relaxation parameter for Jacobi and SOR smoother
Definition: fasp.h:494
AMG_param::smoother
SHORT smoother
smoother type
Definition: fasp.h:482
AMG_param::postsmooth_iter
SHORT postsmooth_iter
number of postsmoothers
Definition: fasp.h:491
AMG_param::presmooth_iter
SHORT presmooth_iter
number of presmoothers
Definition: fasp.h:488
AMG_param::max_row_sum
REAL max_row_sum
maximal row sum parameter
Definition: fasp.h:530
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:157
dCSRmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:163
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:166
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360
iCSRmat
Sparse matrix of INT type in CSR format.
Definition: fasp.h:190
iCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:196
iCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:193
iCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:202
iCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:199
iCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:205
iCSRmat::val
INT * val
nonzero entries of A
Definition: fasp.h:208
ivector
Vector with n entries of INT type.
Definition: fasp.h:368
ivector::row
INT row
number of rows
Definition: fasp.h:371
ivector::val
INT * val
actual vector entries
Definition: fasp.h:374