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: /*@
10: MatFDColoringSetF - Cache the current function value used by the finite-difference coloring context to
11: avoid recomputing `F(x)` during a Jacobian evaluation.
13: Logically Collective
15: Input Parameters:
16: + fd - the `MatFDColoring` context
17: - F - the current function value `F(x)`, or `NULL` to invalidate any cached value
19: Level: advanced
21: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringApply()`
22: @*/
23: PetscErrorCode MatFDColoringSetF(MatFDColoring fd, Vec F)
24: {
25: PetscFunctionBegin;
26: if (F) {
27: PetscCall(VecCopy(F, fd->w1));
28: fd->fset = PETSC_TRUE;
29: } else {
30: fd->fset = PETSC_FALSE;
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: #include <petscdraw.h>
36: static PetscErrorCode MatFDColoringView_Draw_Zoom(PetscDraw draw, void *Aa)
37: {
38: MatFDColoring fd = (MatFDColoring)Aa;
39: PetscMPIInt i, j, nz;
40: PetscInt row;
41: PetscReal x, y;
42: MatEntry *Jentry = fd->matentry;
44: PetscFunctionBegin;
45: /* loop over colors */
46: nz = 0;
47: for (i = 0; i < fd->ncolors; i++) {
48: for (j = 0; j < fd->nrows[i]; j++) {
49: row = Jentry[nz].row;
50: y = fd->M - row - fd->rstart;
51: x = (PetscReal)Jentry[nz++].col;
52: PetscCall(PetscDrawRectangle(draw, x, y, x + 1, y + 1, i + 1, i + 1, i + 1, i + 1));
53: }
54: }
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd, PetscViewer viewer)
59: {
60: PetscBool isnull;
61: PetscDraw draw;
62: PetscReal xr, yr, xl, yl, h, w;
64: PetscFunctionBegin;
65: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
66: PetscCall(PetscDrawIsNull(draw, &isnull));
67: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
69: xr = fd->N;
70: yr = fd->M;
71: h = yr / 10.0;
72: w = xr / 10.0;
73: xr += w;
74: yr += h;
75: xl = -w;
76: yl = -h;
77: PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
78: PetscCall(PetscObjectCompose((PetscObject)fd, "Zoomviewer", (PetscObject)viewer));
79: PetscCall(PetscDrawZoom(draw, MatFDColoringView_Draw_Zoom, fd));
80: PetscCall(PetscObjectCompose((PetscObject)fd, "Zoomviewer", NULL));
81: PetscCall(PetscDrawSave(draw));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@
86: MatFDColoringView - Views a finite difference coloring context.
88: Collective
90: Input Parameters:
91: + c - the coloring context
92: - viewer - visualization context
94: Level: intermediate
96: Notes:
97: The available visualization contexts include
98: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
99: . `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
100: output where only the first processor opens
101: the file. All other processors send their
102: data to the first processor to print.
103: - `PETSC_VIEWER_DRAW_WORLD` - graphical display of nonzero structure
105: Since PETSc uses only a small number of basic colors (currently 33), if the coloring
106: involves more than 33 then some seemingly identical colors are displayed making it look
107: like an illegal coloring. This is just a graphical artifact.
109: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`
110: @*/
111: PetscErrorCode MatFDColoringView(MatFDColoring c, PetscViewer viewer)
112: {
113: PetscInt i, j;
114: PetscBool isdraw, isascii;
115: PetscViewerFormat format;
117: PetscFunctionBegin;
119: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c), &viewer));
121: PetscCheckSameComm(c, 1, viewer, 2);
123: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
124: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
125: if (isdraw) {
126: PetscCall(MatFDColoringView_Draw(c, viewer));
127: } else if (isascii) {
128: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)c, viewer));
129: PetscCall(PetscViewerASCIIPrintf(viewer, " Error tolerance=%g\n", (double)c->error_rel));
130: PetscCall(PetscViewerASCIIPrintf(viewer, " Umin=%g\n", (double)c->umin));
131: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of colors=%" PetscInt_FMT "\n", c->ncolors));
133: PetscCall(PetscViewerGetFormat(viewer, &format));
134: if (format != PETSC_VIEWER_ASCII_INFO) {
135: PetscInt row, col, nz;
136: nz = 0;
137: for (i = 0; i < c->ncolors; i++) {
138: PetscCall(PetscViewerASCIIPrintf(viewer, " Information for color %" PetscInt_FMT "\n", i));
139: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of columns %" PetscInt_FMT "\n", c->ncolumns[i]));
140: for (j = 0; j < c->ncolumns[i]; j++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT "\n", c->columns[i][j]));
141: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of rows %" PetscInt_FMT "\n", c->nrows[i]));
142: if (c->matentry) {
143: for (j = 0; j < c->nrows[i]; j++) {
144: row = c->matentry[nz].row;
145: col = c->matentry[nz++].col;
146: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT " \n", row, col));
147: }
148: }
149: }
150: }
151: PetscCall(PetscViewerFlush(viewer));
152: }
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: /*@
157: MatFDColoringSetParameters - Sets the parameters for the approximation of
158: a sparse Jacobian matrix using finite differences and matrix coloring
160: Logically Collective
162: Input Parameters:
163: + matfd - the coloring context
164: . error - relative error
165: - umin - minimum allowable u-value magnitude
167: Level: advanced
169: Note:
170: The Jacobian is estimated with the differencing approximation
171: .vb
172: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
173: htype = 'ds':
174: h = error_rel*u[i] if abs(u[i]) > umin
175: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
176: dx_{i} = (0, ... 1, .... 0)
178: htype = 'wp':
179: h = error_rel * sqrt(1 + ||u||)
180: .ve
182: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringSetFromOptions()`
183: @*/
184: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd, PetscReal error, PetscReal umin)
185: {
186: PetscFunctionBegin;
190: if (error != (PetscReal)PETSC_DEFAULT) matfd->error_rel = error;
191: if (umin != (PetscReal)PETSC_DEFAULT) matfd->umin = umin;
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: /*@
196: MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.
198: Logically Collective
200: Input Parameters:
201: + matfd - the coloring context
202: . brows - number of rows in the block
203: - bcols - number of columns in the block
205: Level: intermediate
207: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringSetFromOptions()`
208: @*/
209: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd, PetscInt brows, PetscInt bcols)
210: {
211: PetscFunctionBegin;
215: if (brows != PETSC_DEFAULT) matfd->brows = brows;
216: if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: /*@
221: MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.
223: Collective
225: Input Parameters:
226: + mat - the matrix containing the nonzero structure of the Jacobian
227: . iscoloring - the coloring of the matrix; usually obtained with `MatGetColoring()` or `DMCreateColoring()`
228: - color - the matrix coloring context
230: Level: beginner
232: Notes:
233: When the coloring type is `IS_COLORING_LOCAL` the coloring is in the local ordering of the unknowns.
235: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringDestroy()`
236: @*/
237: PetscErrorCode MatFDColoringSetUp(Mat mat, ISColoring iscoloring, MatFDColoring color)
238: {
239: PetscBool eq;
241: PetscFunctionBegin;
244: if (color->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
245: PetscCall(PetscObjectCompareId((PetscObject)mat, color->matid, &eq));
246: PetscCheck(eq, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONG, "Matrix used with MatFDColoringSetUp() must be that used with MatFDColoringCreate()");
248: PetscCall(PetscLogEventBegin(MAT_FDColoringSetUp, mat, 0, 0, 0));
249: PetscUseTypeMethod(mat, fdcoloringsetup, iscoloring, color);
251: color->setupcalled = PETSC_TRUE;
252: PetscCall(PetscLogEventEnd(MAT_FDColoringSetUp, mat, 0, 0, 0));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@C
257: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
259: Not Collective
261: Input Parameter:
262: . matfd - the coloring context
264: Output Parameters:
265: + f - the function, see `MatFDColoringFn` for the calling sequence
266: - fctx - the optional user-defined function context
268: Level: intermediate
270: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`, `MatFDColoringSetFromOptions()`, `MatFDColoringFn`
271: @*/
272: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd, MatFDColoringFn **f, void **fctx)
273: {
274: PetscFunctionBegin;
276: if (f) *f = matfd->f;
277: if (fctx) *fctx = matfd->fctx;
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: /*@C
282: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
284: Logically Collective
286: Input Parameters:
287: + matfd - the coloring context
288: . f - the function, see `MatFDColoringFn` for the calling sequence
289: - fctx - the optional user-defined 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()`, `MatFDColoringFn`
303: @*/
304: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd, MatFDColoringFn *f, 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: . matfd - 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|ds) - See `MATMFFD_WP` and `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: /*@
373: MatFDColoringSetType - Sets the approach for computing the finite difference parameter
375: Collective
377: Input Parameters:
378: + matfd - 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 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: PetscCall(PetscObjectChangeTypeName((PetscObject)matfd, type));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: static PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd, const char prefix[], const char optionname[])
409: {
410: PetscBool flg;
411: PetscViewer viewer;
412: PetscViewerFormat format;
414: PetscFunctionBegin;
415: if (prefix) {
416: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)fd), ((PetscObject)fd)->options, prefix, optionname, &viewer, &format, &flg));
417: } else {
418: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)fd), ((PetscObject)fd)->options, ((PetscObject)fd)->prefix, optionname, &viewer, &format, &flg));
419: }
420: if (flg) {
421: PetscCall(PetscViewerPushFormat(viewer, format));
422: PetscCall(MatFDColoringView(fd, viewer));
423: PetscCall(PetscViewerPopFormat(viewer));
424: PetscCall(PetscViewerDestroy(&viewer));
425: }
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: /*@
430: MatFDColoringCreate - Creates a matrix coloring context for finite difference
431: computation of Jacobians.
433: Collective
435: Input Parameters:
436: + mat - the matrix containing the nonzero structure of the Jacobian
437: - iscoloring - the coloring of the matrix; usually obtained with `MatColoringCreate()` or `DMCreateColoring()`
439: Output Parameter:
440: . color - the new coloring context
442: Level: intermediate
444: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringDestroy()`, `SNESComputeJacobianDefaultColor()`, `ISColoringCreate()`,
445: `MatFDColoringSetFunction()`, `MatFDColoringSetFromOptions()`, `MatFDColoringApply()`,
446: `MatFDColoringView()`, `MatFDColoringSetParameters()`, `MatColoringCreate()`, `DMCreateColoring()`, `MatFDColoringSetValues()`
447: @*/
448: PetscErrorCode MatFDColoringCreate(Mat mat, ISColoring iscoloring, MatFDColoring *color)
449: {
450: MatFDColoring c;
451: MPI_Comm comm;
452: PetscInt M, N;
454: PetscFunctionBegin;
456: PetscAssertPointer(color, 3);
457: PetscCheck(mat->assembled, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Matrix must be assembled by calls to MatAssemblyBegin/End();");
458: PetscCall(PetscLogEventBegin(MAT_FDColoringCreate, mat, 0, 0, 0));
459: PetscCall(MatGetSize(mat, &M, &N));
460: PetscCheck(M == N, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Only for square matrices");
461: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
462: PetscCall(PetscHeaderCreate(c, MAT_FDCOLORING_CLASSID, "MatFDColoring", "Jacobian computation via finite differences with coloring", "Mat", comm, MatFDColoringDestroy, MatFDColoringView));
464: c->ctype = iscoloring->ctype;
465: PetscCall(PetscObjectGetId((PetscObject)mat, &c->matid));
467: PetscUseTypeMethod(mat, fdcoloringcreate, iscoloring, c);
469: PetscCall(MatCreateVecs(mat, NULL, &c->w1));
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->w1, PETSC_TRUE));
472: PetscCall(VecDuplicate(c->w1, &c->w2));
473: /* Vec is used intensively in particular piece of scalar CPU code; won't benefit from bouncing back and forth to the GPU */
474: PetscCall(VecBindToCPU(c->w2, PETSC_TRUE));
476: c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
477: c->umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
478: c->currentcolor = -1;
479: c->htype = "wp";
480: c->fset = PETSC_FALSE;
481: c->setupcalled = PETSC_FALSE;
483: *color = c;
484: PetscCall(PetscObjectCompose((PetscObject)mat, "SNESMatFDColoring", (PetscObject)c));
485: PetscCall(PetscLogEventEnd(MAT_FDColoringCreate, mat, 0, 0, 0));
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: /*@
490: MatFDColoringDestroy - Destroys a matrix coloring context that was created
491: via `MatFDColoringCreate()`.
493: Collective
495: Input Parameter:
496: . c - coloring context
498: Level: intermediate
500: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`
501: @*/
502: PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
503: {
504: PetscInt i;
505: MatFDColoring color = *c;
507: PetscFunctionBegin;
508: if (!*c) PetscFunctionReturn(PETSC_SUCCESS);
509: if (--((PetscObject)color)->refct > 0) {
510: *c = NULL;
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: /* we do not free the column arrays since their entries are owned by the ISs in color->isa */
515: for (i = 0; i < color->ncolors; i++) PetscCall(ISDestroy(&color->isa[i]));
516: PetscCall(PetscFree(color->isa));
517: PetscCall(PetscFree2(color->ncolumns, color->columns));
518: PetscCall(PetscFree(color->nrows));
519: if (color->htype[0] == 'w') {
520: PetscCall(PetscFree(color->matentry2));
521: } else {
522: PetscCall(PetscFree(color->matentry));
523: }
524: PetscCall(PetscFree(color->dy));
525: PetscCall(VecDestroy(&color->vscale));
526: PetscCall(VecDestroy(&color->w1));
527: PetscCall(VecDestroy(&color->w2));
528: PetscCall(VecDestroy(&color->w3));
529: PetscCall(PetscHeaderDestroy(c));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /*@C
534: MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
535: that are currently being perturbed.
537: Not Collective
539: Input Parameter:
540: . coloring - coloring context created with `MatFDColoringCreate()`
542: Output Parameters:
543: + n - the number of local columns being perturbed
544: - cols - the column indices, in global numbering
546: Level: advanced
548: Note:
549: IF the matrix type is `MATBAIJ`, then the block column indices are returned
551: Fortran Note:
552: .vb
553: PetscInt, pointer :: cols(:)
554: .ve
555: Use `PETSC_NULL_INTEGER` if `n` is not needed
557: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringDestroy()`, `MatFDColoringView()`, `MatFDColoringApply()`
558: @*/
559: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring, PetscInt *n, const PetscInt *cols[])
560: {
561: PetscFunctionBegin;
562: if (coloring->currentcolor >= 0) {
563: *n = coloring->ncolumns[coloring->currentcolor];
564: *cols = coloring->columns[coloring->currentcolor];
565: } else {
566: *n = 0;
567: }
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: /*@
572: MatFDColoringApply - Given a matrix for which a `MatFDColoring` context
573: has been created, computes the Jacobian for a function via finite differences.
575: Collective
577: Input Parameters:
578: + J - matrix to store Jacobian entries into
579: . coloring - coloring context created with `MatFDColoringCreate()`
580: . x1 - location at which Jacobian is to be computed
581: - sctx - context required by function, if this is being used with the `SNES` solver then it is `SNES` object, otherwise it is `NULL`
583: Options Database Keys:
584: + -mat_fd_type - "wp" or "ds" (see `MATMFFD_WP` or `MATMFFD_DS`)
585: . -mat_fd_coloring_view - Activates basic viewing or coloring
586: . -mat_fd_coloring_view draw - Activates drawing of coloring
587: - -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info
589: Level: intermediate
591: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringDestroy()`, `MatFDColoringView()`, `MatFDColoringSetFunction()`, `MatFDColoringSetValues()`
592: @*/
593: PetscErrorCode MatFDColoringApply(Mat J, MatFDColoring coloring, Vec x1, void *sctx)
594: {
595: PetscBool eq;
597: PetscFunctionBegin;
601: PetscCall(PetscObjectCompareId((PetscObject)J, coloring->matid, &eq));
602: PetscCheck(eq, PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONG, "Matrix used with MatFDColoringApply() must be that used with MatFDColoringCreate()");
603: PetscCheck(coloring->f, PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONGSTATE, "Must call MatFDColoringSetFunction()");
604: PetscCheck(coloring->setupcalled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call MatFDColoringSetUp()");
606: PetscCall(MatSetUnfactored(J));
607: PetscCall(PetscLogEventBegin(MAT_FDColoringApply, coloring, J, x1, 0));
608: PetscUseTypeMethod(J, fdcoloringapply, coloring, x1, sctx);
609: PetscCall(PetscLogEventEnd(MAT_FDColoringApply, coloring, J, x1, 0));
610: if (!coloring->viewed) {
611: PetscCall(MatFDColoringViewFromOptions(coloring, NULL, "-mat_fd_coloring_view"));
612: coloring->viewed = PETSC_TRUE;
613: }
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }