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