1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2 #include <petscsf.h>
3
4 PETSC_EXTERN PetscErrorCode MatColoringCreateBipartiteGraph(MatColoring, PetscSF *, PetscSF *);
5
MatColoringTest(MatColoring mc,ISColoring coloring)6 PETSC_EXTERN PetscErrorCode MatColoringTest(MatColoring mc, ISColoring coloring)
7 {
8 Mat m = mc->mat;
9 PetscSF etor, etoc;
10 PetscInt s, e;
11 PetscInt ncolors, nrows, ncols;
12 IS *colors;
13 PetscInt i, j, k, l;
14 PetscInt *staterow, *statecol, *statespread;
15 PetscInt nindices;
16 const PetscInt *indices;
17 PetscInt dist = mc->dist;
18 const PetscInt *degrees;
19 PetscInt *stateleafrow, *stateleafcol, nleafrows, nleafcols, idx, nentries, maxcolors;
20 MPI_Datatype itype = MPIU_INT;
21
22 PetscFunctionBegin;
23 PetscCall(MatColoringGetMaxColors(mc, &maxcolors));
24 /* get the communication structures and the colors */
25 PetscCall(MatColoringCreateBipartiteGraph(mc, &etoc, &etor));
26 PetscCall(ISColoringGetIS(coloring, PETSC_USE_POINTER, &ncolors, &colors));
27 PetscCall(PetscSFGetGraph(etor, &nrows, &nleafrows, NULL, NULL));
28 PetscCall(PetscSFGetGraph(etoc, &ncols, &nleafcols, NULL, NULL));
29 PetscCall(MatGetOwnershipRangeColumn(m, &s, &e));
30 PetscCall(PetscMalloc1(ncols, &statecol));
31 PetscCall(PetscMalloc1(nrows, &staterow));
32 PetscCall(PetscMalloc1(nleafcols, &stateleafcol));
33 PetscCall(PetscMalloc1(nleafrows, &stateleafrow));
34
35 for (l = 0; l < ncolors; l++) {
36 if (l > maxcolors) break;
37 for (k = 0; k < ncols; k++) statecol[k] = -1;
38 PetscCall(ISGetLocalSize(colors[l], &nindices));
39 PetscCall(ISGetIndices(colors[l], &indices));
40 for (k = 0; k < nindices; k++) statecol[indices[k] - s] = indices[k];
41 PetscCall(ISRestoreIndices(colors[l], &indices));
42 statespread = statecol;
43 for (k = 0; k < dist; k++) {
44 if (k % 2 == 1) {
45 PetscCall(PetscSFComputeDegreeBegin(etor, °rees));
46 PetscCall(PetscSFComputeDegreeEnd(etor, °rees));
47 nentries = 0;
48 for (i = 0; i < nrows; i++) nentries += degrees[i];
49 idx = 0;
50 for (i = 0; i < nrows; i++) {
51 for (j = 0; j < degrees[i]; j++) {
52 stateleafrow[idx] = staterow[i];
53 idx++;
54 }
55 statecol[i] = 0.;
56 }
57 PetscCheck(idx == nentries, PetscObjectComm((PetscObject)mc), PETSC_ERR_NOT_CONVERGED, "Bad number of entries %" PetscInt_FMT " vs %" PetscInt_FMT, idx, nentries);
58 PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
59 PetscCall(PetscSFReduceBegin(etoc, itype, stateleafrow, statecol, MPI_MAX));
60 PetscCall(PetscSFReduceEnd(etoc, itype, stateleafrow, statecol, MPI_MAX));
61 PetscCall(PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0));
62 statespread = statecol;
63 } else {
64 PetscCall(PetscSFComputeDegreeBegin(etoc, °rees));
65 PetscCall(PetscSFComputeDegreeEnd(etoc, °rees));
66 nentries = 0;
67 for (i = 0; i < ncols; i++) nentries += degrees[i];
68 idx = 0;
69 for (i = 0; i < ncols; i++) {
70 for (j = 0; j < degrees[i]; j++) {
71 stateleafcol[idx] = statecol[i];
72 idx++;
73 }
74 staterow[i] = 0.;
75 }
76 PetscCheck(idx == nentries, PetscObjectComm((PetscObject)mc), PETSC_ERR_NOT_CONVERGED, "Bad number of entries %" PetscInt_FMT " vs %" PetscInt_FMT, idx, nentries);
77 PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
78 PetscCall(PetscSFReduceBegin(etor, itype, stateleafcol, staterow, MPI_MAX));
79 PetscCall(PetscSFReduceEnd(etor, itype, stateleafcol, staterow, MPI_MAX));
80 PetscCall(PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0));
81 statespread = staterow;
82 }
83 }
84 for (k = 0; k < nindices; k++) {
85 if (statespread[indices[k] - s] != indices[k]) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mc), "%" PetscInt_FMT " of color %" PetscInt_FMT " conflicts with %" PetscInt_FMT "\n", indices[k], l, statespread[indices[k] - s]));
86 }
87 PetscCall(ISRestoreIndices(colors[l], &indices));
88 }
89 PetscCall(PetscFree(statecol));
90 PetscCall(PetscFree(staterow));
91 PetscCall(PetscFree(stateleafcol));
92 PetscCall(PetscFree(stateleafrow));
93 PetscCall(PetscSFDestroy(&etor));
94 PetscCall(PetscSFDestroy(&etoc));
95 PetscFunctionReturn(PETSC_SUCCESS);
96 }
97
MatISColoringTest(Mat A,ISColoring iscoloring)98 PETSC_EXTERN PetscErrorCode MatISColoringTest(Mat A, ISColoring iscoloring)
99 {
100 PetscInt nn, c, i, j, M, N, nc, nnz, col, row;
101 const PetscInt *cia, *cja, *cols;
102 IS *isis;
103 MPI_Comm comm;
104 PetscMPIInt size;
105 PetscBool done;
106 PetscBT table;
107
108 PetscFunctionBegin;
109 PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &nn, &isis));
110
111 PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
112 PetscCallMPI(MPI_Comm_size(comm, &size));
113 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support sequential matrix");
114
115 PetscCall(MatGetColumnIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &N, &cia, &cja, &done));
116 PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ordering requires IJ");
117
118 PetscCall(MatGetSize(A, &M, NULL));
119 PetscCall(PetscBTCreate(M, &table));
120 for (c = 0; c < nn; c++) { /* for each color */
121 PetscCall(ISGetSize(isis[c], &nc));
122 if (nc <= 1) continue;
123
124 PetscCall(PetscBTMemzero(M, table));
125 PetscCall(ISGetIndices(isis[c], &cols));
126 for (j = 0; j < nc; j++) { /* for each column */
127 col = cols[j];
128 nnz = cia[col + 1] - cia[col];
129 for (i = 0; i < nnz; i++) {
130 row = cja[cia[col] + i];
131 PetscCheck(!PetscBTLookupSet(table, row), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "color %" PetscInt_FMT ", col %" PetscInt_FMT ": row %" PetscInt_FMT " already in this color", c, col, row);
132 }
133 }
134 PetscCall(ISRestoreIndices(isis[c], &cols));
135 }
136 PetscCall(PetscBTDestroy(&table));
137
138 PetscCall(MatRestoreColumnIJ(A, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done));
139 PetscFunctionReturn(PETSC_SUCCESS);
140 }
141