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: }