Actual source code: fdmatrix.c

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

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

  9: PetscErrorCode MatFDColoringSetF(MatFDColoring fd, Vec F)
 10: {
 11:   PetscFunctionBegin;
 12:   if (F) {
 13:     PetscCall(VecCopy(F, fd->w1));
 14:     fd->fset = PETSC_TRUE;
 15:   } else {
 16:     fd->fset = PETSC_FALSE;
 17:   }
 18:   PetscFunctionReturn(PETSC_SUCCESS);
 19: }

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

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

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

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

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

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

 73:   Collective

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

 79:   Level: intermediate

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

 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: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`
 95: @*/
 96: PetscErrorCode MatFDColoringView(MatFDColoring c, PetscViewer viewer)
 97: {
 98:   PetscInt          i, j;
 99:   PetscBool         isdraw, iascii;
100:   PetscViewerFormat format;

102:   PetscFunctionBegin;
104:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c), &viewer));
106:   PetscCheckSameComm(c, 1, viewer, 2);

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

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

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

145:   Logically Collective

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

152:   Level: advanced

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

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

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

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

183:   Logically Collective

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

190:   Level: intermediate

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

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

208:   Collective

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

215:   Level: beginner

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

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

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

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

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

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

244:   Not Collective

246:   Input Parameter:
247: . matfd - the coloring context

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

253:   Level: intermediate

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

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

269:   Logically Collective

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

276:   Level: advanced

278:   Note:
279:   `f` has two possible calling configurations\:
280: $ PetscErrorCode f(SNES snes, Vec in, Vec out, void *fctx)
281: + snes - the nonlinear solver `SNES` object
282: . in   - the location where the Jacobian is to be computed
283: . out  - the location to put the computed function value
284: - fctx - the function context

286:   and
287: $ PetscErrorCode f(void *dummy, Vec in, Vec out, void *fctx)
288: + dummy - an unused parameter
289: . in    - the location where the Jacobian is to be computed
290: . out   - the location to put the computed function value
291: - fctx  - the function context

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: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringGetFunction()`, `MatFDColoringSetFromOptions()`
302: @*/
303: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd, PetscErrorCode (*f)(void), void *fctx)
304: {
305:   PetscFunctionBegin;
307:   matfd->f    = f;
308:   matfd->fctx = fctx;
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

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

316:   Collective

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

326:   Input Parameter:
327: . matfd - the coloring context

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

337:   Level: intermediate

339: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringView()`, `MatFDColoringSetParameters()`
340: @*/
341: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
342: {
343:   PetscBool flg;
344:   char      value[3];

346:   PetscFunctionBegin;

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

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

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

374:   Collective

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

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

383:   Level: intermediate

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

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

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

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

431:   Collective

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

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

440:   Level: intermediate

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

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

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

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

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

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

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

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

490:   Collective

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

495:   Level: intermediate

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

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

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

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

534:   Not Collective

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

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

543:   Level: advanced

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

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

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

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

579:   Collective

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

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

593:   Level: intermediate

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

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

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