Actual source code: snesj.c
2: #include <petsc/private/snesimpl.h>
3: #include <petscdm.h>
5: /*@C
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_test_err - Square root of function error tolerance, default square root of machine
22: epsilon (1.e-8 in double, 3.e-4 in single)
23: - -mat_fd_type - Either wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
25: Notes:
26: This routine is slow and expensive, and is not currently optimized
27: to take advantage of sparsity in the problem. Although
28: `SNESComputeJacobianDefault()` is not recommended for general use
29: in large-scale applications, It can be useful in checking the
30: correctness of a user-provided Jacobian.
32: An alternative routine that uses coloring to exploit matrix sparsity is
33: `SNESComputeJacobianDefaultColor()`.
35: This routine ignores the maximum number of function evaluations set with `SNESSetTolerances()` and the function
36: evaluations it performs are not counted in what is returned by of `SNESGetNumberFunctionEvals()`.
38: This function can be provided to `SNESSetJacobian()` along with a dense matrix to hold the Jacobian
40: Level: intermediate
42: .seealso: `SNES`, `SNESSetJacobian()`, `SNESSetJacobian()`, `SNESComputeJacobianDefaultColor()`, `MatCreateSNESMF()`
43: @*/
44: PetscErrorCode SNESComputeJacobianDefault(SNES snes, Vec x1, Mat J, Mat B, void *ctx)
45: {
46: Vec j1a, j2a, x2;
47: PetscInt i, N, start, end, j, value, root, max_funcs = snes->max_funcs;
48: PetscScalar dx, *y, wscale;
49: const PetscScalar *xx;
50: PetscReal amax, epsilon = PETSC_SQRT_MACHINE_EPSILON;
51: PetscReal dx_min = 1.e-16, dx_par = 1.e-1, unorm;
52: MPI_Comm comm;
53: PetscBool assembled, use_wp = PETSC_TRUE, flg;
54: const char *list[2] = {"ds", "wp"};
55: PetscMPIInt size;
56: const PetscInt *ranges;
57: DM dm;
58: DMSNES dms;
60: snes->max_funcs = PETSC_MAX_INT;
61: /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */
62: MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
63: PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_test_err", &epsilon, NULL);
65: PetscObjectGetComm((PetscObject)x1, &comm);
66: MPI_Comm_size(comm, &size);
67: MatAssembled(B, &assembled);
68: if (assembled) MatZeroEntries(B);
69: if (!snes->nvwork) {
70: if (snes->dm) {
71: DMGetGlobalVector(snes->dm, &j1a);
72: DMGetGlobalVector(snes->dm, &j2a);
73: DMGetGlobalVector(snes->dm, &x2);
74: } else {
75: snes->nvwork = 3;
76: VecDuplicateVecs(x1, snes->nvwork, &snes->vwork);
77: j1a = snes->vwork[0];
78: j2a = snes->vwork[1];
79: x2 = snes->vwork[2];
80: }
81: }
83: VecGetSize(x1, &N);
84: VecGetOwnershipRange(x1, &start, &end);
85: SNESGetDM(snes, &dm);
86: DMGetDMSNES(dm, &dms);
87: if (dms->ops->computemffunction) {
88: SNESComputeMFFunction(snes, x1, j1a);
89: } else {
90: SNESComputeFunction(snes, x1, j1a);
91: }
93: PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing options", "SNES");
94: PetscOptionsEList("-mat_fd_type", "Algorithm to compute difference parameter", "SNESComputeJacobianDefault", list, 2, "wp", &value, &flg);
95: PetscOptionsEnd();
96: if (flg && !value) use_wp = PETSC_FALSE;
98: if (use_wp) VecNorm(x1, NORM_2, &unorm);
99: /* Compute Jacobian approximation, 1 column at a time.
100: x1 = current iterate, j1a = F(x1)
101: x2 = perturbed iterate, j2a = F(x2)
102: */
103: for (i = 0; i < N; i++) {
104: VecCopy(x1, x2);
105: if (i >= start && i < end) {
106: VecGetArrayRead(x1, &xx);
107: if (use_wp) dx = PetscSqrtReal(1.0 + unorm);
108: else dx = xx[i - start];
109: VecRestoreArrayRead(x1, &xx);
110: if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par;
111: dx *= epsilon;
112: wscale = 1.0 / dx;
113: VecSetValues(x2, 1, &i, &dx, ADD_VALUES);
114: } else {
115: wscale = 0.0;
116: }
117: VecAssemblyBegin(x2);
118: VecAssemblyEnd(x2);
119: if (dms->ops->computemffunction) {
120: SNESComputeMFFunction(snes, x2, j2a);
121: } else {
122: SNESComputeFunction(snes, x2, j2a);
123: }
124: VecAXPY(j2a, -1.0, j1a);
125: /* Communicate scale=1/dx_i to all processors */
126: VecGetOwnershipRanges(x1, &ranges);
127: root = size;
128: for (j = size - 1; j > -1; j--) {
129: root--;
130: if (i >= ranges[j]) break;
131: }
132: MPI_Bcast(&wscale, 1, MPIU_SCALAR, root, comm);
133: VecScale(j2a, wscale);
134: VecNorm(j2a, NORM_INFINITY, &amax);
135: amax *= 1.e-14;
136: VecGetArray(j2a, &y);
137: for (j = start; j < end; j++) {
138: if (PetscAbsScalar(y[j - start]) > amax || j == i) MatSetValues(B, 1, &j, 1, &i, y + j - start, INSERT_VALUES);
139: }
140: VecRestoreArray(j2a, &y);
141: }
142: if (snes->dm) {
143: DMRestoreGlobalVector(snes->dm, &j1a);
144: DMRestoreGlobalVector(snes->dm, &j2a);
145: DMRestoreGlobalVector(snes->dm, &x2);
146: }
147: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
148: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
149: if (B != J) {
150: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
151: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
152: }
153: snes->max_funcs = max_funcs;
154: snes->nfuncs -= N;
155: return 0;
156: }