25#include "fasp_functs.h"
31#include "PreAMGUtil.inl"
88 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
92 printf(
"### DEBUG: Step 1. Find strong connections ......\n");
99 strong_couplings(A, S, param);
102 printf(
"### DEBUG: Step 2. C/F splitting ......\n");
105 switch ( coarse_type ) {
108 col = cfsplitting_clsp(A, S, vertices);
break;
111 col = cfsplitting_agg(A, S, vertices, agg_path);
break;
120 ordering1(S, &order);
121 col = cfsplitting_mis(S, vertices, &order);
127 col = cfsplitting_cls(A, S, vertices);
132 printf(
"### DEBUG: col = %d\n", col);
137 printf(
"### DEBUG: Step 3. Find support of C points ......\n");
140 switch ( interp_type ) {
144 col = clean_ff_couplings(S, vertices, row, col);
145 form_P_pattern_dir(P, S, vertices, row, col);
150 form_P_pattern_std(P, S, vertices, row, col);
break;
158 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
186static void strong_couplings (
dCSRmat *A,
193 const INT row = A->
row, col = A->
col, row1 = row+1;
200 INT i, j, begin_row, end_row;
201 REAL row_scl, row_sum;
208 nthreads = fasp_get_num_threads();
216 S->
row = row; S->
col = col; S->
nnz = nnz; S->
val = NULL;
226 INT mybegin, myend, myid;
228#pragma omp parallel for private(myid, mybegin,myend,i,row_scl,row_sum,begin_row,end_row,j)
230 for ( myid = 0; myid < nthreads; myid++ ) {
232 for ( i = mybegin; i < myend; i++) {
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]);
243 for ( j = begin_row; j < end_row; j++ ) {
244 if ( ja[j] == i ) { S->
JA[j] = -1;
break; }
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;
252 for ( j = begin_row; j < end_row; j++) {
256 if ( A->
val[j] >= epsilon_str*row_scl ) S->
JA[j] = -1;
267 for ( i = 0; i < row; ++i ) {
270 row_scl = row_sum = 0.0;
271 begin_row = ia[i]; end_row = ia[i+1];
273 for ( j = begin_row; j < end_row; j++ ) {
278 row_sum +=
ABS(aj[j]);
283 if ( ja[j] != i ) row_scl =
MAX(row_scl,
ABS(aj[j]));
288 row_scl *= epsilon_str;
291 for ( j = begin_row; j < end_row; j++ ) {
292 if ( ja[j] == i ) { S->
JA[j] = -1;
break; }
299 if ( row_sum < (2 - max_row_sum) *
ABS(diag.
val[i]) ) {
301 for ( j = begin_row; j < end_row; j++ ) S->
JA[j] = -1;
306 switch ( coarse_type ) {
309 for ( j = begin_row; j < end_row; j++ ) {
310 if (
ABS(A->
val[j]) <= row_scl ) S->
JA[j] = -1;
315 for ( j = begin_row; j < end_row; j++ ) {
316 if ( -A->
val[j] <= row_scl ) S->
JA[j] = -1;
349 INT index, i, j, begin_row, end_row;
352 for ( index = i = 0; i < row; ++i ) {
354 begin_row = ia[i]; end_row = ia[i+1];
357 for ( j = begin_row; j < end_row; j++ ) {
358 if ( S->
JA[j] > -1 ) S->
JA[index++] = S->
JA[j];
363 S->
nnz = S->
IA[row] = index;
385static void rem_positive_ff (
dCSRmat *A,
390 INT *ia = A->
IA, *vec = vertices->
val;
392 REAL row_scl, max_entry;
393 INT i, j, ji, max_index;
395 for ( i = 0; i < row; ++i ) {
397 if ( vec[i] !=
FGPT )
continue;
400 for ( ji = ia[i]; ji < ia[i+1]; ++ji ) {
402 if ( j == i )
continue;
408 max_index = -1; max_entry = 0.0;
409 for ( ji = ia[i]; ji < ia[i+1]; ++ji ) {
411 if ( j == i )
continue;
412 if ( vec[j] !=
FGPT )
continue;
413 if ( A->
val[ji] > row_scl ) {
415 if ( A->
val[ji] > max_entry ) {
416 max_entry = A->
val[ji];
423 if ( max_index != -1 ) vec[max_index] =
CGPT;
458 INT maxmeas, maxnode, num_left = 0;
459 INT measure, newmeas;
461 INT myid, mybegin, myend;
464 INT *lists = work, *where = lists+row, *lambda = where+row;
468 INT jkeep = 0, cnt, index;
469 INT row_end_S, ji, row_end_S_nabor, jj;
470 INT *graph_array = lambda;
475 LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
480 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
486 nthreads = fasp_get_num_threads();
492 if ( col < 0 )
goto FINISHED;
499#pragma omp parallel for private(myid, mybegin, myend, i)
501 for ( myid = 0; myid < nthreads; myid++ ) {
503 for ( i = mybegin; i < myend; i++ ) lambda[i] = ST.
IA[i+1] - ST.
IA[i];
507 for ( i = 0; i < row; ++i ) lambda[i] = ST.
IA[i+1] - ST.
IA[i];
515#pragma omp parallel for reduction(+:num_left) private(myid, mybegin, myend, i)
517 for ( myid = 0; myid < nthreads; myid++ ) {
519 for ( i = mybegin; i < myend; i++ ) {
521 if ( S->
IA[i+1] == S->
IA[i] )
523 if ( (ia[i+1]-ia[i]) <= 1 )
540 for ( i = 0; i < row; ++i ) {
543 if ( S->
IA[i+1] == S->
IA[i] )
545 if ( (ia[i+1]-ia[i]) <= 1 )
560 for ( i = 0; i < row; ++i ) {
562 if ( vec[i] ==
ISPT )
continue;
567 enter_list(&LoL_head, &LoL_tail, lambda[i], i, lists, where);
571 if ( measure < 0 ) printf(
"### WARNING: Negative lambda[%d]!\n", i);
578 for ( k = S->
IA[i]; k < S->IA[i+1]; ++k ) {
580 if ( vec[j] ==
ISPT )
continue;
584 remove_node(&LoL_head, &LoL_tail, newmeas, j, lists, where);
586 newmeas = ++(lambda[j]);
587 enter_list(&LoL_head, &LoL_tail, newmeas, j, lists, where);
590 newmeas = ++(lambda[j]);
599 while ( num_left > 0 ) {
602 maxnode = LoL_head->head;
603 maxmeas = lambda[maxnode];
605 printf(
"### WARNING: Head of the list has measure 0!\n");
610 remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
614 for ( i = ST.
IA[maxnode]; i < ST.
IA[maxnode+1]; ++i ) {
618 if ( vec[j] !=
UNPT )
continue;
621 remove_node(&LoL_head, &LoL_tail, lambda[j], j, lists, where);
625 for ( l = S->
IA[j]; l < S->IA[j+1]; l++ ) {
627 if ( vec[k] ==
UNPT ) {
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);
637 for ( i = S->
IA[maxnode]; i < S->IA[maxnode+1]; ++i ) {
641 if ( vec[j] !=
UNPT )
continue;
644 remove_node(&LoL_head, &LoL_tail, measure, j, lists, where);
645 lambda[j] = --measure;
648 enter_list(&LoL_head, &LoL_tail, measure, j, lists, where);
655 for ( l = S->
IA[j]; l < S->IA[j+1]; l++ ) {
657 if ( vec[k] ==
UNPT ) {
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);
673 for (i = 0; i < row; i ++)
677 row_end_S = S->
IA[i+1];
678 for (ji = S->
IA[i]; ji < row_end_S; ji ++)
687 for (ji = S->
IA[i]; ji < row_end_S; ji ++)
693 row_end_S_nabor = S->
IA[j+1];
694 for (jj = S->
IA[j]; jj < row_end_S_nabor; jj ++)
697 if (graph_array[index] == i)
731 LoL_head->prev_node = NULL;
732 LoL_head->next_node = NULL;
733 LoL_head = list_ptr->next_node;
741 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
774 INT maxmeas, maxnode, num_left = 0;
775 INT measure, newmeas;
777 INT myid, mybegin, myend;
779 INT *ia = A->
IA, *vec = vertices->
val;
781 INT *lists = work, *where = lists+row, *lambda = where+row;
783 LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
788 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
794 nthreads = fasp_get_num_threads();
806 if ( compress_S(S) < 0 )
goto FINISHED;
813#pragma omp parallel for private(myid, mybegin,myend,i)
815 for ( myid = 0; myid < nthreads; myid++ ) {
817 for ( i = mybegin; i < myend; i++ ) lambda[i] = ST.
IA[i+1] - ST.
IA[i];
821 for ( i = 0; i < row; ++i ) lambda[i] = ST.
IA[i+1] - ST.
IA[i];
829#pragma omp parallel for reduction(+:num_left) private(myid, mybegin, myend, i)
831 for ( myid = 0; myid < nthreads; myid++ ) {
833 for ( i = mybegin; i < myend; i++ ) {
834 if ( (ia[i+1]-ia[i]) <= 1 ) {
848 for ( i = 0; i < row; ++i ) {
849 if ( (ia[i+1]-ia[i]) <= 1 ) {
862 for ( i = 0; i < row; ++i ) {
864 if ( vec[i] ==
ISPT )
continue;
869 enter_list(&LoL_head, &LoL_tail, lambda[i], i, lists, where);
873 if ( measure < 0 ) printf(
"### WARNING: Negative lambda[%d]!\n", i);
880 for ( k = S->
IA[i]; k < S->IA[i+1]; ++k ) {
883 if ( vec[j] ==
ISPT )
continue;
888 remove_node(&LoL_head, &LoL_tail, newmeas, j, lists, where);
890 newmeas = ++(lambda[j]);
891 enter_list(&LoL_head, &LoL_tail, newmeas, j, lists, where);
894 newmeas = ++(lambda[j]);
904 while ( num_left > 0 ) {
907 maxnode = LoL_head->head;
908 maxmeas = lambda[maxnode];
910 printf(
"### WARNING: Head of the list has measure 0!\n");
915 remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
919 for ( i = ST.
IA[maxnode]; i < ST.
IA[maxnode+1]; ++i ) {
923 if ( vec[j] !=
UNPT )
continue;
926 remove_node(&LoL_head, &LoL_tail, lambda[j], j, lists, where);
930 for ( l = S->
IA[j]; l < S->IA[j+1]; l++ ) {
932 if ( vec[k] ==
UNPT ) {
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);
942 for ( i = S->
IA[maxnode]; i < S->IA[maxnode+1]; ++i ) {
946 if ( vec[j] !=
UNPT )
continue;
949 remove_node(&LoL_head, &LoL_tail, measure, j, lists, where);
950 lambda[j] = --measure;
953 enter_list(&LoL_head, &LoL_tail, measure, j, lists, where);
960 for ( l = S->
IA[j]; l < S->IA[j+1]; l++ ) {
962 if ( vec[k] ==
UNPT ) {
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);
978 LoL_head->prev_node = NULL;
979 LoL_head->next_node = NULL;
980 LoL_head = list_ptr->next_node;
988 rem_positive_ff(A, &Stemp, vertices);
990 if ( compress_S(&Stemp) < 0 )
goto FINISHED;
1003 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1029static void strong_couplings_agg1 (
dCSRmat *A,
1040 INT num_c, count, ci, cj, ck, fj, cck;
1041 INT *cp_index, *cp_rindex, *visited;
1042 INT *vec = vertices->
val;
1045 for ( num_c = i = 0; i < row; i++ ) {
1046 if ( vec[i] ==
CGPT ) num_c++;
1051 cp_rindex = CGPT_rindex->
val;
1055 cp_index = CGPT_index->
val;
1056 for ( j = i = 0; i < row; i++ ) {
1057 if ( vec[i] ==
CGPT ) {
1065 Sh->
row = Sh->
col = num_c;
1066 Sh->
val = Sh->
JA = NULL;
1079 for ( ci = 0; ci < Sh->
row; ci++ ) {
1087 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1091 if ( vec[fj] ==
CGPT && fj != i ) {
1093 if ( visited[cj] != ci ) {
1100 else if ( vec[fj] ==
FGPT ) {
1103 for ( k = S->
IA[fj]; k < S->IA[fj+1]; k++ ) {
1105 if ( vec[ck] ==
CGPT && ck != i ) {
1106 if ( cp_rindex[ck] >= num_c ) {
1107 printf(
"### ERROR: ck=%d, num_c=%d, out of bound!\n",
1111 cck = cp_rindex[ck];
1113 if ( visited[cck] != ci ) {
1124 Sh->
IA[ci+1] = Sh->
IA[ci] + count;
1137 for ( ci = 0; ci < Sh->
row; ci++ ) {
1143 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1147 if ( vec[fj] ==
CGPT && fj != i ) {
1149 if ( visited[cj] != ci ) {
1155 else if ( vec[fj] ==
FGPT ) {
1157 for ( k = S->
IA[fj]; k < S->IA[fj+1]; k++ ) {
1159 if ( vec[ck] ==
CGPT && ck != i ) {
1160 cck = cp_rindex[ck];
1161 if ( visited[cck] != ci ) {
1163 Sh->
JA[count] = cck;
1172 if ( count != Sh->
IA[ci+1] ) {
1173 printf(
"### WARNING: Inconsistent numbers of nonzeros!\n ");
1207static void strong_couplings_agg2 (
dCSRmat *A,
1218 INT num_c, count, ci, cj, ck, fj, cck;
1219 INT *cp_index, *cp_rindex, *visited;
1220 INT *vec = vertices->
val;
1223 for ( num_c = i = 0; i < row; i++ ) {
1224 if ( vec[i] ==
CGPT ) num_c++;
1229 cp_rindex = CGPT_rindex->
val;
1233 cp_index = CGPT_index->
val;
1234 for ( j = i = 0; i < row; i++ ) {
1235 if ( vec[i] ==
CGPT ) {
1243 Sh->
row = Sh->
col = num_c;
1244 Sh->
val = Sh->
JA = NULL;
1249 memset(visited, 0,
sizeof(
INT)*num_c);
1257 for ( ci = 0; ci < Sh->
row; ci++ ) {
1265 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1269 if ( vec[fj] ==
CGPT && fj != i ) {
1271 if ( visited[cj] != ci+1 ) {
1277 else if ( vec[fj] ==
FGPT ) {
1280 for ( k = S->
IA[fj]; k < S->IA[fj+1]; k++ ) {
1284 if ( vec[ck] ==
CGPT && ck != i ) {
1285 if ( cp_rindex[ck] >= num_c ) {
1286 printf(
"### ERROR: ck=%d, num_c=%d, out of bound!\n",
1290 cck = cp_rindex[ck];
1292 if ( visited[cck] == ci+1 ) {
1295 else if ( visited[cck] == -ci-1 ) {
1296 visited[cck] = ci+1;
1300 visited[cck] = -ci-1;
1311 Sh->
IA[ci+1] = Sh->
IA[ci] + count;
1319 memset(visited, 0,
sizeof(
INT)*num_c);
1324 for ( ci = 0; ci < Sh->
row; ci++ ) {
1330 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1334 if ( vec[fj] ==
CGPT && fj != i ) {
1336 if ( visited[cj] != ci+1 ) {
1343 else if ( vec[fj] ==
FGPT ) {
1346 for ( k = S->
IA[fj]; k < S->IA[fj+1]; k++ ) {
1350 if ( vec[ck] ==
CGPT && ck != i ) {
1351 cck = cp_rindex[ck];
1352 if ( visited[cck] == ci+1 ) {
1355 else if ( visited[cck] == -ci-1 ) {
1356 visited[cck] = ci+1;
1357 Sh->
JA[count] = cck;
1361 visited[cck] = -ci-1;
1371 if ( count != Sh->
IA[ci+1] ) {
1372 printf(
"### WARNING: Inconsistent numbers of nonzeros!\n ");
1404 INT aggressive_path)
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;
1417 INT *lists = work, *where = lists+row, *lambda = where+row;
1419 ivector CGPT_index, CGPT_rindex;
1420 LinkList LoL_head = NULL, LoL_tail = NULL, list_ptr = NULL;
1428 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1435 num_c = cfsplitting_cls(A, S, vertices);
1443 if ( aggressive_path < 2 )
1444 strong_couplings_agg1(A, S, &Sh, vertices, &CGPT_index, &CGPT_rindex);
1446 strong_couplings_agg2(A, S, &Sh, vertices, &CGPT_index, &CGPT_rindex);
1450 CGPT_index.
row = num_c;
1451 CGPT_rindex.
row = row;
1452 cp_index = CGPT_index.
val;
1456#pragma omp parallel for if(num_c>OPENMP_HOLDS)
1458 for ( ci = 0; ci < num_c; ++ci ) lambda[ci] = ShT.
IA[ci+1]-ShT.
IA[ci];
1461 for ( ci = 0; ci < num_c; ++ci ) {
1464 measure = lambda[ci];
1466 if ( vec[i] ==
ISPT )
continue;
1468 if ( measure > 0 ) {
1469 enter_list(&LoL_head, &LoL_tail, lambda[ci], ci, lists, where);
1473 if ( measure < 0) printf(
"### WARNING: Negative lambda[%d]!\n", i);
1478 for ( ck = Sh.
IA[ci]; ck < Sh.
IA[ci+1]; ++ck ) {
1483 if ( vec[j] ==
ISPT )
continue;
1486 newmeas = lambda[cj];
1487 if ( newmeas > 0 ) {
1488 remove_node(&LoL_head, &LoL_tail, newmeas, cj, lists, where);
1491 newmeas = ++(lambda[cj]);
1492 enter_list(&LoL_head, &LoL_tail, newmeas, cj, lists, where);
1496 newmeas = ++(lambda[cj]);
1506 while ( num_left > 0 ) {
1509 maxnode = LoL_head->head;
1510 maxmeas = lambda[maxnode];
1511 if ( maxmeas == 0 ) printf(
"### WARNING: Head of the list has measure 0!\n");
1514 vec[cp_index[maxnode]] = 3;
1516 remove_node(&LoL_head, &LoL_tail, maxmeas, maxnode, lists, where);
1517 lambda[maxnode] = 0;
1521 for ( ci = ShT.
IA[maxnode]; ci < ShT.
IA[maxnode+1]; ++ci ) {
1526 if ( vec[j] !=
CGPT )
continue;
1529 remove_node(&LoL_head, &LoL_tail, lambda[cj], cj, lists, where);
1533 for ( cl = Sh.
IA[cj]; cl < Sh.
IA[cj+1]; cl++ ) {
1536 if ( vec[k] ==
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);
1546 for ( ci = Sh.
IA[maxnode]; ci < Sh.
IA[maxnode+1]; ++ci ) {
1551 if ( vec[j] !=
CGPT )
continue;
1553 measure = lambda[cj];
1554 remove_node(&LoL_head, &LoL_tail, measure, cj, lists, where);
1555 lambda[cj] = --measure;
1557 if ( measure > 0 ) {
1558 enter_list(&LoL_head, &LoL_tail,measure, cj, lists, where);
1563 for ( cl = Sh.
IA[cj]; cl < Sh.
IA[cj+1]; cl++ ) {
1566 if ( vec[k] ==
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);
1581#pragma omp parallel for if(row>OPENMP_HOLDS)
1583 for ( i = 0; i < row; i++ ) {
1584 if ( vec[i] ==
CGPT || vec[i] == 4 ) vec[i] =
FGPT;
1588#pragma omp parallel for if(row>OPENMP_HOLDS)
1590 for ( i = 0; i < row; i++ ) {
1591 if ( vec[i] == 3 ) vec[i] =
CGPT;
1600 for ( i = 0; i < row; i++ ) {
1602 if ( vec[i] !=
FGPT )
continue;
1606 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1608 if ( IS_CNEIGH )
break;
1612 if ( vec[k] ==
CGPT ) {
1615 else if ( vec[k] ==
FGPT ) {
1616 for ( l = S->
IA[k]; l < S->IA[k+1]; l++ ) {
1618 if ( vec[m] ==
CGPT ) {
1619 IS_CNEIGH =
TRUE;
break;
1628 vec[i] =
CGPT; col++;
1634 list_ptr = LoL_head;
1635 LoL_head->prev_node = NULL;
1636 LoL_head->next_node = NULL;
1637 LoL_head = list_ptr->next_node;
1649 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1684 INT *vec = vertices->
val;
1687 INT ci_tilde = -1, ci_tilde_mark = -1;
1688 INT ji, jj, i, j, index;
1692 for ( i = 0; i < row; ++i ) {
1694 if ( vec[i] !=
FGPT )
continue;
1696 for ( ji = S->
IA[i]; ji < S->IA[i+1]; ++ji ) {
1698 if ( vec[j] ==
CGPT ) cindex[j] = i;
1699 else cindex[j] = -1;
1702 if ( ci_tilde_mark != i ) ci_tilde = -1;
1704 for ( ji = S->
IA[i]; ji < S->IA[i+1]; ++ji ) {
1708 if ( vec[j] !=
FGPT )
continue;
1712 for ( jj = S->
IA[j]; jj < S->IA[j+1]; ++jj ) {
1714 if ( cindex[index] == i ) {
1715 set_empty =
FALSE;
break;
1721 if ( C_i_nonempty ) {
1722 vec[i] =
CGPT; col++;
1723 if ( ci_tilde > -1 ) {
1724 vec[ci_tilde] =
FGPT; col--;
1727 C_i_nonempty =
FALSE;
1731 vec[j] =
CGPT; col++;
1734 C_i_nonempty =
TRUE;
1767static void form_P_pattern_dir (
dCSRmat *P,
1775 INT *vec = vertices->
val;
1782 nthreads = fasp_get_num_threads();
1787 P->
row = row; P->
col = col;
1793 INT mybegin, myend, myid;
1795#pragma omp parallel for private(myid, mybegin,myend,i,j,k)
1797 for ( myid = 0; myid < nthreads; myid++ ) {
1799 for ( i = mybegin; i < myend; ++i ) {
1802 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1804 if ( vec[k] ==
CGPT ) P->
IA[i+1]++;
1809 P->
IA[i+1] = 1;
break;
1812 P->
IA[i+1] = 0;
break;
1821 for ( i = 0; i < row; ++i ) {
1824 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1826 if ( vec[k] ==
CGPT ) P->
IA[i+1]++;
1831 P->
IA[i+1] = 1;
break;
1834 P->
IA[i+1] = 0;
break;
1841 for ( i = 0; i < P->
row; ++i ) P->
IA[i+1] += P->
IA[i];
1848 for ( index = i = 0; i < row; ++i ) {
1849 if ( vec[i] ==
FGPT ) {
1850 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1852 if ( vec[k] ==
CGPT ) P->
JA[index++] = k;
1855 else if ( vec[i] ==
CGPT ) {
1880static void form_P_pattern_std (
dCSRmat *P,
1887 INT i, j, k, l, h, index;
1888 INT *vec = vertices->
val;
1893 P->
row = row; P->
col = col;
1899 for ( i = 0; i < row; ++i ) {
1901 if ( vec[i] ==
FGPT ) {
1902 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1907 if ( (vec[k] ==
CGPT) && (visited[k] != i) ) {
1913 else if ( (vec[k] ==
FGPT) && (k != i) ) {
1914 for ( l = S->
IA[k]; l < S->IA[k+1]; l++ ) {
1916 if ( (vec[h] ==
CGPT) && (visited[h] != i) ) {
1926 else if ( vec[i] ==
CGPT ) {
1937 for ( i = 0; i < P->
row; ++i ) P->
IA[i+1] += P->
IA[i];
1946 for ( i = 0; i < row; ++i ) {
1948 if ( vec[i] ==
FGPT ) {
1952 for ( j = S->
IA[i]; j < S->IA[i+1]; j++ ) {
1957 if ( (vec[k] ==
CGPT) && (visited[k] != i) ) {
1959 P->
JA[P->
IA[i]+index] = k;
1964 else if ( (vec[k] ==
FGPT) && (k != i) ) {
1965 for ( l = S->
IA[k]; l < S->IA[k+1]; l++ ) {
1967 if ( (vec[h] ==
CGPT) && (visited[h] != i) ) {
1969 P->
JA[P->
IA[i]+index] = h;
1980 else if ( vec[i] ==
CGPT ) {
1981 P->
JA[P->
IA[i]] = i;
2011 INT *vec = vertices->
val;
2016 INT row_begin, row_end;
2020 for (i=0; i<n ; i++)
2023 if (vec[ind] ==
UNPT) {
2025 row_begin = IS[ind]; row_end = IS[ind+1];
2026 for (j = row_begin; j<row_end; j++)
2028 if (vec[JS[j]] ==
CGPT ) {
2033 if (vec[ind] ==
CGPT) {
2035 for (j = row_begin; j<row_end; j++)
2059static void ordering1 (
iCSRmat *S,
2062 const INT n = order->
row;
2065 INT maxind, maxdeg, degree;
2068 for (i = 0; i < n; i++) ord[i] = i;
2070 for (maxind = maxdeg = i = 0; i < n; i++)
2072 degree = IS[i+1] - IS[i];
2073 if (degree > maxdeg)
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
void fasp_iarray_cp(const INT n, const INT *x, INT *y)
Copy an array to the other y=x.
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
void fasp_get_start_end(const INT procid, const INT nprocs, const INT n, INT *start, INT *end)
Assign Load to each thread.
void fasp_dvec_free(dvector *u)
Free vector data space of REAL type.
void fasp_ivec_free(ivector *u)
Free vector data space of INT type.
void fasp_ivec_alloc(const INT m, ivector *u)
Create vector data space of INT type.
ivector fasp_ivec_create(const INT m)
Create vector data space of INT type.
void fasp_ivec_set(INT n, ivector *u, const INT m)
Set ivector value to be m.
void fasp_icsr_trans(const iCSRmat *A, iCSRmat *AT)
Find transpose of iCSRmat matrix A.
void fasp_icsr_free(iCSRmat *A)
Free CSR sparse matrix data memory space.
void fasp_dcsr_getdiag(INT n, const dCSRmat *A, dvector *diag)
Get first n diagonal entries of a CSR matrix A.
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 SHORT
FASP integer and floating point numbers.
#define MAX(a, b)
Definition of max, min, abs.
#define ERROR_AMG_INTERP_TYPE
#define FASP_SUCCESS
Definition of return status and error messages.
#define INTERP_DIR
Definition of interpolation types.
#define TRUE
Definition of logic type.
#define ERROR_AMG_COARSEING
Parameters for AMG methods.
REAL strong_threshold
strong connection threshold for coarsening
INT aggressive_path
number of paths use to determine strongly coupled C points
SHORT interpolation_type
interpolation type
SHORT coarsening_type
coarsening type
REAL max_row_sum
maximal row sum parameter
Sparse matrix of REAL type in CSR format.
INT col
column of matrix A, n
REAL * val
nonzero entries of A
INT row
row number of matrix A, m
INT * IA
integer array of row pointers, the size is m+1
INT nnz
number of nonzero entries
INT * JA
integer array of column indexes, the size is nnz
Vector with n entries of REAL type.
REAL * val
actual vector entries
Sparse matrix of INT type in CSR format.
INT col
column of matrix A, n
INT row
row number of matrix A, m
INT * IA
integer array of row pointers, the size is m+1
INT nnz
number of nonzero entries
INT * JA
integer array of column indexes, the size is nnz
INT * val
nonzero entries of A
Vector with n entries of INT type.
INT * val
actual vector entries