Actual source code: color.c
1: #define PETSCMAT_DLL
3: /*
4: Routines that call the kernel minpack coloring subroutines
5: */
7: #include include/private/matimpl.h
8: #include src/mat/color/color.h
10: /*
11: MatFDColoringDegreeSequence_Minpack - Calls the MINPACK routine seqr() that
12: computes the degree sequence required by MINPACK coloring routines.
13: */
16: PetscErrorCode MatFDColoringDegreeSequence_Minpack(PetscInt m,PetscInt *cja, PetscInt *cia, PetscInt *rja, PetscInt *ria, PetscInt **seq)
17: {
18: PetscInt *work;
22: PetscMalloc(m*sizeof(PetscInt),&work);
23: PetscMalloc(m*sizeof(PetscInt),seq);
25: MINPACKdegr(&m,cja,cia,rja,ria,*seq,work);
27: PetscFree(work);
28: return(0);
29: }
31: /*
32: MatFDColoringMinimumNumberofColors_Private - For a given sparse
33: matrix computes the minimum number of colors needed.
35: */
38: PetscErrorCode MatFDColoringMinimumNumberofColors_Private(PetscInt m,PetscInt *ia,PetscInt *minc)
39: {
40: PetscInt i,c = 0;
43: for (i=0; i<m; i++) {
44: c = PetscMax(c,ia[i+1]-ia[i]);
45: }
46: *minc = c;
47: return(0);
48: }
51: /* ----------------------------------------------------------------------------*/
52: /*
53: MatFDColoringSL_Minpack - Uses the smallest-last (SL) coloring of minpack
54: */
57: PetscErrorCode MatFDColoringSL_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
58: {
60: PetscInt *list,*work,clique,*ria,*rja,*cia,*cja,*seq,*coloring,n;
61: PetscInt ncolors,i;
62: PetscTruth done;
65: MatGetRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
66: MatGetColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
67: if (!done) SETERRQ(PETSC_ERR_SUP,"Ordering requires IJ");
69: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
71: PetscMalloc(5*n*sizeof(PetscInt),&list);
72: work = list + n;
74: MINPACKslo(&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
76: PetscMalloc(n*sizeof(PetscInt),&coloring);
77: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
79: PetscFree(list);
80: PetscFree(seq);
81: MatRestoreRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
82: MatRestoreColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
84: /* shift coloring numbers to start at zero and shorten */
85: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_ERR_SUP,"Maximum color size exceeded");
86: {
87: ISColoringValue *s = (ISColoringValue*) coloring;
88: for (i=0; i<n; i++) {
89: s[i] = (ISColoringValue) (coloring[i]-1);
90: }
91: MatColoringPatch(mat,ncolors,n,s,iscoloring);
92: }
93: return(0);
94: }
98: /* ----------------------------------------------------------------------------*/
99: /*
100: MatFDColoringLF_Minpack -
101: */
104: PetscErrorCode MatFDColoringLF_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
105: {
107: PetscInt *list,*work,*ria,*rja,*cia,*cja,*seq,*coloring,n;
108: PetscInt n1, none,ncolors,i;
109: PetscTruth done;
112: MatGetRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
113: MatGetColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
114: if (!done) SETERRQ(PETSC_ERR_SUP,"Ordering requires IJ");
116: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
118: PetscMalloc(5*n*sizeof(PetscInt),&list);
119: work = list + n;
121: n1 = n - 1;
122: none = -1;
123: MINPACKnumsrt(&n,&n1,seq,&none,list,work+2*n,work+n);
124: PetscMalloc(n*sizeof(PetscInt),&coloring);
125: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
127: PetscFree(list);
128: PetscFree(seq);
130: MatRestoreRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
131: MatRestoreColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
133: /* shift coloring numbers to start at zero and shorten */
134: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_ERR_SUP,"Maximum color size exceeded");
135: {
136: ISColoringValue *s = (ISColoringValue*) coloring;
137: for (i=0; i<n; i++) {
138: s[i] = (ISColoringValue) (coloring[i]-1);
139: }
140: MatColoringPatch(mat,ncolors,n,s,iscoloring);
141: }
142: return(0);
143: }
147: /* ----------------------------------------------------------------------------*/
148: /*
149: MatFDColoringID_Minpack -
150: */
153: PetscErrorCode MatFDColoringID_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
154: {
156: PetscInt *list,*work,clique,*ria,*rja,*cia,*cja,*seq,*coloring,n;
157: PetscInt ncolors,i;
158: PetscTruth done;
161: MatGetRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
162: MatGetColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
163: if (!done) SETERRQ(PETSC_ERR_SUP,"Ordering requires IJ");
165: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
167: PetscMalloc(5*n*sizeof(PetscInt),&list);
168: work = list + n;
170: MINPACKido(&n,&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
172: PetscMalloc(n*sizeof(PetscInt),&coloring);
173: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
175: PetscFree(list);
176: PetscFree(seq);
178: MatRestoreRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
179: MatRestoreColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
181: /* shift coloring numbers to start at zero and shorten */
182: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_ERR_SUP,"Maximum color size exceeded");
183: {
184: ISColoringValue *s = (ISColoringValue*) coloring;
185: for (i=0; i<n; i++) {
186: s[i] = (ISColoringValue) (coloring[i]-1);
187: }
188: MatColoringPatch(mat,ncolors,n,s,iscoloring);
189: }
190: return(0);
191: }
195: /*
196: Simplest coloring, each column of the matrix gets its own unique color.
197: */
200: PetscErrorCode MatColoring_Natural(Mat mat,const MatColoringType color, ISColoring *iscoloring)
201: {
202: PetscErrorCode ierr;
203: PetscInt start,end,i;
204: ISColoringValue *colors;
205: MPI_Comm comm;
208: MatGetOwnershipRange(mat,&start,&end);
209: PetscObjectGetComm((PetscObject)mat,&comm);
210: PetscMalloc((end-start+1)*sizeof(PetscInt),&colors);
211: for (i=start; i<end; i++) {
212: colors[i-start] = i;
213: }
214: ISColoringCreate(comm,mat->cmap.N,end-start,colors,iscoloring);
216: return(0);
217: }
219:
220: /* ===========================================================================================*/
222: #include petscsys.h
224: PetscFList MatColoringList = 0;
225: PetscTruth MatColoringRegisterAllCalled = PETSC_FALSE;
229: PetscErrorCode MatColoringRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(Mat,const MatColoringType,ISColoring*))
230: {
232: char fullname[PETSC_MAX_PATH_LEN];
235: PetscFListConcat(path,name,fullname);
236: PetscFListAdd(&MatColoringList,sname,fullname,(void (*)(void))function);
237: return(0);
238: }
242: /*@C
243: MatColoringRegisterDestroy - Frees the list of coloringing routines.
245: Not Collective
247: Level: developer
249: .keywords: matrix, register, destroy
251: .seealso: MatColoringRegisterDynamic(), MatColoringRegisterAll()
252: @*/
253: PetscErrorCode MatColoringRegisterDestroy(void)
254: {
258: if (MatColoringList) {
259: PetscFListDestroy(&MatColoringList);
260: MatColoringList = 0;
261: }
262: return(0);
263: }
267: /*@C
268: MatGetColoring - Gets a coloring for a matrix to reduce the number of function evaluations
269: needed to compute a sparse Jacobian via differencing.
271: Collective on Mat
273: Input Parameters:
274: . mat - the matrix
275: . type - type of coloring, one of the following:
276: $ MATCOLORING_NATURAL - natural (one color for each column, very slow)
277: $ MATCOLORING_SL - smallest-last
278: $ MATCOLORING_LF - largest-first
279: $ MATCOLORING_ID - incidence-degree
281: Output Parameters:
282: . iscoloring - the coloring
284: Options Database Keys:
285: To specify the coloring through the options database, use one of
286: the following
287: $ -mat_coloring_type natural, -mat_coloring_type sl, -mat_coloring_type lf,
288: $ -mat_coloring_type id
289: To see the coloring use
290: $ -mat_coloring_view
292: Level: intermediate
294: Notes:
295: These compute the graph coloring of the graph of A^{T}A. The coloring used
296: for efficient (parallel or thread based) triangular solves etc is NOT yet
297: available.
299: The user can define additional colorings; see MatColoringRegisterDynamic().
301: The sequential colorings SL, LF, and ID are obtained via the Minpack software that was
302: converted to C using f2c.
304: .keywords: matrix, get, coloring
306: .seealso: MatGetColoringTypeFromOptions(), MatColoringRegisterDynamic(), MatFDColoringCreate(),
307: SNESDefaultComputeJacobianColor()
308: @*/
309: PetscErrorCode MatGetColoring(Mat mat,const MatColoringType type,ISColoring *iscoloring)
310: {
311: PetscTruth flag;
312: PetscErrorCode ierr,(*r)(Mat,const MatColoringType,ISColoring *);
313: char tname[PETSC_MAX_PATH_LEN];
318: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
319: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
320: if (!MatColoringRegisterAllCalled) {
321: MatColoringRegisterAll(PETSC_NULL);
322: }
323:
324: /* look for type on command line */
325: PetscOptionsGetString(mat->prefix,"-mat_coloring_type",tname,256,&flag);
326: if (flag) {
327: type = tname;
328: }
331: PetscFListFind(mat->comm, MatColoringList, type,(void (**)(void)) &r);
332: if (!r) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);}
333: (*r)(mat,type,iscoloring);
336: PetscInfo1(mat,"Number of colors %d\n",(*iscoloring)->n);
337: PetscOptionsHasName(PETSC_NULL,"-mat_coloring_view",&flag);
338: if (flag) {
339: PetscViewer viewer;
340: PetscViewerASCIIGetStdout((*iscoloring)->comm,&viewer);
341: ISColoringView(*iscoloring,viewer);
342: }
343: return(0);
344: }
345: