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