Actual source code: valid.c
1: #include <petsc/private/matimpl.h>
2: #include <petscsf.h>
4: PETSC_EXTERN PetscErrorCode MatColoringCreateBipartiteGraph(MatColoring, PetscSF *, PetscSF *);
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;
22: MatColoringGetMaxColors(mc, &maxcolors);
23: /* get the communication structures and the colors */
24: MatColoringCreateBipartiteGraph(mc, &etoc, &etor);
25: ISColoringGetIS(coloring, PETSC_USE_POINTER, &ncolors, &colors);
26: PetscSFGetGraph(etor, &nrows, &nleafrows, NULL, NULL);
27: PetscSFGetGraph(etoc, &ncols, &nleafcols, NULL, NULL);
28: MatGetOwnershipRangeColumn(m, &s, &e);
29: PetscMalloc1(ncols, &statecol);
30: PetscMalloc1(nrows, &staterow);
31: PetscMalloc1(nleafcols, &stateleafcol);
32: PetscMalloc1(nleafrows, &stateleafrow);
34: for (l = 0; l < ncolors; l++) {
35: if (l > maxcolors) break;
36: for (k = 0; k < ncols; k++) statecol[k] = -1;
37: ISGetLocalSize(colors[l], &nindices);
38: ISGetIndices(colors[l], &indices);
39: for (k = 0; k < nindices; k++) statecol[indices[k] - s] = indices[k];
40: ISRestoreIndices(colors[l], &indices);
41: statespread = statecol;
42: for (k = 0; k < dist; k++) {
43: if (k % 2 == 1) {
44: PetscSFComputeDegreeBegin(etor, °rees);
45: PetscSFComputeDegreeEnd(etor, °rees);
46: nentries = 0;
47: for (i = 0; i < nrows; i++) nentries += degrees[i];
48: idx = 0;
49: for (i = 0; i < nrows; i++) {
50: for (j = 0; j < degrees[i]; j++) {
51: stateleafrow[idx] = staterow[i];
52: idx++;
53: }
54: statecol[i] = 0.;
55: }
57: PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0);
58: PetscSFReduceBegin(etoc, itype, stateleafrow, statecol, MPI_MAX);
59: PetscSFReduceEnd(etoc, itype, stateleafrow, statecol, MPI_MAX);
60: PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0);
61: statespread = statecol;
62: } else {
63: PetscSFComputeDegreeBegin(etoc, °rees);
64: PetscSFComputeDegreeEnd(etoc, °rees);
65: nentries = 0;
66: for (i = 0; i < ncols; i++) nentries += degrees[i];
67: idx = 0;
68: for (i = 0; i < ncols; i++) {
69: for (j = 0; j < degrees[i]; j++) {
70: stateleafcol[idx] = statecol[i];
71: idx++;
72: }
73: staterow[i] = 0.;
74: }
76: PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0);
77: PetscSFReduceBegin(etor, itype, stateleafcol, staterow, MPI_MAX);
78: PetscSFReduceEnd(etor, itype, stateleafcol, staterow, MPI_MAX);
79: PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0);
80: statespread = staterow;
81: }
82: }
83: for (k = 0; k < nindices; k++) {
84: if (statespread[indices[k] - s] != indices[k]) PetscPrintf(PetscObjectComm((PetscObject)mc), "%" PetscInt_FMT " of color %" PetscInt_FMT " conflicts with %" PetscInt_FMT "\n", indices[k], l, statespread[indices[k] - s]);
85: }
86: ISRestoreIndices(colors[l], &indices);
87: }
88: PetscFree(statecol);
89: PetscFree(staterow);
90: PetscFree(stateleafcol);
91: PetscFree(stateleafrow);
92: PetscSFDestroy(&etor);
93: PetscSFDestroy(&etoc);
94: return 0;
95: }
97: PETSC_EXTERN PetscErrorCode MatISColoringTest(Mat A, ISColoring iscoloring)
98: {
99: PetscInt nn, c, i, j, M, N, nc, nnz, col, row;
100: const PetscInt *cia, *cja, *cols;
101: IS *isis;
102: MPI_Comm comm;
103: PetscMPIInt size;
104: PetscBool done;
105: PetscBT table;
107: ISColoringGetIS(iscoloring, PETSC_USE_POINTER, &nn, &isis);
109: PetscObjectGetComm((PetscObject)A, &comm);
110: MPI_Comm_size(comm, &size);
113: MatGetColumnIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &N, &cia, &cja, &done);
116: MatGetSize(A, &M, NULL);
117: PetscBTCreate(M, &table);
118: for (c = 0; c < nn; c++) { /* for each color */
119: ISGetSize(isis[c], &nc);
120: if (nc <= 1) continue;
122: PetscBTMemzero(M, table);
123: ISGetIndices(isis[c], &cols);
124: for (j = 0; j < nc; j++) { /* for each column */
125: col = cols[j];
126: nnz = cia[col + 1] - cia[col];
127: for (i = 0; i < nnz; i++) {
128: row = cja[cia[col] + i];
130: }
131: }
132: ISRestoreIndices(isis[c], &cols);
133: }
134: PetscBTDestroy(&table);
136: MatRestoreColumnIJ(A, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done);
137: return 0;
138: }