1 /*
2 Routines that call the kernel minpack coloring subroutines
3 */
4
5 #include <petsc/private/matimpl.h>
6 #include <petsc/private/isimpl.h>
7 #include <../src/mat/graphops/color/impls/minpack/color.h>
8
9 /*
10 MatFDColoringDegreeSequence_Minpack - Calls the MINPACK routine seqr() that
11 computes the degree sequence required by MINPACK coloring routines.
12 */
MatFDColoringDegreeSequence_Minpack(PetscInt m,const PetscInt * cja,const PetscInt * cia,const PetscInt * rja,const PetscInt * ria,PetscInt ** seq)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;
16
17 PetscFunctionBegin;
18 PetscCall(PetscMalloc1(m, &work));
19 PetscCall(PetscMalloc1(m, seq));
20
21 PetscCall(MINPACKdegr(&m, cja, cia, rja, ria, *seq, work));
22
23 PetscCall(PetscFree(work));
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
MatColoringApply_SL(MatColoring mc,ISColoring * iscoloring)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;
41
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));
48
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 }
55
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");
59
60 PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
61
62 PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
63
64 PetscCall(MINPACKslo(&n, cja, cia, rja, ria, seq, list, &clique, work, work + n, work + 2 * n, work + 3 * n));
65
66 PetscCall(PetscMalloc1(n, &coloring));
67 PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
68
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));
73
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 }
81
82 if (size > 1) {
83 PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
84
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 */
90
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 }
101
102 /*MC
103 MATCOLORINGSL - implements the SL (smallest last) coloring routine {cite}`more:coloring`
104
105 Level: beginner
106
107 Notes:
108 Supports only distance two colorings (for computation of Jacobians/Hessians)
109
110 This is a sequential algorithm
111
112 .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
113 M*/
114
MatColoringCreate_SL(MatColoring mc)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 }
126
MatColoringApply_LF(MatColoring mc,ISColoring * iscoloring)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;
141
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));
148
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 }
155
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");
159
160 PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
161
162 PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
163
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));
169
170 PetscCall(PetscFree2(list, work));
171 PetscCall(PetscFree(seq));
172
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));
175
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 }
183
184 if (size > 1) {
185 PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
186
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 */
192
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];
196
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 }
204
205 /*MC
206 MATCOLORINGLF - implements the LF (largest first) coloring routine {cite}`more:coloring`
207
208 Level: beginner
209
210 Notes:
211 Supports only distance two colorings (for computation of Jacobians/Hessians)
212
213 This is a sequential algorithm
214
215 .seealso: `MatColoringTpe`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
216 M*/
217
MatColoringCreate_LF(MatColoring mc)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 }
229
MatColoringApply_ID(MatColoring mc,ISColoring * iscoloring)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;
244
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));
251
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 }
258
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");
262
263 PetscCall(MatFDColoringDegreeSequence_Minpack(n, cja, cia, rja, ria, &seq));
264
265 PetscCall(PetscMalloc2(n, &list, 4 * n, &work));
266
267 PetscCall(MINPACKido(&n, &n, cja, cia, rja, ria, seq, list, &clique, work, work + n, work + 2 * n, work + 3 * n));
268
269 PetscCall(PetscMalloc1(n, &coloring));
270 PetscCall(MINPACKseq(&n, cja, cia, rja, ria, list, coloring, &ncolors, work));
271
272 PetscCall(PetscFree2(list, work));
273 PetscCall(PetscFree(seq));
274
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));
277
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 }
285
286 if (size > 1) {
287 PetscCall(MatDestroySeqNonzeroStructure(&mat_seq));
288
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 */
294
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 }
305
306 /*MC
307 MATCOLORINGID - implements the ID (incidence degree) coloring routine {cite}`more:coloring`
308
309 Level: beginner
310
311 Notes:
312 Supports only distance two colorings (for computation of Jacobians/Hessians)
313
314 This is a sequential algorithm
315
316 .seealso: `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MATCOLORINGGREEDY`, `MatColoringType`
317 M*/
318
MatColoringCreate_ID(MatColoring mc)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 }
330