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    - 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`)

 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, void *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) {
 93:     PetscCall(SNESComputeMFFunction(snes, x1, j1a));
 94:   } else {
 95:     PetscCall(SNESComputeFunction(snes, x1, j1a));
 96:   }

 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;

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));
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: }