Actual source code: color.c
petsc-3.9.4 2018-09-11
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;
20: PetscMalloc1(m,&work);
21: PetscMalloc1(m,seq);
23: MINPACKdegr(&m,cja,cia,rja,ria,*seq,work);
25: PetscFree(work);
26: return(0);
27: }
29: /*
30: MatFDColoringMinimumNumberofColors_Private - For a given sparse
31: matrix computes the minimum number of colors needed.
33: */
34: PetscErrorCode MatFDColoringMinimumNumberofColors_Private(PetscInt m,PetscInt *ia,PetscInt *minc)
35: {
36: PetscInt i,c = 0;
39: for (i=0; i<m; i++) c = PetscMax(c,ia[i+1]-ia[i]);
40: *minc = c;
41: return(0);
42: }
44: static PetscErrorCode MatColoringApply_SL(MatColoring mc,ISColoring *iscoloring)
45: {
46: PetscErrorCode ierr;
47: PetscInt *list,*work,clique,*seq,*coloring,n;
48: const PetscInt *ria,*rja,*cia,*cja;
49: PetscInt ncolors,i;
50: PetscBool done;
51: Mat mat = mc->mat;
52: Mat mat_seq = mat;
53: PetscMPIInt size;
54: MPI_Comm comm;
55: ISColoring iscoloring_seq;
56: PetscInt bs = 1,rstart,rend,N_loc,nc;
57: ISColoringValue *colors_loc;
58: PetscBool flg1,flg2;
61: if (mc->dist != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SL may only do distance 2 coloring");
62: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
63: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
64: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
65: if (flg1 || flg2) {
66: MatGetBlockSize(mat,&bs);
67: }
69: PetscObjectGetComm((PetscObject)mat,&comm);
70: MPI_Comm_size(comm,&size);
71: if (size > 1) {
72: /* create a sequential iscoloring on all processors */
73: MatGetSeqNonzeroStructure(mat,&mat_seq);
74: }
76: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
77: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
78: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
80: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
82: PetscMalloc2(n,&list,4*n,&work);
84: MINPACKslo(&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
86: PetscMalloc1(n,&coloring);
87: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
89: PetscFree2(list,work);
90: PetscFree(seq);
91: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
92: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
94: /* shift coloring numbers to start at zero and shorten */
95: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
96: {
97: ISColoringValue *s = (ISColoringValue*) coloring;
98: for (i=0; i<n; i++) {
99: s[i] = (ISColoringValue) (coloring[i]-1);
100: }
101: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
102: }
104: if (size > 1) {
105: MatDestroySeqNonzeroStructure(&mat_seq);
107: /* convert iscoloring_seq to a parallel iscoloring */
108: iscoloring_seq = *iscoloring;
109: rstart = mat->rmap->rstart/bs;
110: rend = mat->rmap->rend/bs;
111: N_loc = rend - rstart; /* number of local nodes */
113: /* get local colors for each local node */
114: PetscMalloc1(N_loc+1,&colors_loc);
115: for (i=rstart; i<rend; i++) {
116: colors_loc[i-rstart] = iscoloring_seq->colors[i];
117: }
118: /* create a parallel iscoloring */
119: nc = iscoloring_seq->n;
120: ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
121: ISColoringDestroy(&iscoloring_seq);
122: }
123: return(0);
124: }
126: /*MC
127: MATCOLORINGSL - implements the SL (smallest last) coloring routine
129: Level: beginner
131: Notes: Supports only distance two colorings (for computation of Jacobians)
133: This is a sequential algorithm
135: References:
136: . 1. - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
137: pp. 187-209, 1983.
139: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
140: M*/
142: PETSC_EXTERN PetscErrorCode MatColoringCreate_SL(MatColoring mc)
143: {
145: mc->dist = 2;
146: mc->data = NULL;
147: mc->ops->apply = MatColoringApply_SL;
148: mc->ops->view = NULL;
149: mc->ops->destroy = NULL;
150: mc->ops->setfromoptions = NULL;
151: return(0);
152: }
154: static PetscErrorCode MatColoringApply_LF(MatColoring mc,ISColoring *iscoloring)
155: {
156: PetscErrorCode ierr;
157: PetscInt *list,*work,*seq,*coloring,n;
158: const PetscInt *ria,*rja,*cia,*cja;
159: PetscInt n1, none,ncolors,i;
160: PetscBool done;
161: Mat mat = mc->mat;
162: Mat mat_seq = mat;
163: PetscMPIInt size;
164: MPI_Comm comm;
165: ISColoring iscoloring_seq;
166: PetscInt bs = 1,rstart,rend,N_loc,nc;
167: ISColoringValue *colors_loc;
168: PetscBool flg1,flg2;
171: if (mc->dist != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"LF may only do distance 2 coloring");
172: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
173: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
174: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
175: if (flg1 || flg2) {
176: MatGetBlockSize(mat,&bs);
177: }
179: PetscObjectGetComm((PetscObject)mat,&comm);
180: MPI_Comm_size(comm,&size);
181: if (size > 1) {
182: /* create a sequential iscoloring on all processors */
183: MatGetSeqNonzeroStructure(mat,&mat_seq);
184: }
186: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
187: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
188: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
190: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
192: PetscMalloc2(n,&list,4*n,&work);
194: n1 = n - 1;
195: none = -1;
196: MINPACKnumsrt(&n,&n1,seq,&none,list,work+2*n,work+n);
197: PetscMalloc1(n,&coloring);
198: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
200: PetscFree2(list,work);
201: PetscFree(seq);
203: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
204: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
206: /* shift coloring numbers to start at zero and shorten */
207: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
208: {
209: ISColoringValue *s = (ISColoringValue*) coloring;
210: for (i=0; i<n; i++) s[i] = (ISColoringValue) (coloring[i]-1);
211: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
212: }
214: if (size > 1) {
215: MatDestroySeqNonzeroStructure(&mat_seq);
217: /* convert iscoloring_seq to a parallel iscoloring */
218: iscoloring_seq = *iscoloring;
219: rstart = mat->rmap->rstart/bs;
220: rend = mat->rmap->rend/bs;
221: N_loc = rend - rstart; /* number of local nodes */
223: /* get local colors for each local node */
224: PetscMalloc1(N_loc+1,&colors_loc);
225: for (i=rstart; i<rend; i++) colors_loc[i-rstart] = iscoloring_seq->colors[i];
227: /* create a parallel iscoloring */
228: nc = iscoloring_seq->n;
229: ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
230: ISColoringDestroy(&iscoloring_seq);
231: }
232: return(0);
233: }
235: /*MC
236: MATCOLORINGLF - implements the LF (largest first) coloring routine
238: Level: beginner
240: Notes: Supports only distance two colorings (for computation of Jacobians)
242: This is a sequential algorithm
244: References:
245: . 1. - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
246: pp. 187-209, 1983.
248: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
249: M*/
251: PETSC_EXTERN PetscErrorCode MatColoringCreate_LF(MatColoring mc)
252: {
254: mc->dist = 2;
255: mc->data = NULL;
256: mc->ops->apply = MatColoringApply_LF;
257: mc->ops->view = NULL;
258: mc->ops->destroy = NULL;
259: mc->ops->setfromoptions = NULL;
260: return(0);
261: }
263: static PetscErrorCode MatColoringApply_ID(MatColoring mc,ISColoring *iscoloring)
264: {
265: PetscErrorCode ierr;
266: PetscInt *list,*work,clique,*seq,*coloring,n;
267: const PetscInt *ria,*rja,*cia,*cja;
268: PetscInt ncolors,i;
269: PetscBool done;
270: Mat mat = mc->mat;
271: Mat mat_seq = mat;
272: PetscMPIInt size;
273: MPI_Comm comm;
274: ISColoring iscoloring_seq;
275: PetscInt bs = 1,rstart,rend,N_loc,nc;
276: ISColoringValue *colors_loc;
277: PetscBool flg1,flg2;
280: if (mc->dist != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IDO may only do distance 2 coloring");
281: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
282: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
283: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
284: if (flg1 || flg2) {
285: MatGetBlockSize(mat,&bs);
286: }
288: PetscObjectGetComm((PetscObject)mat,&comm);
289: MPI_Comm_size(comm,&size);
290: if (size > 1) {
291: /* create a sequential iscoloring on all processors */
292: MatGetSeqNonzeroStructure(mat,&mat_seq);
293: }
295: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
296: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
297: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
299: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
301: PetscMalloc2(n,&list,4*n,&work);
303: MINPACKido(&n,&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
305: PetscMalloc1(n,&coloring);
306: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
308: PetscFree2(list,work);
309: PetscFree(seq);
311: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
312: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
314: /* shift coloring numbers to start at zero and shorten */
315: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
316: {
317: ISColoringValue *s = (ISColoringValue*) coloring;
318: for (i=0; i<n; i++) {
319: s[i] = (ISColoringValue) (coloring[i]-1);
320: }
321: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
322: }
324: if (size > 1) {
325: MatDestroySeqNonzeroStructure(&mat_seq);
327: /* convert iscoloring_seq to a parallel iscoloring */
328: iscoloring_seq = *iscoloring;
329: rstart = mat->rmap->rstart/bs;
330: rend = mat->rmap->rend/bs;
331: N_loc = rend - rstart; /* number of local nodes */
333: /* get local colors for each local node */
334: PetscMalloc1(N_loc+1,&colors_loc);
335: for (i=rstart; i<rend; i++) {
336: colors_loc[i-rstart] = iscoloring_seq->colors[i];
337: }
338: /* create a parallel iscoloring */
339: nc = iscoloring_seq->n;
340: ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
341: ISColoringDestroy(&iscoloring_seq);
342: }
343: return(0);
344: }
346: /*MC
347: MATCOLORINGID - implements the ID (incidence degree) coloring routine
349: Level: beginner
351: Notes: Supports only distance two colorings (for computation of Jacobians)
353: This is a sequential algorithm
355: References:
356: . 1. - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
357: pp. 187-209, 1983.
359: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
360: M*/
362: PETSC_EXTERN PetscErrorCode MatColoringCreate_ID(MatColoring mc)
363: {
365: mc->dist = 2;
366: mc->data = NULL;
367: mc->ops->apply = MatColoringApply_ID;
368: mc->ops->view = NULL;
369: mc->ops->destroy = NULL;
370: mc->ops->setfromoptions = NULL;
371: return(0);
372: }