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: PetscMalloc1(m, &work);
19: PetscMalloc1(m, seq);
21: MINPACKdegr(&m, cja, cia, rja, ria, *seq, work);
23: PetscFree(work);
24: return 0;
25: }
27: /*
28: MatFDColoringMinimumNumberofColors_Private - For a given sparse
29: matrix computes the minimum number of colors needed.
31: */
32: PetscErrorCode MatFDColoringMinimumNumberofColors_Private(PetscInt m, PetscInt *ia, PetscInt *minc)
33: {
34: PetscInt i, c = 0;
36: for (i = 0; i < m; i++) c = PetscMax(c, ia[i + 1] - ia[i]);
37: *minc = c;
38: return 0;
39: }
41: static PetscErrorCode MatColoringApply_SL(MatColoring mc, ISColoring *iscoloring)
42: {
43: PetscInt *list, *work, clique, *seq, *coloring, n;
44: const PetscInt *ria, *rja, *cia, *cja;
45: PetscInt ncolors, i;
46: PetscBool done;
47: Mat mat = mc->mat;
48: Mat mat_seq = mat;
49: PetscMPIInt size;
50: MPI_Comm comm;
51: ISColoring iscoloring_seq;
52: PetscInt bs = 1, rstart, rend, N_loc, nc;
53: ISColoringValue *colors_loc;
54: PetscBool flg1, flg2;
57: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
58: PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1);
59: PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2);
60: if (flg1 || flg2) MatGetBlockSize(mat, &bs);
62: PetscObjectGetComm((PetscObject)mat, &comm);
63: MPI_Comm_size(comm, &size);
64: if (size > 1) {
65: /* create a sequential iscoloring on all processors */
66: MatGetSeqNonzeroStructure(mat, &mat_seq);
67: }
69: MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done);
70: MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done);
73: MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq);
75: PetscMalloc2(n, &list, 4 * n, &work);
77: MINPACKslo(&n, cja, cia, rja, ria, seq, list, &clique, work, work + n, work + 2 * n, work + 3 * n);
79: PetscMalloc1(n, &coloring);
80: MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work);
82: PetscFree2(list, work);
83: PetscFree(seq);
84: MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done);
85: MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done);
87: /* shift coloring numbers to start at zero and shorten */
89: {
90: ISColoringValue *s = (ISColoringValue *)coloring;
91: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
92: MatColoringPatch(mat_seq, ncolors, n, s, iscoloring);
93: }
95: if (size > 1) {
96: MatDestroySeqNonzeroStructure(&mat_seq);
98: /* convert iscoloring_seq to a parallel iscoloring */
99: iscoloring_seq = *iscoloring;
100: rstart = mat->rmap->rstart / bs;
101: rend = mat->rmap->rend / bs;
102: N_loc = rend - rstart; /* number of local nodes */
104: /* get local colors for each local node */
105: PetscMalloc1(N_loc + 1, &colors_loc);
106: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
107: /* create a parallel iscoloring */
108: nc = iscoloring_seq->n;
109: ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring);
110: ISColoringDestroy(&iscoloring_seq);
111: }
112: return 0;
113: }
115: /*MC
116: MATCOLORINGSL - implements the SL (smallest last) coloring routine
118: Level: beginner
120: Notes:
121: Supports only distance two colorings (for computation of Jacobians)
123: This is a sequential algorithm
125: References:
126: . * - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
127: pp. 187-209, 1983.
129: .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
130: M*/
132: PETSC_EXTERN PetscErrorCode MatColoringCreate_SL(MatColoring mc)
133: {
134: mc->dist = 2;
135: mc->data = NULL;
136: mc->ops->apply = MatColoringApply_SL;
137: mc->ops->view = NULL;
138: mc->ops->destroy = NULL;
139: mc->ops->setfromoptions = NULL;
140: return 0;
141: }
143: static PetscErrorCode MatColoringApply_LF(MatColoring mc, ISColoring *iscoloring)
144: {
145: PetscInt *list, *work, *seq, *coloring, n;
146: const PetscInt *ria, *rja, *cia, *cja;
147: PetscInt n1, none, ncolors, i;
148: PetscBool done;
149: Mat mat = mc->mat;
150: Mat mat_seq = mat;
151: PetscMPIInt size;
152: MPI_Comm comm;
153: ISColoring iscoloring_seq;
154: PetscInt bs = 1, rstart, rend, N_loc, nc;
155: ISColoringValue *colors_loc;
156: PetscBool flg1, flg2;
159: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
160: PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1);
161: PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2);
162: if (flg1 || flg2) MatGetBlockSize(mat, &bs);
164: PetscObjectGetComm((PetscObject)mat, &comm);
165: MPI_Comm_size(comm, &size);
166: if (size > 1) {
167: /* create a sequential iscoloring on all processors */
168: MatGetSeqNonzeroStructure(mat, &mat_seq);
169: }
171: MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done);
172: MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done);
175: MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq);
177: PetscMalloc2(n, &list, 4 * n, &work);
179: n1 = n - 1;
180: none = -1;
181: MINPACKnumsrt(&n, &n1, seq, &none, list, work + 2 * n, work + n);
182: PetscMalloc1(n, &coloring);
183: MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work);
185: PetscFree2(list, work);
186: PetscFree(seq);
188: MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done);
189: MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done);
191: /* shift coloring numbers to start at zero and shorten */
193: {
194: ISColoringValue *s = (ISColoringValue *)coloring;
195: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
196: MatColoringPatch(mat_seq, ncolors, n, s, iscoloring);
197: }
199: if (size > 1) {
200: MatDestroySeqNonzeroStructure(&mat_seq);
202: /* convert iscoloring_seq to a parallel iscoloring */
203: iscoloring_seq = *iscoloring;
204: rstart = mat->rmap->rstart / bs;
205: rend = mat->rmap->rend / bs;
206: N_loc = rend - rstart; /* number of local nodes */
208: /* get local colors for each local node */
209: PetscMalloc1(N_loc + 1, &colors_loc);
210: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
212: /* create a parallel iscoloring */
213: nc = iscoloring_seq->n;
214: ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring);
215: ISColoringDestroy(&iscoloring_seq);
216: }
217: return 0;
218: }
220: /*MC
221: MATCOLORINGLF - implements the LF (largest first) coloring routine
223: Level: beginner
225: Notes:
226: Supports only distance two colorings (for computation of Jacobians)
228: This is a sequential algorithm
230: References:
231: . * - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
232: pp. 187-209, 1983.
234: .seealso: `MatColoringTpe`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
235: M*/
237: PETSC_EXTERN PetscErrorCode MatColoringCreate_LF(MatColoring mc)
238: {
239: mc->dist = 2;
240: mc->data = NULL;
241: mc->ops->apply = MatColoringApply_LF;
242: mc->ops->view = NULL;
243: mc->ops->destroy = NULL;
244: mc->ops->setfromoptions = NULL;
245: return 0;
246: }
248: static PetscErrorCode MatColoringApply_ID(MatColoring mc, ISColoring *iscoloring)
249: {
250: PetscInt *list, *work, clique, *seq, *coloring, n;
251: const PetscInt *ria, *rja, *cia, *cja;
252: PetscInt ncolors, i;
253: PetscBool done;
254: Mat mat = mc->mat;
255: Mat mat_seq = mat;
256: PetscMPIInt size;
257: MPI_Comm comm;
258: ISColoring iscoloring_seq;
259: PetscInt bs = 1, rstart, rend, N_loc, nc;
260: ISColoringValue *colors_loc;
261: PetscBool flg1, flg2;
264: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
265: PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQBAIJ, &flg1);
266: PetscObjectBaseTypeCompare((PetscObject)mat, MATMPIBAIJ, &flg2);
267: if (flg1 || flg2) MatGetBlockSize(mat, &bs);
269: PetscObjectGetComm((PetscObject)mat, &comm);
270: MPI_Comm_size(comm, &size);
271: if (size > 1) {
272: /* create a sequential iscoloring on all processors */
273: MatGetSeqNonzeroStructure(mat, &mat_seq);
274: }
276: MatGetRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &ria, &rja, &done);
277: MatGetColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, &n, &cia, &cja, &done);
280: MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq);
282: PetscMalloc2(n, &list, 4 * n, &work);
284: MINPACKido(&n, &n, cja, cia, rja, ria, seq, list, &clique, work, work + n, work + 2 * n, work + 3 * n);
286: PetscMalloc1(n, &coloring);
287: MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work);
289: PetscFree2(list, work);
290: PetscFree(seq);
292: MatRestoreRowIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &ria, &rja, &done);
293: MatRestoreColumnIJ(mat_seq, 1, PETSC_FALSE, PETSC_TRUE, NULL, &cia, &cja, &done);
295: /* shift coloring numbers to start at zero and shorten */
297: {
298: ISColoringValue *s = (ISColoringValue *)coloring;
299: for (i = 0; i < n; i++) s[i] = (ISColoringValue)(coloring[i] - 1);
300: MatColoringPatch(mat_seq, ncolors, n, s, iscoloring);
301: }
303: if (size > 1) {
304: MatDestroySeqNonzeroStructure(&mat_seq);
306: /* convert iscoloring_seq to a parallel iscoloring */
307: iscoloring_seq = *iscoloring;
308: rstart = mat->rmap->rstart / bs;
309: rend = mat->rmap->rend / bs;
310: N_loc = rend - rstart; /* number of local nodes */
312: /* get local colors for each local node */
313: PetscMalloc1(N_loc + 1, &colors_loc);
314: for (i = rstart; i < rend; i++) colors_loc[i - rstart] = iscoloring_seq->colors[i];
315: /* create a parallel iscoloring */
316: nc = iscoloring_seq->n;
317: ISColoringCreate(comm, nc, N_loc, colors_loc, PETSC_OWN_POINTER, iscoloring);
318: ISColoringDestroy(&iscoloring_seq);
319: }
320: return 0;
321: }
323: /*MC
324: MATCOLORINGID - implements the ID (incidence degree) coloring routine
326: Level: beginner
328: Notes:
329: Supports only distance two colorings (for computation of Jacobians)
331: This is a sequential algorithm
333: References:
334: . * - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
335: pp. 187-209, 1983.
337: .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
338: M*/
340: PETSC_EXTERN PetscErrorCode MatColoringCreate_ID(MatColoring mc)
341: {
342: mc->dist = 2;
343: mc->data = NULL;
344: mc->ops->apply = MatColoringApply_ID;
345: mc->ops->view = NULL;
346: mc->ops->destroy = NULL;
347: mc->ops->setfromoptions = NULL;
348: return 0;
349: }