Actual source code: color.c
2: /*
3: Routines that call the kernel minpack coloring subroutines
4: */
6: #include <petsc/private/matimpl.h>
7: #include <petsc/private/isimpl.h>
8: #include <../src/mat/color/impls/minpack/color.h>
10: /*
11: MatFDColoringDegreeSequence_Minpack - Calls the MINPACK routine seqr() that
12: computes the degree sequence required by MINPACK coloring routines.
13: */
14: PETSC_INTERN PetscErrorCode MatFDColoringDegreeSequence_Minpack(PetscInt m, const PetscInt *cja, const PetscInt *cia, const PetscInt *rja, const PetscInt *ria, PetscInt **seq)
15: {
16: PetscInt *work;
18: PetscFunctionBegin;
19: PetscCall(PetscMalloc1(m, &work));
20: PetscCall(PetscMalloc1(m, seq));
22: PetscCall(MINPACKdegr(&m, cja, cia, rja, ria, *seq, work));
24: PetscCall(PetscFree(work));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: /*
29: MatFDColoringMinimumNumberofColors_Private - For a given sparse
30: matrix computes the minimum number of colors needed.
32: */
33: PetscErrorCode MatFDColoringMinimumNumberofColors_Private(PetscInt m, PetscInt *ia, PetscInt *minc)
34: {
35: PetscInt i, c = 0;
37: PetscFunctionBegin;
38: for (i = 0; i < m; i++) c = PetscMax(c, ia[i + 1] - ia[i]);
39: *minc = c;
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode MatColoringApply_SL(MatColoring mc, ISColoring *iscoloring)
44: {
45: PetscInt *list, *work, clique, *seq, *coloring, n;
46: const PetscInt *ria, *rja, *cia, *cja;
47: PetscInt ncolors, i;
48: PetscBool done;
49: Mat mat = mc->mat;
50: Mat mat_seq = mat;
51: PetscMPIInt size;
52: MPI_Comm comm;
53: ISColoring iscoloring_seq;
54: PetscInt bs = 1, rstart, rend, N_loc, nc;
55: ISColoringValue *colors_loc;
56: PetscBool flg1, flg2;
58: PetscFunctionBegin;
59: PetscCheck(mc->dist == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "SL may only do distance 2 coloring");
60: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
61: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
62: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
63: if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));
65: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
66: PetscCallMPI(MPI_Comm_size(comm, &size));
67: if (size > 1) {
68: /* create a sequential iscoloring on all processors */
69: PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
70: }
72: PetscCall(MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done));
73: PetscCall(MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done));
74: PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ordering requires IJ");
76: PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
78: PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
80: PetscCall(MINPACKslo(&n, cja, cia, rja, ria, seq, list, &clique, work, work + n, work + 2 * n, work + 3 * n));
82: PetscCall(PetscMalloc1(n, &coloring));
83: PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
85: PetscCall(PetscFree2(list, work));
86: PetscCall(PetscFree(seq));
87: PetscCall(MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done));
88: PetscCall(MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done));
90: /* shift coloring numbers to start at zero and shorten */
91: PetscCheck(ncolors <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");
92: {
93: ISColoringValue *s = (ISColoringValue *)coloring;
94: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
95: PetscCall(MatColoringPatch(mat_seq, ncolors, n, s, iscoloring));
96: }
98: if (size > 1) {
99: PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
101: /* convert iscoloring_seq to a parallel iscoloring */
102: iscoloring_seq = *iscoloring;
103: rstart = mat->rmap->rstart / bs;
104: rend = mat->rmap->rend / bs;
105: N_loc = rend - rstart; /* number of local nodes */
107: /* get local colors for each local node */
108: PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
109: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
110: /* create a parallel iscoloring */
111: nc = iscoloring_seq->n;
112: PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
113: PetscCall(ISColoringDestroy(&iscoloring_seq));
114: }
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: /*MC
119: MATCOLORINGSL - implements the SL (smallest last) coloring routine
121: Level: beginner
123: Notes:
124: Supports only distance two colorings (for computation of Jacobians)
126: This is a sequential algorithm
128: References:
129: . * - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
130: pp. 187-209, 1983.
132: .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
133: M*/
135: PETSC_EXTERN PetscErrorCode MatColoringCreate_SL(MatColoring mc)
136: {
137: PetscFunctionBegin;
138: mc->dist = 2;
139: mc->data = NULL;
140: mc->ops->apply = MatColoringApply_SL;
141: mc->ops->view = NULL;
142: mc->ops->destroy = NULL;
143: mc->ops->setfromoptions = NULL;
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: static PetscErrorCode MatColoringApply_LF(MatColoring mc, ISColoring *iscoloring)
148: {
149: PetscInt *list, *work, *seq, *coloring, n;
150: const PetscInt *ria, *rja, *cia, *cja;
151: PetscInt n1, none, ncolors, i;
152: PetscBool done;
153: Mat mat = mc->mat;
154: Mat mat_seq = mat;
155: PetscMPIInt size;
156: MPI_Comm comm;
157: ISColoring iscoloring_seq;
158: PetscInt bs = 1, rstart, rend, N_loc, nc;
159: ISColoringValue *colors_loc;
160: PetscBool flg1, flg2;
162: PetscFunctionBegin;
163: PetscCheck(mc->dist == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "LF may only do distance 2 coloring");
164: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
165: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
166: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
167: if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));
169: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
170: PetscCallMPI(MPI_Comm_size(comm, &size));
171: if (size > 1) {
172: /* create a sequential iscoloring on all processors */
173: PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
174: }
176: PetscCall(MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done));
177: PetscCall(MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done));
178: PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ordering requires IJ");
180: PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
182: PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
184: n1 = n - 1;
185: none = -1;
186: PetscCall(MINPACKnumsrt(&n, &n1, seq, &none, list, work + 2 * n, work + n));
187: PetscCall(PetscMalloc1(n, &coloring));
188: PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
190: PetscCall(PetscFree2(list, work));
191: PetscCall(PetscFree(seq));
193: PetscCall(MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done));
194: PetscCall(MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done));
196: /* shift coloring numbers to start at zero and shorten */
197: PetscCheck(ncolors <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");
198: {
199: ISColoringValue *s = (ISColoringValue *)coloring;
200: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
201: PetscCall(MatColoringPatch(mat_seq, ncolors, n, s, iscoloring));
202: }
204: if (size > 1) {
205: PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
207: /* convert iscoloring_seq to a parallel iscoloring */
208: iscoloring_seq = *iscoloring;
209: rstart = mat->rmap->rstart / bs;
210: rend = mat->rmap->rend / bs;
211: N_loc = rend - rstart; /* number of local nodes */
213: /* get local colors for each local node */
214: PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
215: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
217: /* create a parallel iscoloring */
218: nc = iscoloring_seq->n;
219: PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
220: PetscCall(ISColoringDestroy(&iscoloring_seq));
221: }
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: /*MC
226: MATCOLORINGLF - implements the LF (largest first) coloring routine
228: Level: beginner
230: Notes:
231: Supports only distance two colorings (for computation of Jacobians)
233: This is a sequential algorithm
235: References:
236: . * - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
237: pp. 187-209, 1983.
239: .seealso: `MatColoringTpe`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
240: M*/
242: PETSC_EXTERN PetscErrorCode MatColoringCreate_LF(MatColoring mc)
243: {
244: PetscFunctionBegin;
245: mc->dist = 2;
246: mc->data = NULL;
247: mc->ops->apply = MatColoringApply_LF;
248: mc->ops->view = NULL;
249: mc->ops->destroy = NULL;
250: mc->ops->setfromoptions = NULL;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: static PetscErrorCode MatColoringApply_ID(MatColoring mc, ISColoring *iscoloring)
255: {
256: PetscInt *list, *work, clique, *seq, *coloring, n;
257: const PetscInt *ria, *rja, *cia, *cja;
258: PetscInt ncolors, i;
259: PetscBool done;
260: Mat mat = mc->mat;
261: Mat mat_seq = mat;
262: PetscMPIInt size;
263: MPI_Comm comm;
264: ISColoring iscoloring_seq;
265: PetscInt bs = 1, rstart, rend, N_loc, nc;
266: ISColoringValue *colors_loc;
267: PetscBool flg1, flg2;
269: PetscFunctionBegin;
270: PetscCheck(mc->dist == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "IDO may only do distance 2 coloring");
271: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
272: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
273: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
274: if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));
276: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
277: PetscCallMPI(MPI_Comm_size(comm, &size));
278: if (size > 1) {
279: /* create a sequential iscoloring on all processors */
280: PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
281: }
283: PetscCall(MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done));
284: PetscCall(MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done));
285: PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ordering requires IJ");
287: PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
289: PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
291: PetscCall(MINPACKido(&n, &n, cja, cia, rja, ria, seq, list, &clique, work, work + n, work + 2 * n, work + 3 * n));
293: PetscCall(PetscMalloc1(n, &coloring));
294: PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
296: PetscCall(PetscFree2(list, work));
297: PetscCall(PetscFree(seq));
299: PetscCall(MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done));
300: PetscCall(MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done));
302: /* shift coloring numbers to start at zero and shorten */
303: PetscCheck(ncolors <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");
304: {
305: ISColoringValue *s = (ISColoringValue *)coloring;
306: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
307: PetscCall(MatColoringPatch(mat_seq, ncolors, n, s, iscoloring));
308: }
310: if (size > 1) {
311: PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
313: /* convert iscoloring_seq to a parallel iscoloring */
314: iscoloring_seq = *iscoloring;
315: rstart = mat->rmap->rstart / bs;
316: rend = mat->rmap->rend / bs;
317: N_loc = rend - rstart; /* number of local nodes */
319: /* get local colors for each local node */
320: PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
321: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
322: /* create a parallel iscoloring */
323: nc = iscoloring_seq->n;
324: PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
325: PetscCall(ISColoringDestroy(&iscoloring_seq));
326: }
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: /*MC
331: MATCOLORINGID - implements the ID (incidence degree) coloring routine
333: Level: beginner
335: Notes:
336: Supports only distance two colorings (for computation of Jacobians)
338: This is a sequential algorithm
340: References:
341: . * - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
342: pp. 187-209, 1983.
344: .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
345: M*/
347: PETSC_EXTERN PetscErrorCode MatColoringCreate_ID(MatColoring mc)
348: {
349: PetscFunctionBegin;
350: mc->dist = 2;
351: mc->data = NULL;
352: mc->ops->apply = MatColoringApply_ID;
353: mc->ops->view = NULL;
354: mc->ops->destroy = NULL;
355: mc->ops->setfromoptions = NULL;
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }