Actual source code: fdmatrix.c

      2: /*
      3:    This is where the abstract matrix operations are defined that are
      4:   used for finite difference computations of Jacobians using coloring.
      5: */

      7: #include <petsc/private/matimpl.h>
      8: #include <petsc/private/isimpl.h>

     10: PetscErrorCode  MatFDColoringSetF(MatFDColoring fd,Vec F)
     11: {

     15:   if (F) {
     16:     VecCopy(F,fd->w1);
     17:     fd->fset = PETSC_TRUE;
     18:   } else {
     19:     fd->fset = PETSC_FALSE;
     20:   }
     21:   return(0);
     22: }

     24: #include <petscdraw.h>
     25: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw,void *Aa)
     26: {
     27:   MatFDColoring  fd = (MatFDColoring)Aa;
     29:   PetscInt       i,j,nz,row;
     30:   PetscReal      x,y;
     31:   MatEntry       *Jentry=fd->matentry;

     34:   /* loop over colors  */
     35:   nz = 0;
     36:   for (i=0; i<fd->ncolors; i++) {
     37:     for (j=0; j<fd->nrows[i]; j++) {
     38:       row = Jentry[nz].row;
     39:       y   = fd->M - row - fd->rstart;
     40:       x   = (PetscReal)Jentry[nz++].col;
     41:       PetscDrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1);
     42:     }
     43:   }
     44:   return(0);
     45: }

     47: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd,PetscViewer viewer)
     48: {
     50:   PetscBool      isnull;
     51:   PetscDraw      draw;
     52:   PetscReal      xr,yr,xl,yl,h,w;

     55:   PetscViewerDrawGetDraw(viewer,0,&draw);
     56:   PetscDrawIsNull(draw,&isnull);
     57:   if (isnull) return(0);

     59:   xr   = fd->N; yr  = fd->M; h = yr/10.0; w = xr/10.0;
     60:   xr  += w;     yr += h;    xl = -w;     yl = -h;
     61:   PetscDrawSetCoordinates(draw,xl,yl,xr,yr);
     62:   PetscObjectCompose((PetscObject)fd,"Zoomviewer",(PetscObject)viewer);
     63:   PetscDrawZoom(draw,MatFDColoringView_Draw_Zoom,fd);
     64:   PetscObjectCompose((PetscObject)fd,"Zoomviewer",NULL);
     65:   PetscDrawSave(draw);
     66:   return(0);
     67: }

     69: /*@C
     70:    MatFDColoringView - Views a finite difference coloring context.

     72:    Collective on MatFDColoring

     74:    Input  Parameters:
     75: +  c - the coloring context
     76: -  viewer - visualization context

     78:    Level: intermediate

     80:    Notes:
     81:    The available visualization contexts include
     82: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
     83: .     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
     84:         output where only the first processor opens
     85:         the file.  All other processors send their
     86:         data to the first processor to print.
     87: -     PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure

     89:    Notes:
     90:      Since PETSc uses only a small number of basic colors (currently 33), if the coloring
     91:    involves more than 33 then some seemingly identical colors are displayed making it look
     92:    like an illegal coloring. This is just a graphical artifact.

     94: .seealso: MatFDColoringCreate()

     96: @*/
     97: PetscErrorCode  MatFDColoringView(MatFDColoring c,PetscViewer viewer)
     98: {
     99:   PetscErrorCode    ierr;
    100:   PetscInt          i,j;
    101:   PetscBool         isdraw,iascii;
    102:   PetscViewerFormat format;

    106:   if (!viewer) {
    107:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c),&viewer);
    108:   }

    112:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
    113:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    114:   if (isdraw) {
    115:     MatFDColoringView_Draw(c,viewer);
    116:   } else if (iascii) {
    117:     PetscObjectPrintClassNamePrefixType((PetscObject)c,viewer);
    118:     PetscViewerASCIIPrintf(viewer,"  Error tolerance=%g\n",(double)c->error_rel);
    119:     PetscViewerASCIIPrintf(viewer,"  Umin=%g\n",(double)c->umin);
    120:     PetscViewerASCIIPrintf(viewer,"  Number of colors=%D\n",c->ncolors);

    122:     PetscViewerGetFormat(viewer,&format);
    123:     if (format != PETSC_VIEWER_ASCII_INFO) {
    124:       PetscInt row,col,nz;
    125:       nz = 0;
    126:       for (i=0; i<c->ncolors; i++) {
    127:         PetscViewerASCIIPrintf(viewer,"  Information for color %D\n",i);
    128:         PetscViewerASCIIPrintf(viewer,"    Number of columns %D\n",c->ncolumns[i]);
    129:         for (j=0; j<c->ncolumns[i]; j++) {
    130:           PetscViewerASCIIPrintf(viewer,"      %D\n",c->columns[i][j]);
    131:         }
    132:         PetscViewerASCIIPrintf(viewer,"    Number of rows %D\n",c->nrows[i]);
    133:         if (c->matentry) {
    134:           for (j=0; j<c->nrows[i]; j++) {
    135:             row  = c->matentry[nz].row;
    136:             col  = c->matentry[nz++].col;
    137:             PetscViewerASCIIPrintf(viewer,"      %D %D \n",row,col);
    138:           }
    139:         }
    140:       }
    141:     }
    142:     PetscViewerFlush(viewer);
    143:   }
    144:   return(0);
    145: }

    147: /*@
    148:    MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
    149:    a Jacobian matrix using finite differences.

    151:    Logically Collective on MatFDColoring

    153:    The Jacobian is estimated with the differencing approximation
    154: .vb
    155:        F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
    156:        htype = 'ds':
    157:          h = error_rel*u[i]                 if  abs(u[i]) > umin
    158:            = +/- error_rel*umin             otherwise, with +/- determined by the sign of u[i]
    159:          dx_{i} = (0, ... 1, .... 0)

    161:        htype = 'wp':
    162:          h = error_rel * sqrt(1 + ||u||)
    163: .ve

    165:    Input Parameters:
    166: +  coloring - the coloring context
    167: .  error_rel - relative error
    168: -  umin - minimum allowable u-value magnitude

    170:    Level: advanced

    172: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()

    174: @*/
    175: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd,PetscReal error,PetscReal umin)
    176: {
    181:   if (error != PETSC_DEFAULT) matfd->error_rel = error;
    182:   if (umin != PETSC_DEFAULT)  matfd->umin      = umin;
    183:   return(0);
    184: }

    186: /*@
    187:    MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.

    189:    Logically Collective on MatFDColoring

    191:    Input Parameters:
    192: +  coloring - the coloring context
    193: .  brows - number of rows in the block
    194: -  bcols - number of columns in the block

    196:    Level: intermediate

    198: .seealso: MatFDColoringCreate(), MatFDColoringSetFromOptions()

    200: @*/
    201: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd,PetscInt brows,PetscInt bcols)
    202: {
    207:   if (brows != PETSC_DEFAULT) matfd->brows = brows;
    208:   if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
    209:   return(0);
    210: }

    212: /*@
    213:    MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.

    215:    Collective on Mat

    217:    Input Parameters:
    218: +  mat - the matrix containing the nonzero structure of the Jacobian
    219: .  iscoloring - the coloring of the matrix; usually obtained with MatGetColoring() or DMCreateColoring()
    220: -  color - the matrix coloring context

    222:    Level: beginner

    224:    Notes: When the coloring type is IS_COLORING_LOCAL the coloring is in the local ordering of the unknowns.

    226: .seealso: MatFDColoringCreate(), MatFDColoringDestroy()
    227: @*/
    228: PetscErrorCode MatFDColoringSetUp(Mat mat,ISColoring iscoloring,MatFDColoring color)
    229: {
    231:   PetscBool      eq;

    236:   if (color->setupcalled) return(0);
    237:   PetscObjectCompareId((PetscObject)mat,color->matid,&eq);
    238:   if (!eq) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONG,"Matrix used with MatFDColoringSetUp() must be that used with MatFDColoringCreate()");

    240:   PetscLogEventBegin(MAT_FDColoringSetUp,mat,0,0,0);
    241:   if (mat->ops->fdcoloringsetup) {
    242:     (*mat->ops->fdcoloringsetup)(mat,iscoloring,color);
    243:   } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);

    245:   color->setupcalled = PETSC_TRUE;
    246:    PetscLogEventEnd(MAT_FDColoringSetUp,mat,0,0,0);
    247:   return(0);
    248: }

    250: /*@C
    251:    MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.

    253:    Not Collective

    255:    Input Parameters:
    256: .  coloring - the coloring context

    258:    Output Parameters:
    259: +  f - the function
    260: -  fctx - the optional user-defined function context

    262:    Level: intermediate

    264: .seealso: MatFDColoringCreate(), MatFDColoringSetFunction(), MatFDColoringSetFromOptions()

    266: @*/
    267: PetscErrorCode  MatFDColoringGetFunction(MatFDColoring matfd,PetscErrorCode (**f)(void),void **fctx)
    268: {
    271:   if (f) *f = matfd->f;
    272:   if (fctx) *fctx = matfd->fctx;
    273:   return(0);
    274: }

    276: /*@C
    277:    MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.

    279:    Logically Collective on MatFDColoring

    281:    Input Parameters:
    282: +  coloring - the coloring context
    283: .  f - the function
    284: -  fctx - the optional user-defined function context

    286:    Calling sequence of (*f) function:
    287:     For SNES:    PetscErrorCode (*f)(SNES,Vec,Vec,void*)
    288:     If not using SNES: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored

    290:    Level: advanced

    292:    Notes:
    293:     This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument
    294:      SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by
    295:      calling MatFDColoringApply()

    297:    Fortran Notes:
    298:     In Fortran you must call MatFDColoringSetFunction() for a coloring object to
    299:   be used without SNES or within the SNES solvers.

    301: .seealso: MatFDColoringCreate(), MatFDColoringGetFunction(), MatFDColoringSetFromOptions()

    303: @*/
    304: PetscErrorCode  MatFDColoringSetFunction(MatFDColoring matfd,PetscErrorCode (*f)(void),void *fctx)
    305: {
    308:   matfd->f    = f;
    309:   matfd->fctx = fctx;
    310:   return(0);
    311: }

    313: /*@
    314:    MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
    315:    the options database.

    317:    Collective on MatFDColoring

    319:    The Jacobian, F'(u), is estimated with the differencing approximation
    320: .vb
    321:        F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
    322:        h = error_rel*u[i]                 if  abs(u[i]) > umin
    323:          = +/- error_rel*umin             otherwise, with +/- determined by the sign of u[i]
    324:        dx_{i} = (0, ... 1, .... 0)
    325: .ve

    327:    Input Parameter:
    328: .  coloring - the coloring context

    330:    Options Database Keys:
    331: +  -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
    332: .  -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
    333: .  -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
    334: .  -mat_fd_coloring_view - Activates basic viewing
    335: .  -mat_fd_coloring_view ::ascii_info - Activates viewing info
    336: -  -mat_fd_coloring_view draw - Activates drawing

    338:     Level: intermediate

    340: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()

    342: @*/
    343: PetscErrorCode  MatFDColoringSetFromOptions(MatFDColoring matfd)
    344: {
    346:   PetscBool      flg;
    347:   char           value[3];

    352:   PetscObjectOptionsBegin((PetscObject)matfd);
    353:   PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,NULL);
    354:   PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,NULL);
    355:   PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,sizeof(value),&flg);
    356:   if (flg) {
    357:     if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
    358:     else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
    359:     else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
    360:   }
    361:   PetscOptionsInt("-mat_fd_coloring_brows","Number of block rows","MatFDColoringSetBlockSize",matfd->brows,&matfd->brows,NULL);
    362:   PetscOptionsInt("-mat_fd_coloring_bcols","Number of block columns","MatFDColoringSetBlockSize",matfd->bcols,&matfd->bcols,&flg);
    363:   if (flg && matfd->bcols > matfd->ncolors) {
    364:     /* input bcols cannot be > matfd->ncolors, thus set it as ncolors */
    365:     matfd->bcols = matfd->ncolors;
    366:   }

    368:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
    369:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)matfd);
    370:   PetscOptionsEnd();
    371:   return(0);
    372: }

    374: /*@C
    375:    MatFDColoringSetType - Sets the approach for computing the finite difference parameter

    377:    Collective on MatFDColoring

    379:    Input Parameters:
    380: +  coloring - the coloring context
    381: -  type - either MATMFFD_WP or MATMFFD_DS

    383:    Options Database Keys:
    384: .  -mat_fd_type - "wp" or "ds"

    386:    Note: It is goofy that the argument type is MatMFFDType since the MatFDColoring actually computes the matrix entries
    387:          but the process of computing the entries is the same as as with the MatMFFD operation so we should reuse the names instead of
    388:          introducing another one.

    390:    Level: intermediate

    392: .seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()

    394: @*/
    395: PetscErrorCode  MatFDColoringSetType(MatFDColoring matfd,MatMFFDType type)
    396: {
    399:   /*
    400:      It is goofy to handle the strings this way but currently there is no code to free a dynamically created matfd->htype
    401:      and this function is being provided as patch for a release so it shouldn't change the implementaton
    402:   */
    403:   if (type[0] == 'w' && type[1] == 'p') matfd->htype = "wp";
    404:   else if (type[0] == 'd' && type[1] == 's') matfd->htype = "ds";
    405:   else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",type);
    406:   return(0);
    407: }

    409: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd,const char prefix[],const char optionname[])
    410: {
    411:   PetscErrorCode    ierr;
    412:   PetscBool         flg;
    413:   PetscViewer       viewer;
    414:   PetscViewerFormat format;

    417:   if (prefix) {
    418:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,prefix,optionname,&viewer,&format,&flg);
    419:   } else {
    420:     PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd),((PetscObject)fd)->options,((PetscObject)fd)->prefix,optionname,&viewer,&format,&flg);
    421:   }
    422:   if (flg) {
    423:     PetscViewerPushFormat(viewer,format);
    424:     MatFDColoringView(fd,viewer);
    425:     PetscViewerPopFormat(viewer);
    426:     PetscViewerDestroy(&viewer);
    427:   }
    428:   return(0);
    429: }

    431: /*@
    432:    MatFDColoringCreate - Creates a matrix coloring context for finite difference
    433:    computation of Jacobians.

    435:    Collective on Mat

    437:    Input Parameters:
    438: +  mat - the matrix containing the nonzero structure of the Jacobian
    439: -  iscoloring - the coloring of the matrix; usually obtained with MatColoringCreate() or DMCreateColoring()

    441:     Output Parameter:
    442: .   color - the new coloring context

    444:     Level: intermediate

    446: .seealso: MatFDColoringDestroy(),SNESComputeJacobianDefaultColor(), ISColoringCreate(),
    447:           MatFDColoringSetFunction(), MatFDColoringSetFromOptions(), MatFDColoringApply(),
    448:           MatFDColoringView(), MatFDColoringSetParameters(), MatColoringCreate(), DMCreateColoring(), MatFDColoringSetValues()
    449: @*/
    450: PetscErrorCode  MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color)
    451: {
    452:   MatFDColoring  c;
    453:   MPI_Comm       comm;
    455:   PetscInt       M,N;

    459:   if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();");
    460:   PetscLogEventBegin(MAT_FDColoringCreate,mat,0,0,0);
    461:   MatGetSize(mat,&M,&N);
    462:   if (M != N) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Only for square matrices");
    463:   PetscObjectGetComm((PetscObject)mat,&comm);
    464:   PetscHeaderCreate(c,MAT_FDCOLORING_CLASSID,"MatFDColoring","Jacobian computation via finite differences with coloring","Mat",comm,MatFDColoringDestroy,MatFDColoringView);

    466:   c->ctype = iscoloring->ctype;
    467:   PetscObjectGetId((PetscObject)mat,&c->matid);

    469:   if (mat->ops->fdcoloringcreate) {
    470:     (*mat->ops->fdcoloringcreate)(mat,iscoloring,c);
    471:   } else SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Code not yet written for matrix type %s",((PetscObject)mat)->type_name);

    473:   MatCreateVecs(mat,NULL,&c->w1);
    474:   /* Vec is used intensively in particular piece of scalar CPU code; won't benefit from bouncing back and forth to the GPU */
    475:   VecBindToCPU(c->w1,PETSC_TRUE);
    476:   PetscLogObjectParent((PetscObject)c,(PetscObject)c->w1);
    477:   VecDuplicate(c->w1,&c->w2);
    478:   /* Vec is used intensively in particular piece of scalar CPU code; won't benefit from bouncing back and forth to the GPU */
    479:   VecBindToCPU(c->w2,PETSC_TRUE);
    480:   PetscLogObjectParent((PetscObject)c,(PetscObject)c->w2);

    482:   c->error_rel    = PETSC_SQRT_MACHINE_EPSILON;
    483:   c->umin         = 100.0*PETSC_SQRT_MACHINE_EPSILON;
    484:   c->currentcolor = -1;
    485:   c->htype        = "wp";
    486:   c->fset         = PETSC_FALSE;
    487:   c->setupcalled  = PETSC_FALSE;

    489:   *color = c;
    490:   PetscObjectCompose((PetscObject)mat,"SNESMatFDColoring",(PetscObject)c);
    491:   PetscLogEventEnd(MAT_FDColoringCreate,mat,0,0,0);
    492:   return(0);
    493: }

    495: /*@
    496:     MatFDColoringDestroy - Destroys a matrix coloring context that was created
    497:     via MatFDColoringCreate().

    499:     Collective on MatFDColoring

    501:     Input Parameter:
    502: .   c - coloring context

    504:     Level: intermediate

    506: .seealso: MatFDColoringCreate()
    507: @*/
    508: PetscErrorCode  MatFDColoringDestroy(MatFDColoring *c)
    509: {
    511:   PetscInt       i;
    512:   MatFDColoring  color = *c;

    515:   if (!*c) return(0);
    516:   if (--((PetscObject)color)->refct > 0) {*c = NULL; return(0);}

    518:   /* we do not free the column arrays since their entries are owned by the ISs in color->isa */
    519:   for (i=0; i<color->ncolors; i++) {
    520:     ISDestroy(&color->isa[i]);
    521:   }
    522:   PetscFree(color->isa);
    523:   PetscFree2(color->ncolumns,color->columns);
    524:   PetscFree(color->nrows);
    525:   if (color->htype[0] == 'w') {
    526:     PetscFree(color->matentry2);
    527:   } else {
    528:     PetscFree(color->matentry);
    529:   }
    530:   PetscFree(color->dy);
    531:   if (color->vscale) {VecDestroy(&color->vscale);}
    532:   VecDestroy(&color->w1);
    533:   VecDestroy(&color->w2);
    534:   VecDestroy(&color->w3);
    535:   PetscHeaderDestroy(c);
    536:   return(0);
    537: }

    539: /*@C
    540:     MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
    541:       that are currently being perturbed.

    543:     Not Collective

    545:     Input Parameters:
    546: .   coloring - coloring context created with MatFDColoringCreate()

    548:     Output Parameters:
    549: +   n - the number of local columns being perturbed
    550: -   cols - the column indices, in global numbering

    552:    Level: advanced

    554:    Note: IF the matrix type is BAIJ, then the block column indices are returned

    556:    Fortran Note:
    557:    This routine has a different interface for Fortran
    558: #include <petsc/finclude/petscmat.h>
    559: $          use petscmat
    560: $          PetscInt, pointer :: array(:)
    561: $          PetscErrorCode  ierr
    562: $          MatFDColoring   i
    563: $          call MatFDColoringGetPerturbedColumnsF90(i,array,ierr)
    564: $      use the entries of array ...
    565: $          call MatFDColoringRestorePerturbedColumnsF90(i,array,ierr)

    567: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringApply()

    569: @*/
    570: PetscErrorCode  MatFDColoringGetPerturbedColumns(MatFDColoring coloring,PetscInt *n,const PetscInt *cols[])
    571: {
    573:   if (coloring->currentcolor >= 0) {
    574:     *n    = coloring->ncolumns[coloring->currentcolor];
    575:     *cols = coloring->columns[coloring->currentcolor];
    576:   } else {
    577:     *n = 0;
    578:   }
    579:   return(0);
    580: }

    582: /*@
    583:     MatFDColoringApply - Given a matrix for which a MatFDColoring context
    584:     has been created, computes the Jacobian for a function via finite differences.

    586:     Collective on MatFDColoring

    588:     Input Parameters:
    589: +   mat - location to store Jacobian
    590: .   coloring - coloring context created with MatFDColoringCreate()
    591: .   x1 - location at which Jacobian is to be computed
    592: -   sctx - context required by function, if this is being used with the SNES solver then it is SNES object, otherwise it is null

    594:     Options Database Keys:
    595: +    -mat_fd_type - "wp" or "ds"  (see MATMFFD_WP or MATMFFD_DS)
    596: .    -mat_fd_coloring_view - Activates basic viewing or coloring
    597: .    -mat_fd_coloring_view draw - Activates drawing of coloring
    598: -    -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info

    600:     Level: intermediate

    602: .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView(), MatFDColoringSetFunction(), MatFDColoringSetValues()

    604: @*/
    605: PetscErrorCode  MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
    606: {
    608:   PetscBool      eq;

    614:   PetscObjectCompareId((PetscObject)J,coloring->matid,&eq);
    615:   if (!eq) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Matrix used with MatFDColoringApply() must be that used with MatFDColoringCreate()");
    616:   if (!coloring->f) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetFunction()");
    617:   if (!J->ops->fdcoloringapply) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported for this matrix type %s",((PetscObject)J)->type_name);
    618:   if (!coloring->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatFDColoringSetUp()");

    620:   MatSetUnfactored(J);
    621:   PetscLogEventBegin(MAT_FDColoringApply,coloring,J,x1,0);
    622:   (*J->ops->fdcoloringapply)(J,coloring,x1,sctx);
    623:   PetscLogEventEnd(MAT_FDColoringApply,coloring,J,x1,0);
    624:   if (!coloring->viewed) {
    625:     MatFDColoringViewFromOptions(coloring,NULL,"-mat_fd_coloring_view");
    626:     coloring->viewed = PETSC_TRUE;
    627:   }
    628:   return(0);
    629: }