Actual source code: color.c

petsc-master 2019-04-24
Report Typos and Errors

  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;

 20:   PetscMalloc1(m,&work);
 21:   PetscMalloc1(m,seq);

 23:   MINPACKdegr(&m,cja,cia,rja,ria,*seq,work);

 25:   PetscFree(work);
 26:   return(0);
 27: }

 29: /*
 30:     MatFDColoringMinimumNumberofColors_Private - For a given sparse
 31:         matrix computes the minimum number of colors needed.

 33: */
 34: PetscErrorCode MatFDColoringMinimumNumberofColors_Private(PetscInt m,PetscInt *ia,PetscInt *minc)
 35: {
 36:   PetscInt i,c = 0;

 39:   for (i=0; i<m; i++) c = PetscMax(c,ia[i+1]-ia[i]);
 40:   *minc = c;
 41:   return(0);
 42: }

 44: static PetscErrorCode MatColoringApply_SL(MatColoring mc,ISColoring *iscoloring)
 45: {
 46:   PetscErrorCode  ierr;
 47:   PetscInt        *list,*work,clique,*seq,*coloring,n;
 48:   const PetscInt  *ria,*rja,*cia,*cja;
 49:   PetscInt        ncolors,i;
 50:   PetscBool       done;
 51:   Mat             mat = mc->mat;
 52:   Mat             mat_seq = mat;
 53:   PetscMPIInt     size;
 54:   MPI_Comm        comm;
 55:   ISColoring      iscoloring_seq;
 56:   PetscInt        bs = 1,rstart,rend,N_loc,nc;
 57:   ISColoringValue *colors_loc;
 58:   PetscBool       flg1,flg2;

 61:   if (mc->dist != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SL may only do distance 2 coloring");
 62:   /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
 63:   PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
 64:   PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
 65:   if (flg1 || flg2) {
 66:     MatGetBlockSize(mat,&bs);
 67:   }

 69:   PetscObjectGetComm((PetscObject)mat,&comm);
 70:   MPI_Comm_size(comm,&size);
 71:   if (size > 1) {
 72:     /* create a sequential iscoloring on all processors */
 73:     MatGetSeqNonzeroStructure(mat,&mat_seq);
 74:   }

 76:   MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
 77:   MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
 78:   if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");

 80:   MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);

 82:   PetscMalloc2(n,&list,4*n,&work);

 84:   MINPACKslo(&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);

 86:   PetscMalloc1(n,&coloring);
 87:   MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);

 89:   PetscFree2(list,work);
 90:   PetscFree(seq);
 91:   MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
 92:   MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);

 94:   /* shift coloring numbers to start at zero and shorten */
 95:   if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
 96:   {
 97:     ISColoringValue *s = (ISColoringValue*) coloring;
 98:     for (i=0; i<n; i++) {
 99:       s[i] = (ISColoringValue) (coloring[i]-1);
100:     }
101:     MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
102:   }

104:   if (size > 1) {
105:     MatDestroySeqNonzeroStructure(&mat_seq);

107:     /* convert iscoloring_seq to a parallel iscoloring */
108:     iscoloring_seq = *iscoloring;
109:     rstart         = mat->rmap->rstart/bs;
110:     rend           = mat->rmap->rend/bs;
111:     N_loc          = rend - rstart; /* number of local nodes */

113:     /* get local colors for each local node */
114:     PetscMalloc1(N_loc+1,&colors_loc);
115:     for (i=rstart; i<rend; i++) {
116:       colors_loc[i-rstart] = iscoloring_seq->colors[i];
117:     }
118:     /* create a parallel iscoloring */
119:     nc   = iscoloring_seq->n;
120:     ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
121:     ISColoringDestroy(&iscoloring_seq);
122:   }
123:   return(0);
124: }

126: /*MC
127:   MATCOLORINGSL - implements the SL (smallest last) coloring routine

129:    Level: beginner

131:    Notes:
132:     Supports only distance two colorings (for computation of Jacobians)

134:           This is a sequential algorithm

136:    References:
137: .  1. - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
138:    pp. 187-209, 1983.

140: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
141: M*/

143: PETSC_EXTERN PetscErrorCode MatColoringCreate_SL(MatColoring mc)
144: {
146:     mc->dist                = 2;
147:     mc->data                = NULL;
148:     mc->ops->apply          = MatColoringApply_SL;
149:     mc->ops->view           = NULL;
150:     mc->ops->destroy        = NULL;
151:     mc->ops->setfromoptions = NULL;
152:     return(0);
153: }

155: static PetscErrorCode MatColoringApply_LF(MatColoring mc,ISColoring *iscoloring)
156: {
157:   PetscErrorCode  ierr;
158:   PetscInt        *list,*work,*seq,*coloring,n;
159:   const PetscInt  *ria,*rja,*cia,*cja;
160:   PetscInt        n1, none,ncolors,i;
161:   PetscBool       done;
162:   Mat             mat = mc->mat;
163:   Mat             mat_seq = mat;
164:   PetscMPIInt     size;
165:   MPI_Comm        comm;
166:   ISColoring      iscoloring_seq;
167:   PetscInt        bs = 1,rstart,rend,N_loc,nc;
168:   ISColoringValue *colors_loc;
169:   PetscBool       flg1,flg2;

172:   if (mc->dist != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"LF may only do distance 2 coloring");
173:   /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
174:   PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
175:   PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
176:   if (flg1 || flg2) {
177:     MatGetBlockSize(mat,&bs);
178:   }

180:   PetscObjectGetComm((PetscObject)mat,&comm);
181:   MPI_Comm_size(comm,&size);
182:   if (size > 1) {
183:     /* create a sequential iscoloring on all processors */
184:     MatGetSeqNonzeroStructure(mat,&mat_seq);
185:   }

187:   MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
188:   MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
189:   if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");

191:   MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);

193:   PetscMalloc2(n,&list,4*n,&work);

195:   n1   = n - 1;
196:   none = -1;
197:   MINPACKnumsrt(&n,&n1,seq,&none,list,work+2*n,work+n);
198:   PetscMalloc1(n,&coloring);
199:   MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);

201:   PetscFree2(list,work);
202:   PetscFree(seq);

204:   MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
205:   MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);

207:   /* shift coloring numbers to start at zero and shorten */
208:   if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
209:   {
210:     ISColoringValue *s = (ISColoringValue*) coloring;
211:     for (i=0; i<n; i++) s[i] = (ISColoringValue) (coloring[i]-1);
212:     MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
213:   }

215:   if (size > 1) {
216:     MatDestroySeqNonzeroStructure(&mat_seq);

218:     /* convert iscoloring_seq to a parallel iscoloring */
219:     iscoloring_seq = *iscoloring;
220:     rstart         = mat->rmap->rstart/bs;
221:     rend           = mat->rmap->rend/bs;
222:     N_loc          = rend - rstart; /* number of local nodes */

224:     /* get local colors for each local node */
225:     PetscMalloc1(N_loc+1,&colors_loc);
226:     for (i=rstart; i<rend; i++) colors_loc[i-rstart] = iscoloring_seq->colors[i];

228:     /* create a parallel iscoloring */
229:     nc   = iscoloring_seq->n;
230:     ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
231:     ISColoringDestroy(&iscoloring_seq);
232:   }
233:   return(0);
234: }

236: /*MC
237:   MATCOLORINGLF - implements the LF (largest first) coloring routine

239:    Level: beginner

241:    Notes:
242:     Supports only distance two colorings (for computation of Jacobians)

244:           This is a sequential algorithm

246:    References:
247: .  1. - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
248:    pp. 187-209, 1983.

250: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
251: M*/

253: PETSC_EXTERN PetscErrorCode MatColoringCreate_LF(MatColoring mc)
254: {
256:     mc->dist                = 2;
257:     mc->data                = NULL;
258:     mc->ops->apply          = MatColoringApply_LF;
259:     mc->ops->view           = NULL;
260:     mc->ops->destroy        = NULL;
261:     mc->ops->setfromoptions = NULL;
262:     return(0);
263: }

265: static PetscErrorCode MatColoringApply_ID(MatColoring mc,ISColoring *iscoloring)
266: {
267:   PetscErrorCode  ierr;
268:   PetscInt        *list,*work,clique,*seq,*coloring,n;
269:   const PetscInt  *ria,*rja,*cia,*cja;
270:   PetscInt        ncolors,i;
271:   PetscBool       done;
272:   Mat             mat = mc->mat;
273:   Mat             mat_seq = mat;
274:   PetscMPIInt     size;
275:   MPI_Comm        comm;
276:   ISColoring      iscoloring_seq;
277:   PetscInt        bs = 1,rstart,rend,N_loc,nc;
278:   ISColoringValue *colors_loc;
279:   PetscBool       flg1,flg2;

282:   if (mc->dist != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IDO may only do distance 2 coloring");
283:   /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
284:   PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
285:   PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
286:   if (flg1 || flg2) {
287:     MatGetBlockSize(mat,&bs);
288:   }

290:   PetscObjectGetComm((PetscObject)mat,&comm);
291:   MPI_Comm_size(comm,&size);
292:   if (size > 1) {
293:     /* create a sequential iscoloring on all processors */
294:     MatGetSeqNonzeroStructure(mat,&mat_seq);
295:   }

297:   MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
298:   MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
299:   if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");

301:   MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);

303:   PetscMalloc2(n,&list,4*n,&work);

305:   MINPACKido(&n,&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);

307:   PetscMalloc1(n,&coloring);
308:   MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);

310:   PetscFree2(list,work);
311:   PetscFree(seq);

313:   MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
314:   MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);

316:   /* shift coloring numbers to start at zero and shorten */
317:   if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
318:   {
319:     ISColoringValue *s = (ISColoringValue*) coloring;
320:     for (i=0; i<n; i++) {
321:       s[i] = (ISColoringValue) (coloring[i]-1);
322:     }
323:     MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
324:   }

326:   if (size > 1) {
327:     MatDestroySeqNonzeroStructure(&mat_seq);

329:     /* convert iscoloring_seq to a parallel iscoloring */
330:     iscoloring_seq = *iscoloring;
331:     rstart         = mat->rmap->rstart/bs;
332:     rend           = mat->rmap->rend/bs;
333:     N_loc          = rend - rstart; /* number of local nodes */

335:     /* get local colors for each local node */
336:     PetscMalloc1(N_loc+1,&colors_loc);
337:     for (i=rstart; i<rend; i++) {
338:       colors_loc[i-rstart] = iscoloring_seq->colors[i];
339:     }
340:     /* create a parallel iscoloring */
341:     nc   = iscoloring_seq->n;
342:     ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
343:     ISColoringDestroy(&iscoloring_seq);
344:   }
345:   return(0);
346: }

348: /*MC
349:   MATCOLORINGID - implements the ID (incidence degree) coloring routine

351:    Level: beginner

353:    Notes:
354:     Supports only distance two colorings (for computation of Jacobians)

356:           This is a sequential algorithm

358:    References:
359: .  1. - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
360:    pp. 187-209, 1983.

362: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
363: M*/

365: PETSC_EXTERN PetscErrorCode MatColoringCreate_ID(MatColoring mc)
366: {
368:     mc->dist                = 2;
369:     mc->data                = NULL;
370:     mc->ops->apply          = MatColoringApply_ID;
371:     mc->ops->view           = NULL;
372:     mc->ops->destroy        = NULL;
373:     mc->ops->setfromoptions = NULL;
374:     return(0);
375: }