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