Actual source code: color.c
1: /*
2: Routines that call the kernel minpack coloring subroutines
3: */
5: #include <petsc/private/matimpl.h>
6: #include <petsc/private/isimpl.h>
7: #include <../src/mat/graphops/color/impls/minpack/color.h>
9: /*
10: MatFDColoringDegreeSequence_Minpack - Calls the MINPACK routine seqr() that
11: computes the degree sequence required by MINPACK coloring routines.
12: */
13: PETSC_INTERN PetscErrorCode MatFDColoringDegreeSequence_Minpack(PetscInt m, const PetscInt *cja, const PetscInt *cia, const PetscInt *rja, const PetscInt *ria, PetscInt **seq)
14: {
15: PetscInt *work;
17: PetscFunctionBegin;
18: PetscCall(PetscMalloc1(m, &work));
19: PetscCall(PetscMalloc1(m, seq));
21: PetscCall(MINPACKdegr(&m, cja, cia, rja, ria, *seq, work));
23: PetscCall(PetscFree(work));
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode MatColoringApply_SL(MatColoring mc, ISColoring *iscoloring)
28: {
29: PetscInt *list, *work, clique, *seq, *coloring, n;
30: const PetscInt *ria, *rja, *cia, *cja;
31: PetscInt ncolors, i;
32: PetscBool done;
33: Mat mat = mc->mat;
34: Mat mat_seq = mat;
35: PetscMPIInt size;
36: MPI_Comm comm;
37: ISColoring iscoloring_seq;
38: PetscInt bs = 1, rstart, rend, N_loc, nc;
39: ISColoringValue *colors_loc;
40: PetscBool flg1, flg2;
42: PetscFunctionBegin;
43: PetscCheck(mc->dist == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "SL may only do distance 2 coloring");
44: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
45: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
46: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
47: if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));
49: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
50: PetscCallMPI(MPI_Comm_size(comm, &size));
51: if (size > 1) {
52: /* create a sequential iscoloring on all processors */
53: PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
54: }
56: PetscCall(MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done));
57: PetscCall(MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done));
58: PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ordering requires IJ");
60: PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
62: PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
64: PetscCall(MINPACKslo(&n, cja, cia, rja, ria, seq, list, &clique, work, work + n, work + 2 * n, work + 3 * n));
66: PetscCall(PetscMalloc1(n, &coloring));
67: PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
69: PetscCall(PetscFree2(list, work));
70: PetscCall(PetscFree(seq));
71: PetscCall(MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done));
72: PetscCall(MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done));
74: /* shift coloring numbers to start at zero and shorten */
75: PetscCheck(ncolors <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");
76: {
77: ISColoringValue *s = (ISColoringValue *)coloring;
78: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
79: PetscCall(MatColoringPatch(mat_seq, ncolors, n, s, iscoloring));
80: }
82: if (size > 1) {
83: PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
85: /* convert iscoloring_seq to a parallel iscoloring */
86: iscoloring_seq = *iscoloring;
87: rstart = mat->rmap->rstart / bs;
88: rend = mat->rmap->rend / bs;
89: N_loc = rend - rstart; /* number of local nodes */
91: /* get local colors for each local node */
92: PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
93: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
94: /* create a parallel iscoloring */
95: nc = iscoloring_seq->n;
96: PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
97: PetscCall(ISColoringDestroy(&iscoloring_seq));
98: }
99: PetscFunctionReturn(PETSC_SUCCESS);
100: }
102: /*MC
103: MATCOLORINGSL - implements the SL (smallest last) coloring routine {cite}`more:coloring`
105: Level: beginner
107: Notes:
108: Supports only distance two colorings (for computation of Jacobians/Hessians)
110: This is a sequential algorithm
112: .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
113: M*/
115: PETSC_EXTERN PetscErrorCode MatColoringCreate_SL(MatColoring mc)
116: {
117: PetscFunctionBegin;
118: mc->dist = 2;
119: mc->data = NULL;
120: mc->ops->apply = MatColoringApply_SL;
121: mc->ops->view = NULL;
122: mc->ops->destroy = NULL;
123: mc->ops->setfromoptions = NULL;
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode MatColoringApply_LF(MatColoring mc, ISColoring *iscoloring)
128: {
129: PetscInt *list, *work, *seq, *coloring, n;
130: const PetscInt *ria, *rja, *cia, *cja;
131: PetscInt n1, none, ncolors, i;
132: PetscBool done;
133: Mat mat = mc->mat;
134: Mat mat_seq = mat;
135: PetscMPIInt size;
136: MPI_Comm comm;
137: ISColoring iscoloring_seq;
138: PetscInt bs = 1, rstart, rend, N_loc, nc;
139: ISColoringValue *colors_loc;
140: PetscBool flg1, flg2;
142: PetscFunctionBegin;
143: PetscCheck(mc->dist == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "LF may only do distance 2 coloring");
144: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
145: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
146: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
147: if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));
149: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
150: PetscCallMPI(MPI_Comm_size(comm, &size));
151: if (size > 1) {
152: /* create a sequential iscoloring on all processors */
153: PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
154: }
156: PetscCall(MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done));
157: PetscCall(MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done));
158: PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ordering requires IJ");
160: PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
162: PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
164: n1 = n - 1;
165: none = -1;
166: PetscCall(MINPACKnumsrt(&n, &n1, seq, &none, list, work + 2 * n, work + n));
167: PetscCall(PetscMalloc1(n, &coloring));
168: PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
170: PetscCall(PetscFree2(list, work));
171: PetscCall(PetscFree(seq));
173: PetscCall(MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done));
174: PetscCall(MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done));
176: /* shift coloring numbers to start at zero and shorten */
177: PetscCheck(ncolors <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");
178: {
179: ISColoringValue *s = (ISColoringValue *)coloring;
180: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
181: PetscCall(MatColoringPatch(mat_seq, ncolors, n, s, iscoloring));
182: }
184: if (size > 1) {
185: PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
187: /* convert iscoloring_seq to a parallel iscoloring */
188: iscoloring_seq = *iscoloring;
189: rstart = mat->rmap->rstart / bs;
190: rend = mat->rmap->rend / bs;
191: N_loc = rend - rstart; /* number of local nodes */
193: /* get local colors for each local node */
194: PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
195: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
197: /* create a parallel iscoloring */
198: nc = iscoloring_seq->n;
199: PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
200: PetscCall(ISColoringDestroy(&iscoloring_seq));
201: }
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: /*MC
206: MATCOLORINGLF - implements the LF (largest first) coloring routine {cite}`more:coloring`
208: Level: beginner
210: Notes:
211: Supports only distance two colorings (for computation of Jacobians/Hessians)
213: This is a sequential algorithm
215: .seealso: `MatColoringTpe`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
216: M*/
218: PETSC_EXTERN PetscErrorCode MatColoringCreate_LF(MatColoring mc)
219: {
220: PetscFunctionBegin;
221: mc->dist = 2;
222: mc->data = NULL;
223: mc->ops->apply = MatColoringApply_LF;
224: mc->ops->view = NULL;
225: mc->ops->destroy = NULL;
226: mc->ops->setfromoptions = NULL;
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: static PetscErrorCode MatColoringApply_ID(MatColoring mc, ISColoring *iscoloring)
231: {
232: PetscInt *list, *work, clique, *seq, *coloring, n;
233: const PetscInt *ria, *rja, *cia, *cja;
234: PetscInt ncolors, i;
235: PetscBool done;
236: Mat mat = mc->mat;
237: Mat mat_seq = mat;
238: PetscMPIInt size;
239: MPI_Comm comm;
240: ISColoring iscoloring_seq;
241: PetscInt bs = 1, rstart, rend, N_loc, nc;
242: ISColoringValue *colors_loc;
243: PetscBool flg1, flg2;
245: PetscFunctionBegin;
246: PetscCheck(mc->dist == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "IDO may only do distance 2 coloring");
247: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
248: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
249: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
250: if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));
252: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
253: PetscCallMPI(MPI_Comm_size(comm, &size));
254: if (size > 1) {
255: /* create a sequential iscoloring on all processors */
256: PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
257: }
259: PetscCall(MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done));
260: PetscCall(MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done));
261: PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ordering requires IJ");
263: PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
265: PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
267: PetscCall(MINPACKido(&n, &n, cja, cia, rja, ria, seq, list, &clique, work, work + n, work + 2 * n, work + 3 * n));
269: PetscCall(PetscMalloc1(n, &coloring));
270: PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
272: PetscCall(PetscFree2(list, work));
273: PetscCall(PetscFree(seq));
275: PetscCall(MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done));
276: PetscCall(MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done));
278: /* shift coloring numbers to start at zero and shorten */
279: PetscCheck(ncolors <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");
280: {
281: ISColoringValue *s = (ISColoringValue *)coloring;
282: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
283: PetscCall(MatColoringPatch(mat_seq, ncolors, n, s, iscoloring));
284: }
286: if (size > 1) {
287: PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
289: /* convert iscoloring_seq to a parallel iscoloring */
290: iscoloring_seq = *iscoloring;
291: rstart = mat->rmap->rstart / bs;
292: rend = mat->rmap->rend / bs;
293: N_loc = rend - rstart; /* number of local nodes */
295: /* get local colors for each local node */
296: PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
297: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
298: /* create a parallel iscoloring */
299: nc = iscoloring_seq->n;
300: PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
301: PetscCall(ISColoringDestroy(&iscoloring_seq));
302: }
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*MC
307: MATCOLORINGID - implements the ID (incidence degree) coloring routine {cite}`more:coloring`
309: Level: beginner
311: Notes:
312: Supports only distance two colorings (for computation of Jacobians/Hessians)
314: This is a sequential algorithm
316: .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
317: M*/
319: PETSC_EXTERN PetscErrorCode MatColoringCreate_ID(MatColoring mc)
320: {
321: PetscFunctionBegin;
322: mc->dist = 2;
323: mc->data = NULL;
324: mc->ops->apply = MatColoringApply_ID;
325: mc->ops->view = NULL;
326: mc->ops->destroy = NULL;
327: mc->ops->setfromoptions = NULL;
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }