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: }