Actual source code: color.c
1: /*$Id: color.c,v 1.60 2001/06/21 21:17:33 bsmith Exp $*/
2:
3: /*
4: Routines that call the kernel minpack coloring subroutines
5: */
7: #include src/mat/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: */
14: #undef __FUNCT__
16: int MatFDColoringDegreeSequence_Minpack(int m,int *cja, int *cia, int *rja, int *ria, int **seq)
17: {
18: int *work;
22: PetscMalloc(m*sizeof(int),&work);
23: PetscMalloc(m*sizeof(int),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: */
36: #undef __FUNCT__
38: int MatFDColoringMinimumNumberofColors_Private(int m,int *ia,int *minc)
39: {
40: int 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: }
50: EXTERN_C_BEGIN
51: /* ----------------------------------------------------------------------------*/
52: /*
53: MatFDColoringSL_Minpack - Uses the smallest-last (SL) coloring of minpack
54: */
55: #undef __FUNCT__
57: int MatFDColoringSL_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
58: {
59: int *list,*work,clique,ierr,*ria,*rja,*cia,*cja,*seq,*coloring,n;
60: int ncolors,i;
61: PetscTruth done;
64: MatGetRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
65: MatGetColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
66: if (!done) SETERRQ(PETSC_ERR_SUP,"Ordering requires IJ");
68: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
70: PetscMalloc(5*n*sizeof(int),&list);
71: work = list + n;
73: MINPACKslo(&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
75: PetscMalloc(n*sizeof(int),&coloring);
76: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
78: PetscFree(list);
79: PetscFree(seq);
80: MatRestoreRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
81: MatRestoreColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
83: /* shift coloring numbers to start at zero */
84: for (i=0; i<n; i++) coloring[i]--;
85: MatColoringPatch(mat,n,ncolors,coloring,iscoloring);
86: return(0);
87: }
88: EXTERN_C_END
90: EXTERN_C_BEGIN
91: /* ----------------------------------------------------------------------------*/
92: /*
93: MatFDColoringLF_Minpack -
94: */
95: #undef __FUNCT__
97: int MatFDColoringLF_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
98: {
99: int *list,*work,ierr,*ria,*rja,*cia,*cja,*seq,*coloring,n;
100: int n1, none,ncolors,i;
101: PetscTruth done;
104: MatGetRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
105: MatGetColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
106: if (!done) SETERRQ(PETSC_ERR_SUP,"Ordering requires IJ");
108: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
110: PetscMalloc(5*n*sizeof(int),&list);
111: work = list + n;
113: n1 = n - 1;
114: none = -1;
115: MINPACKnumsrt(&n,&n1,seq,&none,list,work+2*n,work+n);
116: PetscMalloc(n*sizeof(int),&coloring);
117: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
119: PetscFree(list);
120: PetscFree(seq);
122: MatRestoreRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
123: MatRestoreColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
125: /* shift coloring numbers to start at zero */
126: for (i=0; i<n; i++) coloring[i]--;
127: MatColoringPatch(mat,n,ncolors,coloring,iscoloring);
128: return(0);
129: }
130: EXTERN_C_END
132: EXTERN_C_BEGIN
133: /* ----------------------------------------------------------------------------*/
134: /*
135: MatFDColoringID_Minpack -
136: */
137: #undef __FUNCT__
139: int MatFDColoringID_Minpack(Mat mat,MatColoringType name,ISColoring *iscoloring)
140: {
141: int *list,*work,clique,ierr,*ria,*rja,*cia,*cja,*seq,*coloring,n;
142: int ncolors,i;
143: PetscTruth done;
146: MatGetRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
147: MatGetColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
148: if (!done) SETERRQ(PETSC_ERR_SUP,"Ordering requires IJ");
150: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
152: PetscMalloc(5*n*sizeof(int),&list);
153: work = list + n;
155: MINPACKido(&n,&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
157: PetscMalloc(n*sizeof(int),&coloring);
158: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
160: PetscFree(list);
161: PetscFree(seq);
163: MatRestoreRowIJ(mat,1,PETSC_FALSE,&n,&ria,&rja,&done);
164: MatRestoreColumnIJ(mat,1,PETSC_FALSE,&n,&cia,&cja,&done);
166: /* shift coloring numbers to start at zero */
167: for (i=0; i<n; i++) coloring[i]--;
168: MatColoringPatch(mat,n,ncolors,coloring,iscoloring);
169: return(0);
170: }
171: EXTERN_C_END
173: EXTERN_C_BEGIN
174: /*
175: Simplest coloring, each column of the matrix gets its own unique color.
176: */
177: #undef __FUNCT__
179: int MatColoring_Natural(Mat mat,MatColoringType color, ISColoring *iscoloring)
180: {
181: int start,end,ierr,i,*colors;
182: MPI_Comm comm;
185: MatGetOwnershipRange(mat,&start,&end);
186: PetscObjectGetComm((PetscObject)mat,&comm);
187: PetscMalloc((end-start+1)*sizeof(int),&colors);
188: for (i=start; i<end; i++) {
189: colors[i-start] = i;
190: }
191: ISColoringCreate(comm,end-start,colors,iscoloring);
193: return(0);
194: }
195: EXTERN_C_END
196:
197: /* ===========================================================================================*/
199: #include petscsys.h
201: PetscFList MatColoringList = 0;
202: PetscTruth MatColoringRegisterAllCalled = PETSC_FALSE;
204: /*MC
205: MatColoringRegisterDynamic - Adds a new sparse matrix coloring to the
206: matrix package.
208: Synopsis:
209: int MatColoringRegisterDynamic(char *name_coloring,char *path,char *name_create,int (*routine_create)(MatColoring))
211: Not Collective
213: Input Parameters:
214: + sname - name of Coloring (for example MATCOLORING_SL)
215: . path - location of library where creation routine is
216: . name - name of function that creates the Coloring type, a string
217: - function - function pointer that creates the coloring
219: Level: developer
221: If dynamic libraries are used, then the fourth input argument (function)
222: is ignored.
224: Sample usage:
225: .vb
226: MatColoringRegisterDynamic("my_color",/home/username/my_lib/lib/libO/solaris/mylib.a,
227: "MyColor",MyColor);
228: .ve
230: Then, your partitioner can be chosen with the procedural interface via
231: $ MatColoringSetType(part,"my_color")
232: or at runtime via the option
233: $ -mat_coloring_type my_color
235: $PETSC_ARCH and $BOPT occuring in pathname will be replaced with appropriate values.
237: .keywords: matrix, Coloring, register
239: .seealso: MatColoringRegisterDestroy(), MatColoringRegisterAll()
240: M*/
242: #undef __FUNCT__
244: int MatColoringRegister(char *sname,char *path,char *name,int (*function)(Mat,MatColoringType,ISColoring*))
245: {
246: int ierr;
247: char fullname[256];
250: PetscFListConcat(path,name,fullname);
251: PetscFListAdd(&MatColoringList,sname,fullname,(void (*)(void))function);
252: return(0);
253: }
255: #undef __FUNCT__
257: /*@C
258: MatColoringRegisterDestroy - Frees the list of coloringing routines.
260: Not Collective
262: Level: developer
264: .keywords: matrix, register, destroy
266: .seealso: MatColoringRegisterDynamic(), MatColoringRegisterAll()
267: @*/
268: int MatColoringRegisterDestroy(void)
269: {
273: if (MatColoringList) {
274: PetscFListDestroy(&MatColoringList);
275: MatColoringList = 0;
276: }
277: return(0);
278: }
280: EXTERN int MatAdjustForInodes(Mat,IS *,IS *);
282: #undef __FUNCT__
284: /*@C
285: MatGetColoring - Gets a coloring for a matrix to reduce the number of function evaluations
286: needed to compute a sparse Jacobian via differencing.
288: Collective on Mat
290: Input Parameters:
291: . mat - the matrix
292: . type - type of coloring, one of the following:
293: $ MATCOLORING_NATURAL - natural (one color for each column, very slow)
294: $ MATCOLORING_SL - smallest-last
295: $ MATCOLORING_LF - largest-first
296: $ MATCOLORING_ID - incidence-degree
298: Output Parameters:
299: . iscoloring - the coloring
301: Options Database Keys:
302: To specify the coloring through the options database, use one of
303: the following
304: $ -mat_coloring_type natural, -mat_coloring_type sl, -mat_coloring_type lf,
305: $ -mat_coloring_type id
306: To see the coloring use
307: $ -mat_coloring_view
309: Level: intermediate
311: Notes:
312: These compute the graph coloring of the graph of A^{T}A. The coloring used
313: for efficient (parallel or thread based) triangular solves etc is NOT yet
314: available.
316: The user can define additional colorings; see MatColoringRegisterDynamic().
318: The sequential colorings SL, LF, and ID are obtained via the Minpack software that was
319: converted to C using f2c.
321: .keywords: matrix, get, coloring
323: .seealso: MatGetColoringTypeFromOptions(), MatColoringRegisterDynamic(), MatFDColoringCreate(),
324: SNESDefaultComputeJacobianColor()
325: @*/
326: int MatGetColoring(Mat mat,MatColoringType type,ISColoring *iscoloring)
327: {
328: PetscTruth flag;
329: int ierr,(*r)(Mat,MatColoringType,ISColoring *);
330: char tname[256];
334: if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
335: if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
336: if (!MatColoringRegisterAllCalled) {
337: MatColoringRegisterAll(PETSC_NULL);
338: }
339:
340: /* look for type on command line */
341: PetscOptionsGetString(mat->prefix,"-mat_coloring_type",tname,256,&flag);
342: if (flag) {
343: type = tname;
344: }
346: PetscLogEventBegin(MAT_GetColoring,mat,0,0,0);
347: PetscFListFind(mat->comm, MatColoringList, type,(void (**)(void)) &r);
348: if (!r) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);}
349: (*r)(mat,type,iscoloring);
350: PetscLogEventEnd(MAT_GetColoring,mat,0,0,0);
352: PetscLogInfo((PetscObject)mat,"MatGetColoring:Number of colors %dn",(*iscoloring)->n);
353: PetscOptionsHasName(PETSC_NULL,"-mat_coloring_view",&flag);
354: if (flag) {
355: ISColoringView(*iscoloring,PETSC_VIEWER_STDOUT_((*iscoloring)->comm));
356: }
357: return(0);
358: }
359: