Actual source code: color.c
petsc-3.4.5 2014-06-29
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,const PetscInt *cja,const PetscInt *cia,const PetscInt *rja,const 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++) c = PetscMax(c,ia[i+1]-ia[i]);
43: *minc = c;
44: return(0);
45: }
47: /* ----------------------------------------------------------------------------*/
48: /*
49: MatGetColoring_SL_Minpack - Uses the smallest-last (SL) coloring of minpack
50: */
53: PETSC_EXTERN PetscErrorCode MatGetColoring_SL_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
54: {
55: PetscErrorCode ierr;
56: PetscInt *list,*work,clique,*seq,*coloring,n;
57: const PetscInt *ria,*rja,*cia,*cja;
58: PetscInt ncolors,i;
59: PetscBool done;
60: Mat mat_seq = mat;
61: PetscMPIInt size;
62: MPI_Comm comm;
63: ISColoring iscoloring_seq;
64: PetscInt bs = 1,rstart,rend,N_loc,nc;
65: ISColoringValue *colors_loc;
66: PetscBool flg1,flg2;
69: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
70: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
71: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
72: if (flg1 || flg2) {
73: MatGetBlockSize(mat,&bs);
74: }
76: PetscObjectGetComm((PetscObject)mat,&comm);
77: MPI_Comm_size(comm,&size);
78: if (size > 1) {
79: /* create a sequential iscoloring on all processors */
80: MatGetSeqNonzeroStructure(mat,&mat_seq);
81: }
83: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
84: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
85: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
87: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
89: PetscMalloc2(n,PetscInt,&list,4*n,PetscInt,&work);
91: MINPACKslo(&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
93: PetscMalloc(n*sizeof(PetscInt),&coloring);
94: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
96: PetscFree2(list,work);
97: PetscFree(seq);
98: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
99: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
101: /* shift coloring numbers to start at zero and shorten */
102: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
103: {
104: ISColoringValue *s = (ISColoringValue*) coloring;
105: for (i=0; i<n; i++) {
106: s[i] = (ISColoringValue) (coloring[i]-1);
107: }
108: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
109: }
111: if (size > 1) {
112: MatDestroySeqNonzeroStructure(&mat_seq);
114: /* convert iscoloring_seq to a parallel iscoloring */
115: iscoloring_seq = *iscoloring;
116: rstart = mat->rmap->rstart/bs;
117: rend = mat->rmap->rend/bs;
118: N_loc = rend - rstart; /* number of local nodes */
120: /* get local colors for each local node */
121: PetscMalloc((N_loc+1)*sizeof(ISColoringValue),&colors_loc);
122: for (i=rstart; i<rend; i++) {
123: colors_loc[i-rstart] = iscoloring_seq->colors[i];
124: }
125: /* create a parallel iscoloring */
126: nc = iscoloring_seq->n;
127: ISColoringCreate(comm,nc,N_loc,colors_loc,iscoloring);
128: ISColoringDestroy(&iscoloring_seq);
129: }
130: return(0);
131: }
133: /* ----------------------------------------------------------------------------*/
134: /*
135: MatGetColoring_LF_Minpack -
136: */
139: PETSC_EXTERN PetscErrorCode MatGetColoring_LF_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
140: {
141: PetscErrorCode ierr;
142: PetscInt *list,*work,*seq,*coloring,n;
143: const PetscInt *ria,*rja,*cia,*cja;
144: PetscInt n1, none,ncolors,i;
145: PetscBool done;
146: Mat mat_seq = mat;
147: PetscMPIInt size;
148: MPI_Comm comm;
149: ISColoring iscoloring_seq;
150: PetscInt bs = 1,rstart,rend,N_loc,nc;
151: ISColoringValue *colors_loc;
152: PetscBool flg1,flg2;
155: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
156: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
157: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
158: if (flg1 || flg2) {
159: MatGetBlockSize(mat,&bs);
160: }
162: PetscObjectGetComm((PetscObject)mat,&comm);
163: MPI_Comm_size(comm,&size);
164: if (size > 1) {
165: /* create a sequential iscoloring on all processors */
166: MatGetSeqNonzeroStructure(mat,&mat_seq);
167: }
169: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
170: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
171: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
173: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
175: PetscMalloc2(n,PetscInt,&list,4*n,PetscInt,&work);
177: n1 = n - 1;
178: none = -1;
179: MINPACKnumsrt(&n,&n1,seq,&none,list,work+2*n,work+n);
180: PetscMalloc(n*sizeof(PetscInt),&coloring);
181: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
183: PetscFree2(list,work);
184: PetscFree(seq);
186: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
187: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
189: /* shift coloring numbers to start at zero and shorten */
190: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
191: {
192: ISColoringValue *s = (ISColoringValue*) coloring;
193: for (i=0; i<n; i++) s[i] = (ISColoringValue) (coloring[i]-1);
194: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
195: }
197: if (size > 1) {
198: MatDestroySeqNonzeroStructure(&mat_seq);
200: /* convert iscoloring_seq to a parallel iscoloring */
201: iscoloring_seq = *iscoloring;
202: rstart = mat->rmap->rstart/bs;
203: rend = mat->rmap->rend/bs;
204: N_loc = rend - rstart; /* number of local nodes */
206: /* get local colors for each local node */
207: PetscMalloc((N_loc+1)*sizeof(ISColoringValue),&colors_loc);
208: for (i=rstart; i<rend; i++) colors_loc[i-rstart] = iscoloring_seq->colors[i];
210: /* create a parallel iscoloring */
211: nc = iscoloring_seq->n;
212: ISColoringCreate(comm,nc,N_loc,colors_loc,iscoloring);
213: ISColoringDestroy(&iscoloring_seq);
214: }
215: return(0);
216: }
218: /* ----------------------------------------------------------------------------*/
219: /*
220: MatGetColoring_ID_Minpack -
221: */
224: PETSC_EXTERN PetscErrorCode MatGetColoring_ID_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
225: {
226: PetscErrorCode ierr;
227: PetscInt *list,*work,clique,*seq,*coloring,n;
228: const PetscInt *ria,*rja,*cia,*cja;
229: PetscInt ncolors,i;
230: PetscBool done;
231: Mat mat_seq = mat;
232: PetscMPIInt size;
233: MPI_Comm comm;
234: ISColoring iscoloring_seq;
235: PetscInt bs = 1,rstart,rend,N_loc,nc;
236: ISColoringValue *colors_loc;
237: PetscBool flg1,flg2;
240: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
241: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
242: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
243: if (flg1 || flg2) {
244: MatGetBlockSize(mat,&bs);
245: }
247: PetscObjectGetComm((PetscObject)mat,&comm);
248: MPI_Comm_size(comm,&size);
249: if (size > 1) {
250: /* create a sequential iscoloring on all processors */
251: MatGetSeqNonzeroStructure(mat,&mat_seq);
252: }
254: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
255: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
256: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
258: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
260: PetscMalloc2(n,PetscInt,&list,4*n,PetscInt,&work);
262: MINPACKido(&n,&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
264: PetscMalloc(n*sizeof(PetscInt),&coloring);
265: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
267: PetscFree2(list,work);
268: PetscFree(seq);
270: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
271: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
273: /* shift coloring numbers to start at zero and shorten */
274: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
275: {
276: ISColoringValue *s = (ISColoringValue*) coloring;
277: for (i=0; i<n; i++) {
278: s[i] = (ISColoringValue) (coloring[i]-1);
279: }
280: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
281: }
283: if (size > 1) {
284: MatDestroySeqNonzeroStructure(&mat_seq);
286: /* convert iscoloring_seq to a parallel iscoloring */
287: iscoloring_seq = *iscoloring;
288: rstart = mat->rmap->rstart/bs;
289: rend = mat->rmap->rend/bs;
290: N_loc = rend - rstart; /* number of local nodes */
292: /* get local colors for each local node */
293: PetscMalloc((N_loc+1)*sizeof(ISColoringValue),&colors_loc);
294: for (i=rstart; i<rend; i++) {
295: colors_loc[i-rstart] = iscoloring_seq->colors[i];
296: }
297: /* create a parallel iscoloring */
298: nc = iscoloring_seq->n;
299: ISColoringCreate(comm,nc,N_loc,colors_loc,iscoloring);
300: ISColoringDestroy(&iscoloring_seq);
301: }
302: return(0);
303: }
305: /*
306: Simplest coloring, each column of the matrix gets its own unique color.
307: */
310: PETSC_EXTERN PetscErrorCode MatGetColoring_Natural(Mat mat,MatColoringType color, ISColoring *iscoloring)
311: {
312: PetscErrorCode ierr;
313: PetscInt start,end,i,bs = 1,n;
314: ISColoringValue *colors;
315: MPI_Comm comm;
316: PetscBool flg1,flg2;
317: Mat mat_seq = mat;
318: PetscMPIInt size;
319: ISColoring iscoloring_seq;
320: ISColoringValue *colors_loc;
321: PetscInt rstart,rend,N_loc,nc;
324: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
325: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
326: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
327: if (flg1 || flg2) {
328: MatGetBlockSize(mat,&bs);
329: }
331: PetscObjectGetComm((PetscObject)mat,&comm);
332: MPI_Comm_size(comm,&size);
333: if (size > 1) {
334: /* create a sequential iscoloring on all processors */
335: MatGetSeqNonzeroStructure(mat,&mat_seq);
336: }
338: MatGetSize(mat_seq,NULL,&n);
339: MatGetOwnershipRange(mat_seq,&start,&end);
340: n = n/bs;
341: if (n > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
343: start = start/bs;
344: end = end/bs;
345: PetscMalloc((end-start+1)*sizeof(PetscInt),&colors);
346: for (i=start; i<end; i++) {
347: colors[i-start] = (ISColoringValue)i;
348: }
349: ISColoringCreate(comm,n,end-start,colors,iscoloring);
351: if (size > 1) {
352: MatDestroySeqNonzeroStructure(&mat_seq);
354: /* convert iscoloring_seq to a parallel iscoloring */
355: iscoloring_seq = *iscoloring;
356: rstart = mat->rmap->rstart/bs;
357: rend = mat->rmap->rend/bs;
358: N_loc = rend - rstart; /* number of local nodes */
360: /* get local colors for each local node */
361: PetscMalloc((N_loc+1)*sizeof(ISColoringValue),&colors_loc);
362: for (i=rstart; i<rend; i++) {
363: colors_loc[i-rstart] = iscoloring_seq->colors[i];
364: }
365: /* create a parallel iscoloring */
366: nc = iscoloring_seq->n;
367: ISColoringCreate(comm,nc,N_loc,colors_loc,iscoloring);
368: ISColoringDestroy(&iscoloring_seq);
369: }
370: return(0);
371: }
373: /* ===========================================================================================*/
375: PetscFunctionList MatColoringList = 0;
376: PetscBool MatColoringRegisterAllCalled = PETSC_FALSE;
380: /*@C
381: MatColoringRegister - Adds a new sparse matrix coloring to the matrix package.
383: Not Collective
385: Input Parameters:
386: + sname - name of Coloring (for example MATCOLORINGSL)
387: - function - function pointer that creates the coloring
389: Level: developer
391: Sample usage:
392: .vb
393: MatColoringRegister("my_color",MyColor);
394: .ve
396: Then, your partitioner can be chosen with the procedural interface via
397: $ MatColoringSetType(part,"my_color")
398: or at runtime via the option
399: $ -mat_coloring_type my_color
401: .keywords: matrix, Coloring, register
403: .seealso: MatColoringRegisterDestroy(), MatColoringRegisterAll()
404: @*/
405: PetscErrorCode MatColoringRegister(const char sname[],PetscErrorCode (*function)(Mat,MatColoringType,ISColoring*))
406: {
410: PetscFunctionListAdd(&MatColoringList,sname,function);
411: return(0);
412: }
416: /*@C
417: MatGetColoring - Gets a coloring for a matrix, from its sparsity structure,
418: to reduce the number of function evaluations needed to compute a sparse Jacobian via differencing.
420: Collective on Mat
422: Input Parameters:
423: . mat - the matrix
424: . type - type of coloring, one of the following:
425: $ MATCOLORINGNATURAL - natural (one color for each column, very slow)
426: $ MATCOLORINGSL - smallest-last
427: $ MATCOLORINGLF - largest-first
428: $ MATCOLORINGID - incidence-degree
430: Output Parameters:
431: . iscoloring - the coloring
433: Options Database Keys:
434: To specify the coloring through the options database, use one of
435: the following
436: $ -mat_coloring_type natural, -mat_coloring_type sl, -mat_coloring_type lf,
437: $ -mat_coloring_type id
438: To see the coloring use
439: $ -mat_coloring_view
441: Level: intermediate
443: Notes:
444: $ A graph coloring C(A) is a division of vertices so that two vertices of the same color do not share any common edges.
445: $ A suitable coloring for a smoother is simply C(A).
446: $ 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.
447: $ This corresponds to C(A^{T} A). This is what MatGetColoring() computes.
449: The user can define additional colorings; see MatColoringRegister().
451: For parallel matrices currently converts to sequential matrix and uses the sequential coloring
452: on that.
454: The colorings SL, LF, and ID are obtained via the Minpack software that was
455: converted to C using f2c.
457: For BAIJ matrices this colors the blocks. The true number of colors would be block size times the number of colors
458: returned here.
460: References:
461: $ Thomas F. Coleman and Jorge J. More, Estimation of Sparse {J}acobian Matrices and Graph Coloring Problems,
462: $ SIAM Journal on Numerical Analysis, 1983, pages 187-209, volume 20
463: $ Jorge J. Mor\'{e} and Danny C. Sorenson and Burton S. Garbow and Kenneth E. Hillstrom, The {MINPACK} Project,
464: $ Sources and Development of Mathematical Software, Wayne R. Cowell editor, 1984, pages 88-111
466: .keywords: matrix, get, coloring
468: .seealso: MatGetColoringTypeFromOptions(), MatColoringRegister(), MatFDColoringCreate(),
469: SNESComputeJacobianDefaultColor()
470: @*/
471: PetscErrorCode MatGetColoring(Mat mat,MatColoringType type,ISColoring *iscoloring)
472: {
473: PetscBool flag;
474: PetscErrorCode ierr,(*r)(Mat,MatColoringType,ISColoring*);
475: char tname[PETSC_MAX_PATH_LEN];
476: MPI_Comm comm;
481: if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
482: if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
484: /* look for type on command line */
485: if (!MatColoringRegisterAllCalled) {MatColoringRegisterAll();}
486: PetscOptionsGetString(((PetscObject)mat)->prefix,"-mat_coloring_type",tname,256,&flag);
487: if (flag) type = tname;
489: PetscObjectGetComm((PetscObject)mat,&comm);
490: PetscFunctionListFind(MatColoringList,type,&r);
491: if (!r) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);
493: PetscLogEventBegin(MAT_GetColoring,mat,0,0,0);
494: (*r)(mat,type,iscoloring);
495: PetscLogEventEnd(MAT_GetColoring,mat,0,0,0);
497: PetscInfo1(mat,"Number of colors %d\n",(*iscoloring)->n);
498: flag = PETSC_FALSE;
499: PetscOptionsGetBool(NULL,"-mat_coloring_view",&flag,NULL);
500: if (flag) {
501: PetscViewer viewer;
502: PetscViewerASCIIGetStdout((*iscoloring)->comm,&viewer);
503: ISColoringView(*iscoloring,viewer);
504: }
505: return(0);
506: }