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 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(PetscOptionsRestoreViewer(&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: }