Actual source code: color.c
petsc-3.10.5 2019-03-28
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:
132: Supports only distance two colorings (for computation of Jacobians)
134: This is a sequential algorithm
136: References:
137: . 1. - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
138: pp. 187-209, 1983.
140: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
141: M*/
143: PETSC_EXTERN PetscErrorCode MatColoringCreate_SL(MatColoring mc)
144: {
146: mc->dist = 2;
147: mc->data = NULL;
148: mc->ops->apply = MatColoringApply_SL;
149: mc->ops->view = NULL;
150: mc->ops->destroy = NULL;
151: mc->ops->setfromoptions = NULL;
152: return(0);
153: }
155: static PetscErrorCode MatColoringApply_LF(MatColoring mc,ISColoring *iscoloring)
156: {
157: PetscErrorCode ierr;
158: PetscInt *list,*work,*seq,*coloring,n;
159: const PetscInt *ria,*rja,*cia,*cja;
160: PetscInt n1, none,ncolors,i;
161: PetscBool done;
162: Mat mat = mc->mat;
163: Mat mat_seq = mat;
164: PetscMPIInt size;
165: MPI_Comm comm;
166: ISColoring iscoloring_seq;
167: PetscInt bs = 1,rstart,rend,N_loc,nc;
168: ISColoringValue *colors_loc;
169: PetscBool flg1,flg2;
172: if (mc->dist != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"LF may only do distance 2 coloring");
173: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
174: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
175: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
176: if (flg1 || flg2) {
177: MatGetBlockSize(mat,&bs);
178: }
180: PetscObjectGetComm((PetscObject)mat,&comm);
181: MPI_Comm_size(comm,&size);
182: if (size > 1) {
183: /* create a sequential iscoloring on all processors */
184: MatGetSeqNonzeroStructure(mat,&mat_seq);
185: }
187: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
188: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
189: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
191: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
193: PetscMalloc2(n,&list,4*n,&work);
195: n1 = n - 1;
196: none = -1;
197: MINPACKnumsrt(&n,&n1,seq,&none,list,work+2*n,work+n);
198: PetscMalloc1(n,&coloring);
199: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
201: PetscFree2(list,work);
202: PetscFree(seq);
204: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
205: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
207: /* shift coloring numbers to start at zero and shorten */
208: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
209: {
210: ISColoringValue *s = (ISColoringValue*) coloring;
211: for (i=0; i<n; i++) s[i] = (ISColoringValue) (coloring[i]-1);
212: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
213: }
215: if (size > 1) {
216: MatDestroySeqNonzeroStructure(&mat_seq);
218: /* convert iscoloring_seq to a parallel iscoloring */
219: iscoloring_seq = *iscoloring;
220: rstart = mat->rmap->rstart/bs;
221: rend = mat->rmap->rend/bs;
222: N_loc = rend - rstart; /* number of local nodes */
224: /* get local colors for each local node */
225: PetscMalloc1(N_loc+1,&colors_loc);
226: for (i=rstart; i<rend; i++) colors_loc[i-rstart] = iscoloring_seq->colors[i];
228: /* create a parallel iscoloring */
229: nc = iscoloring_seq->n;
230: ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
231: ISColoringDestroy(&iscoloring_seq);
232: }
233: return(0);
234: }
236: /*MC
237: MATCOLORINGLF - implements the LF (largest first) coloring routine
239: Level: beginner
241: Notes:
242: Supports only distance two colorings (for computation of Jacobians)
244: This is a sequential algorithm
246: References:
247: . 1. - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
248: pp. 187-209, 1983.
250: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
251: M*/
253: PETSC_EXTERN PetscErrorCode MatColoringCreate_LF(MatColoring mc)
254: {
256: mc->dist = 2;
257: mc->data = NULL;
258: mc->ops->apply = MatColoringApply_LF;
259: mc->ops->view = NULL;
260: mc->ops->destroy = NULL;
261: mc->ops->setfromoptions = NULL;
262: return(0);
263: }
265: static PetscErrorCode MatColoringApply_ID(MatColoring mc,ISColoring *iscoloring)
266: {
267: PetscErrorCode ierr;
268: PetscInt *list,*work,clique,*seq,*coloring,n;
269: const PetscInt *ria,*rja,*cia,*cja;
270: PetscInt ncolors,i;
271: PetscBool done;
272: Mat mat = mc->mat;
273: Mat mat_seq = mat;
274: PetscMPIInt size;
275: MPI_Comm comm;
276: ISColoring iscoloring_seq;
277: PetscInt bs = 1,rstart,rend,N_loc,nc;
278: ISColoringValue *colors_loc;
279: PetscBool flg1,flg2;
282: if (mc->dist != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IDO may only do distance 2 coloring");
283: /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
284: PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
285: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);
286: if (flg1 || flg2) {
287: MatGetBlockSize(mat,&bs);
288: }
290: PetscObjectGetComm((PetscObject)mat,&comm);
291: MPI_Comm_size(comm,&size);
292: if (size > 1) {
293: /* create a sequential iscoloring on all processors */
294: MatGetSeqNonzeroStructure(mat,&mat_seq);
295: }
297: MatGetRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&ria,&rja,&done);
298: MatGetColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,&n,&cia,&cja,&done);
299: if (!done) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Ordering requires IJ");
301: MatFDColoringDegreeSequence_Minpack(n,cja,cia,rja,ria,&seq);
303: PetscMalloc2(n,&list,4*n,&work);
305: MINPACKido(&n,&n,cja,cia,rja,ria,seq,list,&clique,work,work+n,work+2*n,work+3*n);
307: PetscMalloc1(n,&coloring);
308: MINPACKseq(&n,cja,cia,rja,ria,list,coloring,&ncolors,work);
310: PetscFree2(list,work);
311: PetscFree(seq);
313: MatRestoreRowIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&ria,&rja,&done);
314: MatRestoreColumnIJ(mat_seq,1,PETSC_FALSE,PETSC_TRUE,NULL,&cia,&cja,&done);
316: /* shift coloring numbers to start at zero and shorten */
317: if (ncolors > IS_COLORING_MAX-1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Maximum color size exceeded");
318: {
319: ISColoringValue *s = (ISColoringValue*) coloring;
320: for (i=0; i<n; i++) {
321: s[i] = (ISColoringValue) (coloring[i]-1);
322: }
323: MatColoringPatch(mat_seq,ncolors,n,s,iscoloring);
324: }
326: if (size > 1) {
327: MatDestroySeqNonzeroStructure(&mat_seq);
329: /* convert iscoloring_seq to a parallel iscoloring */
330: iscoloring_seq = *iscoloring;
331: rstart = mat->rmap->rstart/bs;
332: rend = mat->rmap->rend/bs;
333: N_loc = rend - rstart; /* number of local nodes */
335: /* get local colors for each local node */
336: PetscMalloc1(N_loc+1,&colors_loc);
337: for (i=rstart; i<rend; i++) {
338: colors_loc[i-rstart] = iscoloring_seq->colors[i];
339: }
340: /* create a parallel iscoloring */
341: nc = iscoloring_seq->n;
342: ISColoringCreate(comm,nc,N_loc,colors_loc,PETSC_OWN_POINTER,iscoloring);
343: ISColoringDestroy(&iscoloring_seq);
344: }
345: return(0);
346: }
348: /*MC
349: MATCOLORINGID - implements the ID (incidence degree) coloring routine
351: Level: beginner
353: Notes:
354: Supports only distance two colorings (for computation of Jacobians)
356: This is a sequential algorithm
358: References:
359: . 1. - TF Coleman and J More, "Estimation of sparse Jacobian matrices and graph coloring," SIAM Journal on Numerical Analysis, vol. 20, no. 1,
360: pp. 187-209, 1983.
362: .seealso: MatColoringCreate(), MatColoring, MatColoringSetType(), MATCOLORINGGREEDY, MatColoringType
363: M*/
365: PETSC_EXTERN PetscErrorCode MatColoringCreate_ID(MatColoring mc)
366: {
368: mc->dist = 2;
369: mc->data = NULL;
370: mc->ops->apply = MatColoringApply_ID;
371: mc->ops->view = NULL;
372: mc->ops->destroy = NULL;
373: mc->ops->setfromoptions = NULL;
374: return(0);
375: }