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