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: