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) {
 61:     MatGetBlockSize(mat,&bs);
 62:   }

 64:   PetscObjectGetComm((PetscObject)mat,&comm);
 65:   MPI_Comm_size(comm,&size);
 66:   if (size > 1) {
 67:     /* create a sequential iscoloring on all processors */
 68:     MatGetSeqNonzeroStructure(mat,&mat_seq);
 69:   }

 71:   MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
 72:   MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);

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

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

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

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

 84:   PetscFree2(list,work);
 85:   PetscFree(seq);
 86:   MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
 87:   MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);

 89:   /* shift coloring numbers to start at zero and shorten */
 91:   {
 92:     ISColoringValue *s = (ISColoringValue*) coloring;
 93:     for (i=0; i<n; i++) {
 94:       s[i] = (ISColoringValue) (coloring[i]-1);
 95:     }
 96:     MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
 97:   }

 99:   if (size > 1) {
100:     MatDestroySeqNonzeroStructure(&mat_seq);

102:     /* convert iscoloring_seq to a parallel iscoloring */
103:     iscoloring_seq = *iscoloring;
104:     rstart         = mat->rmap->rstart/bs;
105:     rend           = mat->rmap->rend/bs;
106:     N_loc          = rend - rstart; /* number of local nodes */

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

121: /*MC
122:   MATCOLORINGSL - implements the SL (smallest last) coloring routine

124:    Level: beginner

126:    Notes:
127:     Supports only distance two colorings (for computation of Jacobians)

129:           This is a sequential algorithm

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

135: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
136: M*/

138: PETSC_EXTERN PetscErrorCode MatColoringCreate_SL(MatColoring mc)
139: {
140:     mc->dist                = 2;
141:     mc->data                = NULL;
142:     mc->ops->apply          = MatColoringApply_SL;
143:     mc->ops->view           = NULL;
144:     mc->ops->destroy        = NULL;
145:     mc->ops->setfromoptions = NULL;
146:     return 0;
147: }

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

165:   /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
166:   PetscObjectBaseTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
167:   PetscObjectBaseTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
168:   if (flg1 || flg2) {
169:     MatGetBlockSize(mat,&bs);
170:   }

172:   PetscObjectGetComm((PetscObject)mat,&comm);
173:   MPI_Comm_size(comm,&size);
174:   if (size > 1) {
175:     /* create a sequential iscoloring on all processors */
176:     MatGetSeqNonzeroStructure(mat,&mat_seq);
177:   }

179:   MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
180:   MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);

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

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

187:   n1   = n - 1;
188:   none = -1;
189:   MINPACKnumsrt(&n,&n1,seq,&none,list,work+2*n,work+n);
190:   PetscMalloc1(n,&coloring);
191:   MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);

193:   PetscFree2(list,work);
194:   PetscFree(seq);

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

199:   /* shift coloring numbers to start at zero and shorten */
201:   {
202:     ISColoringValue *s = (ISColoringValue*) coloring;
203:     for (i=0; i<n; i++) s[i] = (ISColoringValue) (coloring[i]-1);
204:     MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
205:   }

207:   if (size > 1) {
208:     MatDestroySeqNonzeroStructure(&mat_seq);

210:     /* convert iscoloring_seq to a parallel iscoloring */
211:     iscoloring_seq = *iscoloring;
212:     rstart         = mat->rmap->rstart/bs;
213:     rend           = mat->rmap->rend/bs;
214:     N_loc          = rend - rstart; /* number of local nodes */

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

220:     /* create a parallel iscoloring */
221:     nc   = iscoloring_seq->n;
222:     ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
223:     ISColoringDestroy(&iscoloring_seq);
224:   }
225:   return 0;
226: }

228: /*MC
229:   MATCOLORINGLF - implements the LF (largest first) coloring routine

231:    Level: beginner

233:    Notes:
234:     Supports only distance two colorings (for computation of Jacobians)

236:           This is a sequential algorithm

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

242: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
243: M*/

245: PETSC_EXTERN PetscErrorCode MatColoringCreate_LF(MatColoring mc)
246: {
247:     mc->dist                = 2;
248:     mc->data                = NULL;
249:     mc->ops->apply          = MatColoringApply_LF;
250:     mc->ops->view           = NULL;
251:     mc->ops->destroy        = NULL;
252:     mc->ops->setfromoptions = NULL;
253:     return 0;
254: }

256: static PetscErrorCode MatColoringApply_ID(MatColoring mc,ISColoring *iscoloring)
257: {
258:   PetscInt        *list,*work,clique,*seq,*coloring,n;
259:   const PetscInt  *ria,*rja,*cia,*cja;
260:   PetscInt        ncolors,i;
261:   PetscBool       done;
262:   Mat             mat = mc->mat;
263:   Mat             mat_seq = mat;
264:   PetscMPIInt     size;
265:   MPI_Comm        comm;
266:   ISColoring      iscoloring_seq;
267:   PetscInt        bs = 1,rstart,rend,N_loc,nc;
268:   ISColoringValue *colors_loc;
269:   PetscBool       flg1,flg2;

272:   /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
273:   PetscObjectBaseTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
274:   PetscObjectBaseTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
275:   if (flg1 || flg2) {
276:     MatGetBlockSize(mat,&bs);
277:   }

279:   PetscObjectGetComm((PetscObject)mat,&comm);
280:   MPI_Comm_size(comm,&size);
281:   if (size > 1) {
282:     /* create a sequential iscoloring on all processors */
283:     MatGetSeqNonzeroStructure(mat,&mat_seq);
284:   }

286:   MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
287:   MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);

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

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

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

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

299:   PetscFree2(list,work);
300:   PetscFree(seq);

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

305:   /* shift coloring numbers to start at zero and shorten */
307:   {
308:     ISColoringValue *s = (ISColoringValue*) coloring;
309:     for (i=0; i<n; i++) {
310:       s[i] = (ISColoringValue) (coloring[i]-1);
311:     }
312:     MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
313:   }

315:   if (size > 1) {
316:     MatDestroySeqNonzeroStructure(&mat_seq);

318:     /* convert iscoloring_seq to a parallel iscoloring */
319:     iscoloring_seq = *iscoloring;
320:     rstart         = mat->rmap->rstart/bs;
321:     rend           = mat->rmap->rend/bs;
322:     N_loc          = rend - rstart; /* number of local nodes */

324:     /* get local colors for each local node */
325:     PetscMalloc1(N_loc+1,&colors_loc);
326:     for (i=rstart; i<rend; i++) {
327:       colors_loc[i-rstart] = iscoloring_seq->colors[i];
328:     }
329:     /* create a parallel iscoloring */
330:     nc   = iscoloring_seq->n;
331:     ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
332:     ISColoringDestroy(&iscoloring_seq);
333:   }
334:   return 0;
335: }

337: /*MC
338:   MATCOLORINGID - implements the ID (incidence degree) coloring routine

340:    Level: beginner

342:    Notes:
343:     Supports only distance two colorings (for computation of Jacobians)

345:           This is a sequential algorithm

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

351: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
352: M*/

354: PETSC_EXTERN PetscErrorCode MatColoringCreate_ID(MatColoring mc)
355: {
356:     mc->dist                = 2;
357:     mc->data                = NULL;
358:     mc->ops->apply          = MatColoringApply_ID;
359:     mc->ops->view           = NULL;
360:     mc->ops->destroy        = NULL;
361:     mc->ops->setfromoptions = NULL;
362:     return 0;
363: }