Actual source code: color.c

      1: #define PETSCMAT_DLL

      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: */
     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 PETSCMAT_DLLEXPORT 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,n,ncolors,s,iscoloring);
     92:   }
     93:   return(0);
     94: }

     98: /* ----------------------------------------------------------------------------*/
     99: /*
    100:     MatFDColoringLF_Minpack -
    101: */
    104: PetscErrorCode PETSCMAT_DLLEXPORT 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,n,ncolors,s,iscoloring);
    141:   }
    142:   return(0);
    143: }

    147: /* ----------------------------------------------------------------------------*/
    148: /*
    149:     MatFDColoringID_Minpack -
    150: */
    153: PetscErrorCode PETSCMAT_DLLEXPORT 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,n,ncolors,s,iscoloring);
    189:   }
    190:   return(0);
    191: }

    195: /*
    196:    Simplest coloring, each column of the matrix gets its own unique color.
    197: */
    200: PetscErrorCode PETSCMAT_DLLEXPORT 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,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 PETSCMAT_DLLEXPORT 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 PETSCMAT_DLLEXPORT 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 PETSCMAT_DLLEXPORT 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:   }

    330:   PetscLogEventBegin(MAT_GetColoring,mat,0,0,0);
    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);
    334:   PetscLogEventEnd(MAT_GetColoring,mat,0,0,0);

    336:   PetscLogInfo(((PetscObject)mat,"MatGetColoring:Number of colors %d\n",(*iscoloring)->n));
    337:   PetscOptionsHasName(PETSC_NULL,"-mat_coloring_view",&flag);
    338:   if (flag) {
    339:     ISColoringView(*iscoloring,PETSC_VIEWER_STDOUT_((*iscoloring)->comm));
    340:   }
    341:   return(0);
    342: }
    343: