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: }