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: {
 12:   PetscFunctionBegin;
 13:   if (F) {
 14:     PetscCall(VecCopy(F, fd->w1));
 15:     fd->fset = PETSC_TRUE;
 16:   } else {
 17:     fd->fset = PETSC_FALSE;
 18:   }
 19:   PetscFunctionReturn(PETSC_SUCCESS);
 20: }

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

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

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

 50:   PetscFunctionBegin;
 51:   PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
 52:   PetscCall(PetscDrawIsNull(draw, &isnull));
 53:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

 55:   xr = fd->N;
 56:   yr = fd->M;
 57:   h  = yr / 10.0;
 58:   w  = xr / 10.0;
 59:   xr += w;
 60:   yr += h;
 61:   xl = -w;
 62:   yl = -h;
 63:   PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
 64:   PetscCall(PetscObjectCompose((PetscObject)fd, "Zoomviewer", (PetscObject)viewer));
 65:   PetscCall(PetscDrawZoom(draw, MatFDColoringView_Draw_Zoom, fd));
 66:   PetscCall(PetscObjectCompose((PetscObject)fd, "Zoomviewer", NULL));
 67:   PetscCall(PetscDrawSave(draw));
 68:   PetscFunctionReturn(PETSC_SUCCESS);
 69: }

 71: /*@C
 72:    MatFDColoringView - Views a finite difference coloring context.

 74:    Collective

 76:    Input Parameters:
 77: +  c - the coloring context
 78: -  viewer - visualization context

 80:    Level: intermediate

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

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

 95: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`
 96: @*/
 97: PetscErrorCode MatFDColoringView(MatFDColoring c, PetscViewer viewer)
 98: {
 99:   PetscInt          i, j;
100:   PetscBool         isdraw, iascii;
101:   PetscViewerFormat format;

103:   PetscFunctionBegin;
105:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c), &viewer));
107:   PetscCheckSameComm(c, 1, viewer, 2);

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

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

142: /*@
143:    MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
144:    a Jacobian matrix using finite differences.

146:    Logically Collective

148:    Input Parameters:
149: +  matfd - the coloring context
150: .  error - relative error
151: -  umin - minimum allowable u-value magnitude

153:    Level: advanced

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

164:        htype = 'wp':
165:          h = error_rel * sqrt(1 + ||u||)
166: .ve

168: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringSetFromOptions()`
169: @*/
170: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd, PetscReal error, PetscReal umin)
171: {
172:   PetscFunctionBegin;
176:   if (error != (PetscReal)PETSC_DEFAULT) matfd->error_rel = error;
177:   if (umin != (PetscReal)PETSC_DEFAULT) matfd->umin = umin;
178:   PetscFunctionReturn(PETSC_SUCCESS);
179: }

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

184:    Logically Collective

186:    Input Parameters:
187: +  coloring - the coloring context
188: .  brows - number of rows in the block
189: -  bcols - number of columns in the block

191:    Level: intermediate

193: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringSetFromOptions()`
194: @*/
195: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd, PetscInt brows, PetscInt bcols)
196: {
197:   PetscFunctionBegin;
201:   if (brows != PETSC_DEFAULT) matfd->brows = brows;
202:   if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
203:   PetscFunctionReturn(PETSC_SUCCESS);
204: }

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

209:    Collective

211:    Input Parameters:
212: +  mat - the matrix containing the nonzero structure of the Jacobian
213: .  iscoloring - the coloring of the matrix; usually obtained with `MatGetColoring()` or `DMCreateColoring()`
214: -  color - the matrix coloring context

216:    Level: beginner

218:    Notes:
219:    When the coloring type is `IS_COLORING_LOCAL` the coloring is in the local ordering of the unknowns.

221: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringDestroy()`
222: @*/
223: PetscErrorCode MatFDColoringSetUp(Mat mat, ISColoring iscoloring, MatFDColoring color)
224: {
225:   PetscBool eq;

227:   PetscFunctionBegin;
230:   if (color->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
231:   PetscCall(PetscObjectCompareId((PetscObject)mat, color->matid, &eq));
232:   PetscCheck(eq, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Matrix used with MatFDColoringSetUp() must be that used with MatFDColoringCreate()");

234:   PetscCall(PetscLogEventBegin(MAT_FDColoringSetUp, mat, 0, 0, 0));
235:   PetscUseTypeMethod(mat, fdcoloringsetup, iscoloring, color);

237:   color->setupcalled = PETSC_TRUE;
238:   PetscCall(PetscLogEventEnd(MAT_FDColoringSetUp, mat, 0, 0, 0));
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

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

245:    Not Collective

247:    Input Parameter:
248: .  coloring - the coloring context

250:    Output Parameters:
251: +  f - the function
252: -  fctx - the optional user-defined function context

254:    Level: intermediate

256: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`, `MatFDColoringSetFromOptions()`
257: @*/
258: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd, PetscErrorCode (**f)(void), void **fctx)
259: {
260:   PetscFunctionBegin;
262:   if (f) *f = matfd->f;
263:   if (fctx) *fctx = matfd->fctx;
264:   PetscFunctionReturn(PETSC_SUCCESS);
265: }

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

270:    Logically Collective

272:    Input Parameters:
273: +  coloring - the coloring context
274: .  f - the function
275: -  fctx - the optional user-defined function context

277:    Calling sequence with `SNES` of `f`:
278: $   PetscErrorCode f(SNES snes, Vec in, Vec out, void *fctx)
279: +  snes - the nonlinear solver `SNES` object
280: .  in - the location where the Jacobian is to be computed
281: .  out - the location to put the computed function value
282: -  fctx - the function context

284:    Calling sequence without `SNES` of `f`:
285: $   PetscErrorCode f(void *dummy, Vec in, Vec out, void *fctx)
286: +  dummy - an unused parameter
287: .  in - the location where the Jacobian is to be computed
288: .  out - the location to put the computed function value
289: -  fctx - the function context

291:    Level: advanced

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

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

302: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringGetFunction()`, `MatFDColoringSetFromOptions()`
303: @*/
304: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd, PetscErrorCode (*f)(void), void *fctx)
305: {
306:   PetscFunctionBegin;
308:   matfd->f    = f;
309:   matfd->fctx = fctx;
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

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

317:    Collective

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: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringView()`, `MatFDColoringSetParameters()`
341: @*/
342: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
343: {
344:   PetscBool flg;
345:   char      value[3];

347:   PetscFunctionBegin;

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

366:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
367:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)matfd, PetscOptionsObject));
368:   PetscOptionsEnd();
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }

372: /*@C
373:    MatFDColoringSetType - Sets the approach for computing the finite difference parameter

375:    Collective

377:    Input Parameters:
378: +  coloring - the coloring context
379: -  type - either `MATMFFD_WP` or `MATMFFD_DS`

381:    Options Database Key:
382: .  -mat_fd_type - "wp" or "ds"

384:    Level: intermediate

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

391: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringView()`, `MatFDColoringSetParameters()`
392: @*/
393: PetscErrorCode MatFDColoringSetType(MatFDColoring matfd, MatMFFDType type)
394: {
395:   PetscFunctionBegin;
397:   /*
398:      It is goofy to handle the strings this way but currently there is no code to free a dynamically created matfd->htype
399:      and this function is being provided as patch for a release so it shouldn't change the implementation
400:   */
401:   if (type[0] == 'w' && type[1] == 'p') matfd->htype = "wp";
402:   else if (type[0] == 'd' && type[1] == 's') matfd->htype = "ds";
403:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown finite differencing type %s", type);
404:   PetscFunctionReturn(PETSC_SUCCESS);
405: }

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

413:   PetscFunctionBegin;
414:   if (prefix) {
415:     PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd), ((PetscObject)fd)->options, prefix, optionname, &viewer, &format, &flg));
416:   } else {
417:     PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd), ((PetscObject)fd)->options, ((PetscObject)fd)->prefix, optionname, &viewer, &format, &flg));
418:   }
419:   if (flg) {
420:     PetscCall(PetscViewerPushFormat(viewer, format));
421:     PetscCall(MatFDColoringView(fd, viewer));
422:     PetscCall(PetscViewerPopFormat(viewer));
423:     PetscCall(PetscViewerDestroy(&viewer));
424:   }
425:   PetscFunctionReturn(PETSC_SUCCESS);
426: }

428: /*@
429:    MatFDColoringCreate - Creates a matrix coloring context for finite difference
430:    computation of Jacobians.

432:    Collective

434:    Input Parameters:
435: +  mat - the matrix containing the nonzero structure of the Jacobian
436: -  iscoloring - the coloring of the matrix; usually obtained with `MatColoringCreate()` or `DMCreateColoring()`

438:     Output Parameter:
439: .   color - the new coloring context

441:     Level: intermediate

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

453:   PetscFunctionBegin;
455:   PetscCheck(mat->assembled, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Matrix must be assembled by calls to MatAssemblyBegin/End();");
456:   PetscCall(PetscLogEventBegin(MAT_FDColoringCreate, mat, 0, 0, 0));
457:   PetscCall(MatGetSize(mat, &M, &N));
458:   PetscCheck(M == N, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Only for square matrices");
459:   PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
460:   PetscCall(PetscHeaderCreate(c, MAT_FDCOLORING_CLASSID, "MatFDColoring", "Jacobian computation via finite differences with coloring", "Mat", comm, MatFDColoringDestroy, MatFDColoringView));

462:   c->ctype = iscoloring->ctype;
463:   PetscCall(PetscObjectGetId((PetscObject)mat, &c->matid));

465:   PetscUseTypeMethod(mat, fdcoloringcreate, iscoloring, c);

467:   PetscCall(MatCreateVecs(mat, NULL, &c->w1));
468:   /* Vec is used intensively in particular piece of scalar CPU code; won't benefit from bouncing back and forth to the GPU */
469:   PetscCall(VecBindToCPU(c->w1, PETSC_TRUE));
470:   PetscCall(VecDuplicate(c->w1, &c->w2));
471:   /* Vec is used intensively in particular piece of scalar CPU code; won't benefit from bouncing back and forth to the GPU */
472:   PetscCall(VecBindToCPU(c->w2, PETSC_TRUE));

474:   c->error_rel    = PETSC_SQRT_MACHINE_EPSILON;
475:   c->umin         = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
476:   c->currentcolor = -1;
477:   c->htype        = "wp";
478:   c->fset         = PETSC_FALSE;
479:   c->setupcalled  = PETSC_FALSE;

481:   *color = c;
482:   PetscCall(PetscObjectCompose((PetscObject)mat, "SNESMatFDColoring", (PetscObject)c));
483:   PetscCall(PetscLogEventEnd(MAT_FDColoringCreate, mat, 0, 0, 0));
484:   PetscFunctionReturn(PETSC_SUCCESS);
485: }

487: /*@
488:     MatFDColoringDestroy - Destroys a matrix coloring context that was created
489:     via `MatFDColoringCreate()`.

491:     Collective

493:     Input Parameter:
494: .   c - coloring context

496:     Level: intermediate

498: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`
499: @*/
500: PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
501: {
502:   PetscInt      i;
503:   MatFDColoring color = *c;

505:   PetscFunctionBegin;
506:   if (!*c) PetscFunctionReturn(PETSC_SUCCESS);
507:   if (--((PetscObject)color)->refct > 0) {
508:     *c = NULL;
509:     PetscFunctionReturn(PETSC_SUCCESS);
510:   }

512:   /* we do not free the column arrays since their entries are owned by the ISs in color->isa */
513:   for (i = 0; i < color->ncolors; i++) PetscCall(ISDestroy(&color->isa[i]));
514:   PetscCall(PetscFree(color->isa));
515:   PetscCall(PetscFree2(color->ncolumns, color->columns));
516:   PetscCall(PetscFree(color->nrows));
517:   if (color->htype[0] == 'w') {
518:     PetscCall(PetscFree(color->matentry2));
519:   } else {
520:     PetscCall(PetscFree(color->matentry));
521:   }
522:   PetscCall(PetscFree(color->dy));
523:   if (color->vscale) PetscCall(VecDestroy(&color->vscale));
524:   PetscCall(VecDestroy(&color->w1));
525:   PetscCall(VecDestroy(&color->w2));
526:   PetscCall(VecDestroy(&color->w3));
527:   PetscCall(PetscHeaderDestroy(c));
528:   PetscFunctionReturn(PETSC_SUCCESS);
529: }

531: /*@C
532:     MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
533:       that are currently being perturbed.

535:     Not Collective

537:     Input Parameter:
538: .   coloring - coloring context created with `MatFDColoringCreate()`

540:     Output Parameters:
541: +   n - the number of local columns being perturbed
542: -   cols - the column indices, in global numbering

544:    Level: advanced

546:    Note:
547:    IF the matrix type is `MATBAIJ`, then the block column indices are returned

549:    Fortran Note:
550:    This routine has a different interface for Fortran
551: .vb
552: #include <petsc/finclude/petscmat.h>
553:           use petscmat
554:           PetscInt, pointer :: array(:)
555:           PetscErrorCode  ierr
556:           MatFDColoring   i
557:           call MatFDColoringGetPerturbedColumnsF90(i,array,ierr)
558:       use the entries of array ...
559:           call MatFDColoringRestorePerturbedColumnsF90(i,array,ierr)
560: .ve

562: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringDestroy()`, `MatFDColoringView()`, `MatFDColoringApply()`
563: @*/
564: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring, PetscInt *n, const PetscInt *cols[])
565: {
566:   PetscFunctionBegin;
567:   if (coloring->currentcolor >= 0) {
568:     *n    = coloring->ncolumns[coloring->currentcolor];
569:     *cols = coloring->columns[coloring->currentcolor];
570:   } else {
571:     *n = 0;
572:   }
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }

576: /*@
577:     MatFDColoringApply - Given a matrix for which a `MatFDColoring` context
578:     has been created, computes the Jacobian for a function via finite differences.

580:     Collective

582:     Input Parameters:
583: +   mat - location to store Jacobian
584: .   coloring - coloring context created with `MatFDColoringCreate()`
585: .   x1 - location at which Jacobian is to be computed
586: -   sctx - context required by function, if this is being used with the SNES solver then it is `SNES` object, otherwise it is null

588:     Options Database Keys:
589: +    -mat_fd_type - "wp" or "ds"  (see `MATMFFD_WP` or `MATMFFD_DS`)
590: .    -mat_fd_coloring_view - Activates basic viewing or coloring
591: .    -mat_fd_coloring_view draw - Activates drawing of coloring
592: -    -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info

594:     Level: intermediate

596: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringDestroy()`, `MatFDColoringView()`, `MatFDColoringSetFunction()`, `MatFDColoringSetValues()`
597: @*/
598: PetscErrorCode MatFDColoringApply(Mat J, MatFDColoring coloring, Vec x1, void *sctx)
599: {
600:   PetscBool eq;

602:   PetscFunctionBegin;
606:   PetscCall(PetscObjectCompareId((PetscObject)J, coloring->matid, &eq));
607:   PetscCheck(eq, PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONG, "Matrix used with MatFDColoringApply() must be that used with MatFDColoringCreate()");
608:   PetscCheck(coloring->f, PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONGSTATE, "Must call MatFDColoringSetFunction()");
609:   PetscCheck(coloring->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call MatFDColoringSetUp()");

611:   PetscCall(MatSetUnfactored(J));
612:   PetscCall(PetscLogEventBegin(MAT_FDColoringApply, coloring, J, x1, 0));
613:   PetscUseTypeMethod(J, fdcoloringapply, coloring, x1, sctx);
614:   PetscCall(PetscLogEventEnd(MAT_FDColoringApply, coloring, J, x1, 0));
615:   if (!coloring->viewed) {
616:     PetscCall(MatFDColoringViewFromOptions(coloring, NULL, "-mat_fd_coloring_view"));
617:     coloring->viewed = PETSC_TRUE;
618:   }
619:   PetscFunctionReturn(PETSC_SUCCESS);
620: }