1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2 #include <petsc/private/vecimpl.h> /* for Vec->ops->setvalues */
3 #include <petscdm.h>
4
5 /*@
6 SNESComputeJacobianDefault - Computes the Jacobian using finite differences.
7
8 Collective
9
10 Input Parameters:
11 + snes - the `SNES` context
12 . x1 - compute Jacobian at this point
13 - ctx - application's function context, as set with `SNESSetFunction()`
14
15 Output Parameters:
16 + J - Jacobian matrix (not altered in this routine)
17 - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
18
19 Options Database Keys:
20 + -snes_fd - Activates `SNESComputeJacobianDefault()`
21 . -snes_fd_coloring - Activates a faster computation that uses a graph coloring of the matrix
22 . -snes_test_err - Square root of function error tolerance, default square root of machine
23 epsilon (1.e-8 in double, 3.e-4 in single)
24 - -mat_fd_type - Either wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
25
26 Level: intermediate
27
28 Notes:
29 This routine is slow and expensive, and is not currently optimized
30 to take advantage of sparsity in the problem. Although
31 `SNESComputeJacobianDefault()` is not recommended for general use
32 in large-scale applications, It can be useful in checking the
33 correctness of a user-provided Jacobian.
34
35 An alternative routine that uses coloring to exploit matrix sparsity is
36 `SNESComputeJacobianDefaultColor()`.
37
38 This routine ignores the maximum number of function evaluations set with `SNESSetTolerances()` and the function
39 evaluations it performs are not counted in what is returned by of `SNESGetNumberFunctionEvals()`.
40
41 This function can be provided to `SNESSetJacobian()` along with a dense matrix to hold the Jacobian
42
43 Developer Note:
44 The function has a poorly chosen name since it does not mention the use of finite differences
45
46 .seealso: [](ch_snes), `SNES`, `SNESSetJacobian()`, `SNESComputeJacobianDefaultColor()`, `MatCreateSNESMF()`
47 @*/
SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,PetscCtx ctx)48 PetscErrorCode SNESComputeJacobianDefault(SNES snes, Vec x1, Mat J, Mat B, PetscCtx ctx)
49 {
50 Vec j1a, j2a, x2;
51 PetscInt i, N, start, end, j, value, max_funcs = snes->max_funcs;
52 PetscScalar dx, *y, wscale;
53 const PetscScalar *xx;
54 PetscReal amax, epsilon = PETSC_SQRT_MACHINE_EPSILON;
55 PetscReal dx_min = 1.e-16, dx_par = 1.e-1, unorm;
56 MPI_Comm comm;
57 PetscBool assembled, use_wp = PETSC_TRUE, flg;
58 const char *list[2] = {"ds", "wp"};
59 PetscMPIInt size, root;
60 const PetscInt *ranges;
61 DM dm;
62 DMSNES dms;
63
64 PetscFunctionBegin;
65 snes->max_funcs = PETSC_INT_MAX;
66 /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */
67 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
68 PetscCall(PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_test_err", &epsilon, NULL));
69
70 PetscCall(PetscObjectGetComm((PetscObject)x1, &comm));
71 PetscCallMPI(MPI_Comm_size(comm, &size));
72 PetscCall(MatAssembled(B, &assembled));
73 if (assembled) PetscCall(MatZeroEntries(B));
74 if (!snes->nvwork) {
75 if (snes->dm) {
76 PetscCall(DMGetGlobalVector(snes->dm, &j1a));
77 PetscCall(DMGetGlobalVector(snes->dm, &j2a));
78 PetscCall(DMGetGlobalVector(snes->dm, &x2));
79 } else {
80 snes->nvwork = 3;
81 PetscCall(VecDuplicateVecs(x1, snes->nvwork, &snes->vwork));
82 j1a = snes->vwork[0];
83 j2a = snes->vwork[1];
84 x2 = snes->vwork[2];
85 }
86 }
87
88 PetscCall(VecGetSize(x1, &N));
89 PetscCall(VecGetOwnershipRange(x1, &start, &end));
90 PetscCall(SNESGetDM(snes, &dm));
91 PetscCall(DMGetDMSNES(dm, &dms));
92 if (dms->ops->computemffunction) {
93 PetscCall(SNESComputeMFFunction(snes, x1, j1a));
94 } else {
95 PetscCall(SNESComputeFunction(snes, x1, j1a)); /* does not handle use of SNESSetFunctionDomainError() correctly */
96 }
97
98 PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing options", "SNES");
99 PetscCall(PetscOptionsEList("-mat_fd_type", "Algorithm to compute difference parameter", "SNESComputeJacobianDefault", list, 2, "wp", &value, &flg));
100 PetscOptionsEnd();
101 if (flg && !value) use_wp = PETSC_FALSE;
102
103 if (use_wp) PetscCall(VecNorm(x1, NORM_2, &unorm));
104 /* Compute Jacobian approximation, 1 column at a time.
105 x1 = current iterate, j1a = F(x1)
106 x2 = perturbed iterate, j2a = F(x2)
107 */
108 for (i = 0; i < N; i++) {
109 PetscCall(VecCopy(x1, x2));
110 if (i >= start && i < end) {
111 PetscCall(VecGetArrayRead(x1, &xx));
112 if (use_wp) dx = PetscSqrtReal(1.0 + unorm);
113 else dx = xx[i - start];
114 PetscCall(VecRestoreArrayRead(x1, &xx));
115 if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par;
116 dx *= epsilon;
117 wscale = 1.0 / dx;
118 if (x2->ops->setvalues) PetscCall(VecSetValues(x2, 1, &i, &dx, ADD_VALUES));
119 else {
120 PetscCall(VecGetArray(x2, &y));
121 y[i - start] += dx;
122 PetscCall(VecRestoreArray(x2, &y));
123 }
124 } else {
125 wscale = 0.0;
126 }
127 PetscCall(VecAssemblyBegin(x2));
128 PetscCall(VecAssemblyEnd(x2));
129 if (dms->ops->computemffunction) {
130 PetscCall(SNESComputeMFFunction(snes, x2, j2a));
131 } else {
132 PetscCall(SNESComputeFunction(snes, x2, j2a)); /* does not handle use of SNESSetFunctionDomainError() correctly */
133 }
134 PetscCall(VecAXPY(j2a, -1.0, j1a));
135 /* Communicate scale=1/dx_i to all processors */
136 PetscCall(VecGetOwnershipRanges(x1, &ranges));
137 root = size;
138 for (j = size - 1; j > -1; j--) {
139 root--;
140 if (i >= ranges[j]) break;
141 }
142 PetscCallMPI(MPI_Bcast(&wscale, 1, MPIU_SCALAR, root, comm));
143 PetscCall(VecScale(j2a, wscale));
144 PetscCall(VecNorm(j2a, NORM_INFINITY, &amax));
145 amax *= 1.e-14;
146 PetscCall(VecGetArray(j2a, &y));
147 for (j = start; j < end; j++) {
148 if (PetscAbsScalar(y[j - start]) > amax || j == i) PetscCall(MatSetValues(B, 1, &j, 1, &i, y + j - start, INSERT_VALUES));
149 }
150 PetscCall(VecRestoreArray(j2a, &y));
151 }
152 if (snes->dm) {
153 PetscCall(DMRestoreGlobalVector(snes->dm, &j1a));
154 PetscCall(DMRestoreGlobalVector(snes->dm, &j2a));
155 PetscCall(DMRestoreGlobalVector(snes->dm, &x2));
156 }
157 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
158 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
159 if (B != J) {
160 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
161 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
162 }
163 snes->max_funcs = max_funcs;
164 snes->nfuncs -= N;
165 PetscFunctionReturn(PETSC_SUCCESS);
166 }
167