Actual source code: snesj.c
1: #include <petsc/private/snesimpl.h>
2: #include <petsc/private/vecimpl.h>
3: #include <petscdm.h>
5: /*@
6: SNESComputeJacobianDefault - Computes the Jacobian using finite differences.
8: Collective
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()`
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`)
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 etol - 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 (wp|ds) - See `MATMFFD_WP` and `MATMFFD_DS`
26: Level: intermediate
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.
35: An alternative routine that uses coloring to exploit matrix sparsity is
36: `SNESComputeJacobianDefaultColor()`.
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()`.
41: This function can be provided to `SNESSetJacobian()` along with a dense matrix to hold the Jacobian
43: Developer Note:
44: The function has a poorly chosen name since it does not mention the use of finite differences
46: .seealso: [](ch_snes), `SNES`, `SNESSetJacobian()`, `SNESComputeJacobianDefaultColor()`, `MatCreateSNESMF()`
47: @*/
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;
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));
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: }
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) PetscCall(SNESComputeMFFunction(snes, x1, j1a));
93: else PetscCall(SNESComputeFunction(snes, x1, j1a)); /* does not handle use of SNESSetFunctionDomainError() correctly */
95: PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing options", "SNES");
96: PetscCall(PetscOptionsEList("-mat_fd_type", "Algorithm to compute difference parameter", "SNESComputeJacobianDefault", list, 2, "wp", &value, &flg));
97: PetscOptionsEnd();
98: if (flg && !value) use_wp = PETSC_FALSE;
100: if (use_wp) PetscCall(VecNorm(x1, NORM_2, &unorm));
101: /* Compute Jacobian approximation, 1 column at a time.
102: x1 = current iterate, j1a = F(x1)
103: x2 = perturbed iterate, j2a = F(x2)
104: */
105: for (i = 0; i < N; i++) {
106: PetscCall(VecCopy(x1, x2));
107: if (i >= start && i < end) {
108: PetscCall(VecGetArrayRead(x1, &xx));
109: if (use_wp) dx = PetscSqrtReal(1.0 + unorm);
110: else dx = xx[i - start];
111: PetscCall(VecRestoreArrayRead(x1, &xx));
112: if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par;
113: dx *= epsilon;
114: wscale = 1.0 / dx;
115: if (x2->ops->setvalues) PetscCall(VecSetValues(x2, 1, &i, &dx, ADD_VALUES));
116: else {
117: PetscCall(VecGetArray(x2, &y));
118: y[i - start] += dx;
119: PetscCall(VecRestoreArray(x2, &y));
120: }
121: } else {
122: wscale = 0.0;
123: }
124: PetscCall(VecAssemblyBegin(x2));
125: PetscCall(VecAssemblyEnd(x2));
126: if (dms->ops->computemffunction) PetscCall(SNESComputeMFFunction(snes, x2, j2a));
127: else PetscCall(SNESComputeFunction(snes, x2, j2a)); /* does not handle use of SNESSetFunctionDomainError() correctly */
128: PetscCall(VecAXPY(j2a, -1.0, j1a));
129: /* Communicate scale=1/dx_i to all processors */
130: PetscCall(VecGetOwnershipRanges(x1, &ranges));
131: root = size;
132: for (j = size - 1; j > -1; j--) {
133: root--;
134: if (i >= ranges[j]) break;
135: }
136: PetscCallMPI(MPI_Bcast(&wscale, 1, MPIU_SCALAR, root, comm));
137: PetscCall(VecScale(j2a, wscale));
138: PetscCall(VecNorm(j2a, NORM_INFINITY, &amax));
139: amax *= 1.e-14;
140: PetscCall(VecGetArray(j2a, &y));
141: for (j = start; j < end; j++) {
142: if (PetscAbsScalar(y[j - start]) > amax || j == i) PetscCall(MatSetValues(B, 1, &j, 1, &i, y + j - start, INSERT_VALUES));
143: }
144: PetscCall(VecRestoreArray(j2a, &y));
145: }
146: if (snes->dm) {
147: PetscCall(DMRestoreGlobalVector(snes->dm, &j1a));
148: PetscCall(DMRestoreGlobalVector(snes->dm, &j2a));
149: PetscCall(DMRestoreGlobalVector(snes->dm, &x2));
150: }
151: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
152: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
153: if (B != J) {
154: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
155: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
156: }
157: snes->max_funcs = max_funcs;
158: snes->nfuncs -= N;
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }