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/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
105: Level: beginner
107: Notes:
108: Supports only distance two colorings (for computation of Jacobians)
110: This is a sequential algorithm
112: References:
113: . * - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
114: pp. 187-209, 1983.
116: .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
117: M*/
119: PETSC_EXTERN PetscErrorCode MatColoringCreate_SL(MatColoring mc)
120: {
121: PetscFunctionBegin;
122: mc->dist = 2;
123: mc->data = NULL;
124: mc->ops->apply = MatColoringApply_SL;
125: mc->ops->view = NULL;
126: mc->ops->destroy = NULL;
127: mc->ops->setfromoptions = NULL;
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode MatColoringApply_LF(MatColoring mc, ISColoring *iscoloring)
132: {
133: PetscInt *list, *work, *seq, *coloring, n;
134: const PetscInt *ria, *rja, *cia, *cja;
135: PetscInt n1, none, ncolors, i;
136: PetscBool done;
137: Mat mat = mc->mat;
138: Mat mat_seq = mat;
139: PetscMPIInt size;
140: MPI_Comm comm;
141: ISColoring iscoloring_seq;
142: PetscInt bs = 1, rstart, rend, N_loc, nc;
143: ISColoringValue *colors_loc;
144: PetscBool flg1, flg2;
146: PetscFunctionBegin;
147: PetscCheck(mc->dist == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "LF may only do distance 2 coloring");
148: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
149: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
150: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
151: if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));
153: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
154: PetscCallMPI(MPI_Comm_size(comm, &size));
155: if (size > 1) {
156: /* create a sequential iscoloring on all processors */
157: PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
158: }
160: PetscCall(MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done));
161: PetscCall(MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done));
162: PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ordering requires IJ");
164: PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
166: PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
168: n1 = n - 1;
169: none = -1;
170: PetscCall(MINPACKnumsrt(&n, &n1, seq, &none, list, work + 2 * n, work + n));
171: PetscCall(PetscMalloc1(n, &coloring));
172: PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
174: PetscCall(PetscFree2(list, work));
175: PetscCall(PetscFree(seq));
177: PetscCall(MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done));
178: PetscCall(MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done));
180: /* shift coloring numbers to start at zero and shorten */
181: PetscCheck(ncolors <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");
182: {
183: ISColoringValue *s = (ISColoringValue *)coloring;
184: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
185: PetscCall(MatColoringPatch(mat_seq, ncolors, n, s, iscoloring));
186: }
188: if (size > 1) {
189: PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
191: /* convert iscoloring_seq to a parallel iscoloring */
192: iscoloring_seq = *iscoloring;
193: rstart = mat->rmap->rstart / bs;
194: rend = mat->rmap->rend / bs;
195: N_loc = rend - rstart; /* number of local nodes */
197: /* get local colors for each local node */
198: PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
199: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
201: /* create a parallel iscoloring */
202: nc = iscoloring_seq->n;
203: PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
204: PetscCall(ISColoringDestroy(&iscoloring_seq));
205: }
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /*MC
210: MATCOLORINGLF - implements the LF (largest first) coloring routine
212: Level: beginner
214: Notes:
215: Supports only distance two colorings (for computation of Jacobians)
217: This is a sequential algorithm
219: References:
220: . * - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
221: pp. 187-209, 1983.
223: .seealso: `MatColoringTpe`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
224: M*/
226: PETSC_EXTERN PetscErrorCode MatColoringCreate_LF(MatColoring mc)
227: {
228: PetscFunctionBegin;
229: mc->dist = 2;
230: mc->data = NULL;
231: mc->ops->apply = MatColoringApply_LF;
232: mc->ops->view = NULL;
233: mc->ops->destroy = NULL;
234: mc->ops->setfromoptions = NULL;
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: static PetscErrorCode MatColoringApply_ID(MatColoring mc, ISColoring *iscoloring)
239: {
240: PetscInt *list, *work, clique, *seq, *coloring, n;
241: const PetscInt *ria, *rja, *cia, *cja;
242: PetscInt ncolors, i;
243: PetscBool done;
244: Mat mat = mc->mat;
245: Mat mat_seq = mat;
246: PetscMPIInt size;
247: MPI_Comm comm;
248: ISColoring iscoloring_seq;
249: PetscInt bs = 1, rstart, rend, N_loc, nc;
250: ISColoringValue *colors_loc;
251: PetscBool flg1, flg2;
253: PetscFunctionBegin;
254: PetscCheck(mc->dist == 2, PETSC_COMM_SELF, PETSC_ERR_SUP, "IDO may only do distance 2 coloring");
255: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
256: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1));
257: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2));
258: if (flg1 || flg2) PetscCall(MatGetBlockSize(mat, &bs));
260: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
261: PetscCallMPI(MPI_Comm_size(comm, &size));
262: if (size > 1) {
263: /* create a sequential iscoloring on all processors */
264: PetscCall(MatGetSeqNonzeroStructure(mat, &mat_seq));
265: }
267: PetscCall(MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done));
268: PetscCall(MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done));
269: PetscCheck(done, PETSC_COMM_SELF, PETSC_ERR_SUP, "Ordering requires IJ");
271: PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
273: PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
275: PetscCall(MINPACKido(&n, &n, cja, cia, rja, ria, seq, list, &clique, work, work + n, work + 2 * n, work + 3 * n));
277: PetscCall(PetscMalloc1(n, &coloring));
278: PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
280: PetscCall(PetscFree2(list, work));
281: PetscCall(PetscFree(seq));
283: PetscCall(MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done));
284: PetscCall(MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done));
286: /* shift coloring numbers to start at zero and shorten */
287: PetscCheck(ncolors <= IS_COLORING_MAX - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Maximum color size exceeded");
288: {
289: ISColoringValue *s = (ISColoringValue *)coloring;
290: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
291: PetscCall(MatColoringPatch(mat_seq, ncolors, n, s, iscoloring));
292: }
294: if (size > 1) {
295: PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
297: /* convert iscoloring_seq to a parallel iscoloring */
298: iscoloring_seq = *iscoloring;
299: rstart = mat->rmap->rstart / bs;
300: rend = mat->rmap->rend / bs;
301: N_loc = rend - rstart; /* number of local nodes */
303: /* get local colors for each local node */
304: PetscCall(PetscMalloc1(N_loc + 1, &colors_loc));
305: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
306: /* create a parallel iscoloring */
307: nc = iscoloring_seq->n;
308: PetscCall(ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring));
309: PetscCall(ISColoringDestroy(&iscoloring_seq));
310: }
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*MC
315: MATCOLORINGID - implements the ID (incidence degree) coloring routine
317: Level: beginner
319: Notes:
320: Supports only distance two colorings (for computation of Jacobians)
322: This is a sequential algorithm
324: References:
325: . * - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
326: pp. 187-209, 1983.
328: .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
329: M*/
331: PETSC_EXTERN PetscErrorCode MatColoringCreate_ID(MatColoring mc)
332: {
333: PetscFunctionBegin;
334: mc->dist = 2;
335: mc->data = NULL;
336: mc->ops->apply = MatColoringApply_ID;
337: mc->ops->view = NULL;
338: mc->ops->destroy = NULL;
339: mc->ops->setfromoptions = NULL;
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }