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: if (F) {
13: VecCopy(F, fd->w1);
14: fd->fset = PETSC_TRUE;
15: } else {
16: fd->fset = PETSC_FALSE;
17: }
18: return 0;
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: /* loop over colors */
30: nz = 0;
31: for (i = 0; i < fd->ncolors; i++) {
32: for (j = 0; j < fd->nrows[i]; j++) {
33: row = Jentry[nz].row;
34: y = fd->M - row - fd->rstart;
35: x = (PetscReal)Jentry[nz++].col;
36: PetscDrawRectangle(draw, x, y, x + 1, y + 1, i + 1, i + 1, i + 1, i + 1);
37: }
38: }
39: return 0;
40: }
42: static PetscErrorCode MatFDColoringView_Draw(MatFDColoring fd, PetscViewer viewer)
43: {
44: PetscBool isnull;
45: PetscDraw draw;
46: PetscReal xr, yr, xl, yl, h, w;
48: PetscViewerDrawGetDraw(viewer, 0, &draw);
49: PetscDrawIsNull(draw, &isnull);
50: if (isnull) return 0;
52: xr = fd->N;
53: yr = fd->M;
54: h = yr / 10.0;
55: w = xr / 10.0;
56: xr += w;
57: yr += h;
58: xl = -w;
59: yl = -h;
60: PetscDrawSetCoordinates(draw, xl, yl, xr, yr);
61: PetscObjectCompose((PetscObject)fd, "Zoomviewer", (PetscObject)viewer);
62: PetscDrawZoom(draw, MatFDColoringView_Draw_Zoom, fd);
63: PetscObjectCompose((PetscObject)fd, "Zoomviewer", NULL);
64: PetscDrawSave(draw);
65: return 0;
66: }
68: /*@C
69: MatFDColoringView - Views a finite difference coloring context.
71: Collective on c
73: Input Parameters:
74: + c - the coloring context
75: - viewer - visualization context
77: Level: intermediate
79: Notes:
80: The available visualization contexts include
81: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
82: . `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
83: output where only the first processor opens
84: the file. All other processors send their
85: data to the first processor to print.
86: - `PETSC_VIEWER_DRAW_WORLD` - graphical display of nonzero structure
88: Since PETSc uses only a small number of basic colors (currently 33), if the coloring
89: involves more than 33 then some seemingly identical colors are displayed making it look
90: like an illegal coloring. This is just a graphical artifact.
92: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`
93: @*/
94: PetscErrorCode MatFDColoringView(MatFDColoring c, PetscViewer viewer)
95: {
96: PetscInt i, j;
97: PetscBool isdraw, iascii;
98: PetscViewerFormat format;
101: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)c), &viewer);
105: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
106: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
107: if (isdraw) {
108: MatFDColoringView_Draw(c, viewer);
109: } else if (iascii) {
110: PetscObjectPrintClassNamePrefixType((PetscObject)c, viewer);
111: PetscViewerASCIIPrintf(viewer, " Error tolerance=%g\n", (double)c->error_rel);
112: PetscViewerASCIIPrintf(viewer, " Umin=%g\n", (double)c->umin);
113: PetscViewerASCIIPrintf(viewer, " Number of colors=%" PetscInt_FMT "\n", c->ncolors);
115: PetscViewerGetFormat(viewer, &format);
116: if (format != PETSC_VIEWER_ASCII_INFO) {
117: PetscInt row, col, nz;
118: nz = 0;
119: for (i = 0; i < c->ncolors; i++) {
120: PetscViewerASCIIPrintf(viewer, " Information for color %" PetscInt_FMT "\n", i);
121: PetscViewerASCIIPrintf(viewer, " Number of columns %" PetscInt_FMT "\n", c->ncolumns[i]);
122: for (j = 0; j < c->ncolumns[i]; j++) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT "\n", c->columns[i][j]);
123: PetscViewerASCIIPrintf(viewer, " Number of rows %" PetscInt_FMT "\n", c->nrows[i]);
124: if (c->matentry) {
125: for (j = 0; j < c->nrows[i]; j++) {
126: row = c->matentry[nz].row;
127: col = c->matentry[nz++].col;
128: PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT " \n", row, col);
129: }
130: }
131: }
132: }
133: PetscViewerFlush(viewer);
134: }
135: return 0;
136: }
138: /*@
139: MatFDColoringSetParameters - Sets the parameters for the sparse approximation of
140: a Jacobian matrix using finite differences.
142: Logically Collective
144: The Jacobian is estimated with the differencing approximation
145: .vb
146: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
147: htype = 'ds':
148: h = error_rel*u[i] if abs(u[i]) > umin
149: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
150: dx_{i} = (0, ... 1, .... 0)
152: htype = 'wp':
153: h = error_rel * sqrt(1 + ||u||)
154: .ve
156: Input Parameters:
157: + matfd - the coloring context
158: . error - relative error
159: - umin - minimum allowable u-value magnitude
161: Level: advanced
163: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringSetFromOptions()`
164: @*/
165: PetscErrorCode MatFDColoringSetParameters(MatFDColoring matfd, PetscReal error, PetscReal umin)
166: {
170: if (error != PETSC_DEFAULT) matfd->error_rel = error;
171: if (umin != PETSC_DEFAULT) matfd->umin = umin;
172: return 0;
173: }
175: /*@
176: MatFDColoringSetBlockSize - Sets block size for efficient inserting entries of Jacobian matrix.
178: Logically Collective
180: Input Parameters:
181: + coloring - the coloring context
182: . brows - number of rows in the block
183: - bcols - number of columns in the block
185: Level: intermediate
187: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringSetFromOptions()`
188: @*/
189: PetscErrorCode MatFDColoringSetBlockSize(MatFDColoring matfd, PetscInt brows, PetscInt bcols)
190: {
194: if (brows != PETSC_DEFAULT) matfd->brows = brows;
195: if (bcols != PETSC_DEFAULT) matfd->bcols = bcols;
196: return 0;
197: }
199: /*@
200: MatFDColoringSetUp - Sets up the internal data structures of matrix coloring context for the later use.
202: Collective
204: Input Parameters:
205: + mat - the matrix containing the nonzero structure of the Jacobian
206: . iscoloring - the coloring of the matrix; usually obtained with `MatGetColoring()` or `DMCreateColoring()`
207: - color - the matrix coloring context
209: Level: beginner
211: Notes:
212: When the coloring type is `IS_COLORING_LOCAL` the coloring is in the local ordering of the unknowns.
214: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringDestroy()`
215: @*/
216: PetscErrorCode MatFDColoringSetUp(Mat mat, ISColoring iscoloring, MatFDColoring color)
217: {
218: PetscBool eq;
222: if (color->setupcalled) return 0;
223: PetscObjectCompareId((PetscObject)mat, color->matid, &eq);
226: PetscLogEventBegin(MAT_FDColoringSetUp, mat, 0, 0, 0);
227: PetscUseTypeMethod(mat, fdcoloringsetup, iscoloring, color);
229: color->setupcalled = PETSC_TRUE;
230: PetscLogEventEnd(MAT_FDColoringSetUp, mat, 0, 0, 0);
231: return 0;
232: }
234: /*@C
235: MatFDColoringGetFunction - Gets the function to use for computing the Jacobian.
237: Not Collective
239: Input Parameter:
240: . coloring - the coloring context
242: Output Parameters:
243: + f - the function
244: - fctx - the optional user-defined function context
246: Level: intermediate
248: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`, `MatFDColoringSetFromOptions()`
249: @*/
250: PetscErrorCode MatFDColoringGetFunction(MatFDColoring matfd, PetscErrorCode (**f)(void), void **fctx)
251: {
253: if (f) *f = matfd->f;
254: if (fctx) *fctx = matfd->fctx;
255: return 0;
256: }
258: /*@C
259: MatFDColoringSetFunction - Sets the function to use for computing the Jacobian.
261: Logically Collective
263: Input Parameters:
264: + coloring - the coloring context
265: . f - the function
266: - fctx - the optional user-defined function context
268: Calling sequence of (*f) function:
269: For `SNES`: PetscErrorCode (*f)(SNES,Vec,Vec,void*)
270: If not using `SNES`: PetscErrorCode (*f)(void *dummy,Vec,Vec,void*) and dummy is ignored
272: Level: advanced
274: Note:
275: This function is usually used automatically by `SNES` (when one uses `SNESSetJacobian()` with the argument
276: `SNESComputeJacobianDefaultColor()`) and only needs to be used by someone computing a matrix via coloring directly by
277: calling `MatFDColoringApply()`
279: Fortran Note:
280: In Fortran you must call `MatFDColoringSetFunction()` for a coloring object to
281: be used without `SNES` or within the `SNES` solvers.
283: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringGetFunction()`, `MatFDColoringSetFromOptions()`
284: @*/
285: PetscErrorCode MatFDColoringSetFunction(MatFDColoring matfd, PetscErrorCode (*f)(void), void *fctx)
286: {
288: matfd->f = f;
289: matfd->fctx = fctx;
290: return 0;
291: }
293: /*@
294: MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
295: the options database.
297: Collective
299: The Jacobian, F'(u), is estimated with the differencing approximation
300: .vb
301: F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
302: h = error_rel*u[i] if abs(u[i]) > umin
303: = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i]
304: dx_{i} = (0, ... 1, .... 0)
305: .ve
307: Input Parameter:
308: . coloring - the coloring context
310: Options Database Keys:
311: + -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
312: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
313: . -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
314: . -mat_fd_coloring_view - Activates basic viewing
315: . -mat_fd_coloring_view ::ascii_info - Activates viewing info
316: - -mat_fd_coloring_view draw - Activates drawing
318: Level: intermediate
320: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringView()`, `MatFDColoringSetParameters()`
321: @*/
322: PetscErrorCode MatFDColoringSetFromOptions(MatFDColoring matfd)
323: {
324: PetscBool flg;
325: char value[3];
329: PetscObjectOptionsBegin((PetscObject)matfd);
330: PetscOptionsReal("-mat_fd_coloring_err", "Square root of relative error in function", "MatFDColoringSetParameters", matfd->error_rel, &matfd->error_rel, NULL);
331: PetscOptionsReal("-mat_fd_coloring_umin", "Minimum allowable u magnitude", "MatFDColoringSetParameters", matfd->umin, &matfd->umin, NULL);
332: PetscOptionsString("-mat_fd_type", "Algorithm to compute h, wp or ds", "MatFDColoringCreate", matfd->htype, value, sizeof(value), &flg);
333: if (flg) {
334: if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
335: else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
336: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown finite differencing type %s", value);
337: }
338: PetscOptionsInt("-mat_fd_coloring_brows", "Number of block rows", "MatFDColoringSetBlockSize", matfd->brows, &matfd->brows, NULL);
339: PetscOptionsInt("-mat_fd_coloring_bcols", "Number of block columns", "MatFDColoringSetBlockSize", matfd->bcols, &matfd->bcols, &flg);
340: if (flg && matfd->bcols > matfd->ncolors) {
341: /* input bcols cannot be > matfd->ncolors, thus set it as ncolors */
342: matfd->bcols = matfd->ncolors;
343: }
345: /* process any options handlers added with PetscObjectAddOptionsHandler() */
346: PetscObjectProcessOptionsHandlers((PetscObject)matfd, PetscOptionsObject);
347: PetscOptionsEnd();
348: return 0;
349: }
351: /*@C
352: MatFDColoringSetType - Sets the approach for computing the finite difference parameter
354: Collective
356: Input Parameters:
357: + coloring - the coloring context
358: - type - either `MATMFFD_WP` or `MATMFFD_DS`
360: Options Database Keys:
361: . -mat_fd_type - "wp" or "ds"
363: Level: intermediate
365: Note:
366: It is goofy that the argument type is `MatMFFDType` since the `MatFDColoring` actually computes the matrix entries
367: but the process of computing the entries is the same as as with the `MATMFFD` operation so we should reuse the names instead of
368: introducing another one.
370: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringView()`, `MatFDColoringSetParameters()`
371: @*/
372: PetscErrorCode MatFDColoringSetType(MatFDColoring matfd, MatMFFDType type)
373: {
375: /*
376: It is goofy to handle the strings this way but currently there is no code to free a dynamically created matfd->htype
377: and this function is being provided as patch for a release so it shouldn't change the implementation
378: */
379: if (type[0] == 'w' && type[1] == 'p') matfd->htype = "wp";
380: else if (type[0] == 'd' && type[1] == 's') matfd->htype = "ds";
381: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown finite differencing type %s", type);
382: return 0;
383: }
385: PetscErrorCode MatFDColoringViewFromOptions(MatFDColoring fd, const char prefix[], const char optionname[])
386: {
387: PetscBool flg;
388: PetscViewer viewer;
389: PetscViewerFormat format;
391: if (prefix) {
392: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd), ((PetscObject)fd)->options, prefix, optionname, &viewer, &format, &flg);
393: } else {
394: PetscOptionsGetViewer(PetscObjectComm((PetscObject)fd), ((PetscObject)fd)->options, ((PetscObject)fd)->prefix, optionname, &viewer, &format, &flg);
395: }
396: if (flg) {
397: PetscViewerPushFormat(viewer, format);
398: MatFDColoringView(fd, viewer);
399: PetscViewerPopFormat(viewer);
400: PetscViewerDestroy(&viewer);
401: }
402: return 0;
403: }
405: /*@
406: MatFDColoringCreate - Creates a matrix coloring context for finite difference
407: computation of Jacobians.
409: Collective
411: Input Parameters:
412: + mat - the matrix containing the nonzero structure of the Jacobian
413: - iscoloring - the coloring of the matrix; usually obtained with `MatColoringCreate()` or `DMCreateColoring()`
415: Output Parameter:
416: . color - the new coloring context
418: Level: intermediate
420: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringDestroy()`, `SNESComputeJacobianDefaultColor()`, `ISColoringCreate()`,
421: `MatFDColoringSetFunction()`, `MatFDColoringSetFromOptions()`, `MatFDColoringApply()`,
422: `MatFDColoringView()`, `MatFDColoringSetParameters()`, `MatColoringCreate()`, `DMCreateColoring()`, `MatFDColoringSetValues()`
423: @*/
424: PetscErrorCode MatFDColoringCreate(Mat mat, ISColoring iscoloring, MatFDColoring *color)
425: {
426: MatFDColoring c;
427: MPI_Comm comm;
428: PetscInt M, N;
432: PetscLogEventBegin(MAT_FDColoringCreate, mat, 0, 0, 0);
433: MatGetSize(mat, &M, &N);
435: PetscObjectGetComm((PetscObject)mat, &comm);
436: PetscHeaderCreate(c, MAT_FDCOLORING_CLASSID, "MatFDColoring", "Jacobian computation via finite differences with coloring", "Mat", comm, MatFDColoringDestroy, MatFDColoringView);
438: c->ctype = iscoloring->ctype;
439: PetscObjectGetId((PetscObject)mat, &c->matid);
441: PetscUseTypeMethod(mat, fdcoloringcreate, iscoloring, c);
443: MatCreateVecs(mat, NULL, &c->w1);
444: /* Vec is used intensively in particular piece of scalar CPU code; won't benefit from bouncing back and forth to the GPU */
445: VecBindToCPU(c->w1, PETSC_TRUE);
446: VecDuplicate(c->w1, &c->w2);
447: /* Vec is used intensively in particular piece of scalar CPU code; won't benefit from bouncing back and forth to the GPU */
448: VecBindToCPU(c->w2, PETSC_TRUE);
450: c->error_rel = PETSC_SQRT_MACHINE_EPSILON;
451: c->umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
452: c->currentcolor = -1;
453: c->htype = "wp";
454: c->fset = PETSC_FALSE;
455: c->setupcalled = PETSC_FALSE;
457: *color = c;
458: PetscObjectCompose((PetscObject)mat, "SNESMatFDColoring", (PetscObject)c);
459: PetscLogEventEnd(MAT_FDColoringCreate, mat, 0, 0, 0);
460: return 0;
461: }
463: /*@
464: MatFDColoringDestroy - Destroys a matrix coloring context that was created
465: via `MatFDColoringCreate()`.
467: Collective
469: Input Parameter:
470: . c - coloring context
472: Level: intermediate
474: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`
475: @*/
476: PetscErrorCode MatFDColoringDestroy(MatFDColoring *c)
477: {
478: PetscInt i;
479: MatFDColoring color = *c;
481: if (!*c) return 0;
482: if (--((PetscObject)color)->refct > 0) {
483: *c = NULL;
484: return 0;
485: }
487: /* we do not free the column arrays since their entries are owned by the ISs in color->isa */
488: for (i = 0; i < color->ncolors; i++) ISDestroy(&color->isa[i]);
489: PetscFree(color->isa);
490: PetscFree2(color->ncolumns, color->columns);
491: PetscFree(color->nrows);
492: if (color->htype[0] == 'w') {
493: PetscFree(color->matentry2);
494: } else {
495: PetscFree(color->matentry);
496: }
497: PetscFree(color->dy);
498: if (color->vscale) VecDestroy(&color->vscale);
499: VecDestroy(&color->w1);
500: VecDestroy(&color->w2);
501: VecDestroy(&color->w3);
502: PetscHeaderDestroy(c);
503: return 0;
504: }
506: /*@C
507: MatFDColoringGetPerturbedColumns - Returns the indices of the columns that
508: that are currently being perturbed.
510: Not Collective
512: Input Parameter:
513: . coloring - coloring context created with `MatFDColoringCreate()`
515: Output Parameters:
516: + n - the number of local columns being perturbed
517: - cols - the column indices, in global numbering
519: Level: advanced
521: Note:
522: IF the matrix type is `MATBAIJ`, then the block column indices are returned
524: Fortran Note:
525: This routine has a different interface for Fortran
526: .vb
527: #include <petsc/finclude/petscmat.h>
528: use petscmat
529: PetscInt, pointer :: array(:)
530: PetscErrorCode ierr
531: MatFDColoring i
532: call MatFDColoringGetPerturbedColumnsF90(i,array,ierr)
533: use the entries of array ...
534: call MatFDColoringRestorePerturbedColumnsF90(i,array,ierr)
535: .ve
537: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringDestroy()`, `MatFDColoringView()`, `MatFDColoringApply()`
538: @*/
539: PetscErrorCode MatFDColoringGetPerturbedColumns(MatFDColoring coloring, PetscInt *n, const PetscInt *cols[])
540: {
541: if (coloring->currentcolor >= 0) {
542: *n = coloring->ncolumns[coloring->currentcolor];
543: *cols = coloring->columns[coloring->currentcolor];
544: } else {
545: *n = 0;
546: }
547: return 0;
548: }
550: /*@
551: MatFDColoringApply - Given a matrix for which a `MatFDColoring` context
552: has been created, computes the Jacobian for a function via finite differences.
554: Collective
556: Input Parameters:
557: + mat - location to store Jacobian
558: . coloring - coloring context created with `MatFDColoringCreate()`
559: . x1 - location at which Jacobian is to be computed
560: - sctx - context required by function, if this is being used with the SNES solver then it is `SNES` object, otherwise it is null
562: Options Database Keys:
563: + -mat_fd_type - "wp" or "ds" (see `MATMFFD_WP` or `MATMFFD_DS`)
564: . -mat_fd_coloring_view - Activates basic viewing or coloring
565: . -mat_fd_coloring_view draw - Activates drawing of coloring
566: - -mat_fd_coloring_view ::ascii_info - Activates viewing of coloring info
568: Level: intermediate
570: .seealso: `Mat`, `MatFDColoring`, `MatFDColoringCreate()`, `MatFDColoringDestroy()`, `MatFDColoringView()`, `MatFDColoringSetFunction()`, `MatFDColoringSetValues()`
571: @*/
572: PetscErrorCode MatFDColoringApply(Mat J, MatFDColoring coloring, Vec x1, void *sctx)
573: {
574: PetscBool eq;
579: PetscObjectCompareId((PetscObject)J, coloring->matid, &eq);
584: MatSetUnfactored(J);
585: PetscLogEventBegin(MAT_FDColoringApply, coloring, J, x1, 0);
586: PetscUseTypeMethod(J, fdcoloringapply, coloring, x1, sctx);
587: PetscLogEventEnd(MAT_FDColoringApply, coloring, J, x1, 0);
588: if (!coloring->viewed) {
589: MatFDColoringViewFromOptions(coloring, NULL, "-mat_fd_coloring_view");
590: coloring->viewed = PETSC_TRUE;
591: }
592: return 0;
593: }