Actual source code: color.c
petsc-3.3-p7 2013-05-11
2: /*
3: Routines that call the kernel minpack coloring subroutines
4: */
6: #include <petsc-private/matimpl.h>
7: #include <../src/mat/color/color.h>
9: /*
10: MatFDColoringDegreeSequence_Minpack - Calls the MINPACK routine seqr() that
11: computes the degree sequence required by MINPACK coloring routines.
12: */
15: PetscErrorCode MatFDColoringDegreeSequence_Minpack(PetscInt m,PetscInt *cja, PetscInt *cia, PetscInt *rja, PetscInt *ria, PetscInt **seq)
16: {
17: PetscInt *work;
21: PetscMalloc(m*sizeof(PetscInt),&work);
22: PetscMalloc(m*sizeof(PetscInt),seq);
24: MINPACKdegr(&m,cja,cia,rja,ria,*seq,work);
26: PetscFree(work);
27: return(0);
28: }
30: /*
31: MatFDColoringMinimumNumberofColors_Private - For a given sparse
32: matrix computes the minimum number of colors needed.
34: */
37: PetscErrorCode MatFDColoringMinimumNumberofColors_Private(PetscInt m,PetscInt *ia,PetscInt *minc)
38: {
39: PetscInt i,c = 0;
42: for (i=0; i<m; i++) {
43: c = PetscMax(c,ia[i+1]-ia[i]);
44: }
45: *minc = c;
46: return(0);
47: }
49: EXTERN_C_BEGIN
50: /* ----------------------------------------------------------------------------*/
51: /*
52: MatGetColoring_SL_Minpack - Uses the smallest-last (SL) coloring of minpack
53: */
56: PetscErrorCode MatGetColoring_SL_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
57: {
58: PetscErrorCode ierr;
59: PetscInt *list,*work,clique,*ria,*rja,*cia,*cja,*seq,*coloring,n;
60: PetscInt ncolors,i;
61: PetscBool done;
62: Mat mat_seq = mat;
63: PetscMPIInt size;
64: MPI_Comm comm;
65: ISColoring iscoloring_seq;
66: PetscInt bs = 1,rstart,rend,N_loc,nc;
67: ISColoringValue *colors_loc;
68: PetscBool flg1,flg2;
71: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
72: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
73: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
74: if (flg1 || flg2) {
75: MatGetBlockSize(mat,&bs);
76: }
78: PetscObjectGetComm((PetscObject)mat,&comm);
79: MPI_Comm_size(comm,&size);
80: if (size > 1){
81: /* create a sequential iscoloring on all processors */
82: MatGetSeqNonzeroStructure(mat,&mat_seq);
83: }
85: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
86: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
87: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
89: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
91: PetscMalloc2(n,PetscInt,&list,4*n,PetscInt,&work);
93: MINPACKslo(&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
95: PetscMalloc(n*sizeof(PetscInt),&coloring);
96: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
98: PetscFree2(list,work);
99: PetscFree(seq);
100: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
101: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
103: /* shift coloring numbers to start at zero and shorten */
104: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
105: {
106: ISColoringValue *s = (ISColoringValue*) coloring;
107: for (i=0; i<n; i++) {
108: s[i] = (ISColoringValue) (coloring[i]-1);
109: }
110: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
111: }
113: if (size > 1) {
114: MatDestroySeqNonzeroStructure(&mat_seq);
116: /* convert iscoloring_seq to a parallel iscoloring */
117: iscoloring_seq = *iscoloring;
118: rstart = mat->rmap->rstart/bs;
119: rend = mat->rmap->rend/bs;
120: N_loc = rend - rstart; /* number of local nodes */
122: /* get local colors for each local node */
123: PetscMalloc((N_loc+1)*sizeof(ISColoringValue),&colors_loc);
124: for (i=rstart; i<rend; i++){
125: colors_loc[i-rstart] = iscoloring_seq->colors[i];
126: }
127: /* create a parallel iscoloring */
128: nc=iscoloring_seq->n;
129: ISColoringCreate(comm,nc,N_loc,colors_loc,iscoloring);
130: ISColoringDestroy(&iscoloring_seq);
131: }
132: return(0);
133: }
134: EXTERN_C_END
136: EXTERN_C_BEGIN
137: /* ----------------------------------------------------------------------------*/
138: /*
139: MatGetColoring_LF_Minpack -
140: */
143: PetscErrorCode MatGetColoring_LF_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
144: {
145: PetscErrorCode ierr;
146: PetscInt *list,*work,*ria,*rja,*cia,*cja,*seq,*coloring,n;
147: PetscInt n1, none,ncolors,i;
148: PetscBool done;
149: Mat mat_seq = mat;
150: PetscMPIInt size;
151: MPI_Comm comm;
152: ISColoring iscoloring_seq;
153: PetscInt bs = 1,rstart,rend,N_loc,nc;
154: ISColoringValue *colors_loc;
155: PetscBool flg1,flg2;
158: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
159: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
160: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
161: if (flg1 || flg2) {
162: MatGetBlockSize(mat,&bs);
163: }
165: PetscObjectGetComm((PetscObject)mat,&comm);
166: MPI_Comm_size(comm,&size);
167: if (size > 1){
168: /* create a sequential iscoloring on all processors */
169: MatGetSeqNonzeroStructure(mat,&mat_seq);
170: }
172: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
173: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
174: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
176: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
178: PetscMalloc2(n,PetscInt,&list,4*n,PetscInt,&work);
180: n1 = n - 1;
181: none = -1;
182: MINPACKnumsrt(&n,&n1,seq,&none,list,work+2*n,work+n);
183: PetscMalloc(n*sizeof(PetscInt),&coloring);
184: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
186: PetscFree2(list,work);
187: PetscFree(seq);
189: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
190: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
192: /* shift coloring numbers to start at zero and shorten */
193: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
194: {
195: ISColoringValue *s = (ISColoringValue*) coloring;
196: for (i=0; i<n; i++) {
197: s[i] = (ISColoringValue) (coloring[i]-1);
198: }
199: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
200: }
202: if (size > 1) {
203: MatDestroySeqNonzeroStructure(&mat_seq);
205: /* convert iscoloring_seq to a parallel iscoloring */
206: iscoloring_seq = *iscoloring;
207: rstart = mat->rmap->rstart/bs;
208: rend = mat->rmap->rend/bs;
209: N_loc = rend - rstart; /* number of local nodes */
211: /* get local colors for each local node */
212: PetscMalloc((N_loc+1)*sizeof(ISColoringValue),&colors_loc);
213: for (i=rstart; i<rend; i++){
214: colors_loc[i-rstart] = iscoloring_seq->colors[i];
215: }
216: /* create a parallel iscoloring */
217: nc=iscoloring_seq->n;
218: ISColoringCreate(comm,nc,N_loc,colors_loc,iscoloring);
219: ISColoringDestroy(&iscoloring_seq);
220: }
221: return(0);
222: }
223: EXTERN_C_END
225: EXTERN_C_BEGIN
226: /* ----------------------------------------------------------------------------*/
227: /*
228: MatGetColoring_ID_Minpack -
229: */
232: PetscErrorCode MatGetColoring_ID_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
233: {
234: PetscErrorCode ierr;
235: PetscInt *list,*work,clique,*ria,*rja,*cia,*cja,*seq,*coloring,n;
236: PetscInt ncolors,i;
237: PetscBool done;
238: Mat mat_seq = mat;
239: PetscMPIInt size;
240: MPI_Comm comm;
241: ISColoring iscoloring_seq;
242: PetscInt bs = 1,rstart,rend,N_loc,nc;
243: ISColoringValue *colors_loc;
244: PetscBool flg1,flg2;
247: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
248: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
249: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
250: if (flg1 || flg2) {
251: MatGetBlockSize(mat,&bs);
252: }
254: PetscObjectGetComm((PetscObject)mat,&comm);
255: MPI_Comm_size(comm,&size);
256: if (size > 1){
257: /* create a sequential iscoloring on all processors */
258: MatGetSeqNonzeroStructure(mat,&mat_seq);
259: }
261: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
262: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
263: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
265: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
267: PetscMalloc2(n,PetscInt,&list,4*n,PetscInt,&work);
269: MINPACKido(&n,&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
271: PetscMalloc(n*sizeof(PetscInt),&coloring);
272: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
274: PetscFree2(list,work);
275: PetscFree(seq);
277: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
278: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
280: /* shift coloring numbers to start at zero and shorten */
281: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
282: {
283: ISColoringValue *s = (ISColoringValue*) coloring;
284: for (i=0; i<n; i++) {
285: s[i] = (ISColoringValue) (coloring[i]-1);
286: }
287: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
288: }
290: if (size > 1) {
291: MatDestroySeqNonzeroStructure(&mat_seq);
293: /* convert iscoloring_seq to a parallel iscoloring */
294: iscoloring_seq = *iscoloring;
295: rstart = mat->rmap->rstart/bs;
296: rend = mat->rmap->rend/bs;
297: N_loc = rend - rstart; /* number of local nodes */
299: /* get local colors for each local node */
300: PetscMalloc((N_loc+1)*sizeof(ISColoringValue),&colors_loc);
301: for (i=rstart; i<rend; i++){
302: colors_loc[i-rstart] = iscoloring_seq->colors[i];
303: }
304: /* create a parallel iscoloring */
305: nc=iscoloring_seq->n;
306: ISColoringCreate(comm,nc,N_loc,colors_loc,iscoloring);
307: ISColoringDestroy(&iscoloring_seq);
308: }
309: return(0);
310: }
311: EXTERN_C_END
313: EXTERN_C_BEGIN
314: /*
315: Simplest coloring, each column of the matrix gets its own unique color.
316: */
319: PetscErrorCode MatGetColoring_Natural(Mat mat,MatColoringType color, ISColoring *iscoloring)
320: {
321: PetscErrorCode ierr;
322: PetscInt start,end,i,bs = 1,n;
323: ISColoringValue *colors;
324: MPI_Comm comm;
325: PetscBool flg1,flg2;
326: Mat mat_seq = mat;
327: PetscMPIInt size;
328: ISColoring iscoloring_seq;
329: ISColoringValue *colors_loc;
330: PetscInt rstart,rend,N_loc,nc;
333: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
334: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
335: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
336: if (flg1 || flg2) {
337: MatGetBlockSize(mat,&bs);
338: }
340: PetscObjectGetComm((PetscObject)mat,&comm);
341: MPI_Comm_size(comm,&size);
342: if (size > 1){
343: /* create a sequential iscoloring on all processors */
344: MatGetSeqNonzeroStructure(mat,&mat_seq);
345: }
347: MatGetSize(mat_seq,PETSC_NULL,&n);
348: MatGetOwnershipRange(mat_seq,&start,&end);
349: n = n/bs;
350: if (n > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
352: start = start/bs;
353: end = end/bs;
354: PetscMalloc((end-start+1)*sizeof(PetscInt),&colors);
355: for (i=start; i<end; i++) {
356: colors[i-start] = (ISColoringValue)i;
357: }
358: ISColoringCreate(comm,n,end-start,colors,iscoloring);
360: if (size > 1) {
361: MatDestroySeqNonzeroStructure(&mat_seq);
363: /* convert iscoloring_seq to a parallel iscoloring */
364: iscoloring_seq = *iscoloring;
365: rstart = mat->rmap->rstart/bs;
366: rend = mat->rmap->rend/bs;
367: N_loc = rend - rstart; /* number of local nodes */
369: /* get local colors for each local node */
370: PetscMalloc((N_loc+1)*sizeof(ISColoringValue),&colors_loc);
371: for (i=rstart; i<rend; i++){
372: colors_loc[i-rstart] = iscoloring_seq->colors[i];
373: }
374: /* create a parallel iscoloring */
375: nc=iscoloring_seq->n;
376: ISColoringCreate(comm,nc,N_loc,colors_loc,iscoloring);
377: ISColoringDestroy(&iscoloring_seq);
378: }
379: return(0);
380: }
381: EXTERN_C_END
382:
383: /* ===========================================================================================*/
385: PetscFList MatColoringList = 0;
386: PetscBool MatColoringRegisterAllCalled = PETSC_FALSE;
390: PetscErrorCode MatColoringRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(Mat,MatColoringType,ISColoring*))
391: {
393: char fullname[PETSC_MAX_PATH_LEN];
396: PetscFListConcat(path,name,fullname);
397: PetscFListAdd(&MatColoringList,sname,fullname,(void (*)(void))function);
398: return(0);
399: }
403: /*@C
404: MatColoringRegisterDestroy - Frees the list of coloringing routines.
406: Not Collective
408: Level: developer
410: .keywords: matrix, register, destroy
412: .seealso: MatColoringRegisterDynamic(), MatColoringRegisterAll()
413: @*/
414: PetscErrorCode MatColoringRegisterDestroy(void)
415: {
419: PetscFListDestroy(&MatColoringList);
420: MatColoringRegisterAllCalled = PETSC_FALSE;
421: return(0);
422: }
426: /*@C
427: MatGetColoring - Gets a coloring for a matrix, from its sparsity structure,
428: to reduce the number of function evaluations needed to compute a sparse Jacobian via differencing.
430: Collective on Mat
432: Input Parameters:
433: . mat - the matrix
434: . type - type of coloring, one of the following:
435: $ MATCOLORINGNATURAL - natural (one color for each column, very slow)
436: $ MATCOLORINGSL - smallest-last
437: $ MATCOLORINGLF - largest-first
438: $ MATCOLORINGID - incidence-degree
440: Output Parameters:
441: . iscoloring - the coloring
443: Options Database Keys:
444: To specify the coloring through the options database, use one of
445: the following
446: $ -mat_coloring_type natural, -mat_coloring_type sl, -mat_coloring_type lf,
447: $ -mat_coloring_type id
448: To see the coloring use
449: $ -mat_coloring_view
451: Level: intermediate
453: Notes:
454: $ A graph coloring C(A) is a division of vertices so that two vertices of the same color do not share any common edges.
455: $ A suitable coloring for a smoother is simply C(A).
456: $ A suitable coloring for efficient Jacobian computation is a division of the columns so that two columns of the same color do not share any common rows.
457: $ This corresponds to C(A^{T} A). This is what MatGetColoring() computes.
459: The user can define additional colorings; see MatColoringRegisterDynamic().
461: For parallel matrices currently converts to sequential matrix and uses the sequential coloring
462: on that.
464: The colorings SL, LF, and ID are obtained via the Minpack software that was
465: converted to C using f2c.
467: For BAIJ matrices this colors the blocks. The true number of colors would be block size times the number of colors
468: returned here.
470: References:
471: $ Thomas F. Coleman and Jorge J. More, Estimation of Sparse {J}acobian Matrices and Graph Coloring Problems,
472: $ SIAM Journal on Numerical Analysis, 1983, pages 187-209, volume 20
473: $ Jorge J. Mor\'{e} and Danny C. Sorenson and Burton S. Garbow and Kenneth E. Hillstrom, The {MINPACK} Project,
474: $ Sources and Development of Mathematical Software, Wayne R. Cowell editor, 1984, pages 88-111
476: .keywords: matrix, get, coloring
478: .seealso: MatGetColoringTypeFromOptions(), MatColoringRegisterDynamic(), MatFDColoringCreate(),
479: SNESDefaultComputeJacobianColor()
480: @*/
481: PetscErrorCode MatGetColoring(Mat mat,const MatColoringType type,ISColoring *iscoloring)
482: {
483: PetscBool flag;
484: PetscErrorCode ierr,(*r)(Mat,const MatColoringType,ISColoring *);
485: char tname[PETSC_MAX_PATH_LEN];
486: MPI_Comm comm;
491: if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
492: if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
494: /* look for type on command line */
495: if (!MatColoringRegisterAllCalled) {MatColoringRegisterAll(PETSC_NULL);}
496: PetscOptionsGetString(((PetscObject)mat)->prefix,"-mat_coloring_type",tname,256,&flag);
497: if (flag) { type = tname; }
499: PetscObjectGetComm((PetscObject)mat,&comm);
500: PetscFListFind(MatColoringList,comm, type,PETSC_TRUE,(void (**)(void)) &r);
501: if (!r) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);
503: PetscLogEventBegin(MAT_GetColoring,mat,0,0,0);
504: (*r)(mat,type,iscoloring);
505: PetscLogEventEnd(MAT_GetColoring,mat,0,0,0);
507: PetscInfo1(mat,"Number of colors %d\n",(*iscoloring)->n);
508: flag = PETSC_FALSE;
509: PetscOptionsGetBool(PETSC_NULL,"-mat_coloring_view",&flag,PETSC_NULL);
510: if (flag) {
511: PetscViewer viewer;
512: PetscViewerASCIIGetStdout((*iscoloring)->comm,&viewer);
513: ISColoringView(*iscoloring,viewer);
514: }
515: return(0);
516: }
517: