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: