Actual source code: color.c

      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:   PetscObjectBaseTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
     64:   PetscObjectBaseTypeCompare((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:   PetscObjectBaseTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
    175:   PetscObjectBaseTypeCompare((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:   PetscObjectBaseTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);
    285:   PetscObjectBaseTypeCompare((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: }