Actual source code: snesj.c
petsc-3.6.4 2016-04-12
2: #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
6: /*@C
7: SNESComputeJacobianDefault - Computes the Jacobian using finite differences.
9: Collective on SNES
11: Input Parameters:
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 Key:
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: Level: intermediate
37: .keywords: SNES, finite differences, Jacobian
39: .seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF()
40: @*/
41: PetscErrorCode SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
42: {
43: Vec j1a,j2a,x2;
44: PetscErrorCode ierr;
45: PetscInt i,N,start,end,j,value,root;
46: PetscScalar dx,*y,wscale;
47: const PetscScalar *xx;
48: PetscReal amax,epsilon = PETSC_SQRT_MACHINE_EPSILON;
49: PetscReal dx_min = 1.e-16,dx_par = 1.e-1,unorm;
50: MPI_Comm comm;
51: PetscBool assembled,use_wp = PETSC_TRUE,flg;
52: const char *list[2] = {"ds","wp"};
53: PetscMPIInt size;
54: const PetscInt *ranges;
57: PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);
59: PetscObjectGetComm((PetscObject)x1,&comm);
60: MPI_Comm_size(comm,&size);
61: MatAssembled(B,&assembled);
62: if (assembled) {
63: MatZeroEntries(B);
64: }
65: if (!snes->nvwork) {
66: snes->nvwork = 3;
68: VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);
69: PetscLogObjectParents(snes,snes->nvwork,snes->vwork);
70: }
71: j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2];
73: VecGetSize(x1,&N);
74: VecGetOwnershipRange(x1,&start,&end);
75: SNESComputeFunction(snes,x1,j1a);
77: PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES");
78: PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg);
79: PetscOptionsEnd();
80: if (flg && !value) use_wp = PETSC_FALSE;
82: if (use_wp) {
83: VecNorm(x1,NORM_2,&unorm);
84: }
85: /* Compute Jacobian approximation, 1 column at a time.
86: x1 = current iterate, j1a = F(x1)
87: x2 = perturbed iterate, j2a = F(x2)
88: */
89: for (i=0; i<N; i++) {
90: VecCopy(x1,x2);
91: if (i>= start && i<end) {
92: VecGetArrayRead(x1,&xx);
93: if (use_wp) dx = PetscSqrtReal(1.0 + unorm);
94: else dx = xx[i-start];
95: VecRestoreArrayRead(x1,&xx);
96: if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par;
97: dx *= epsilon;
98: wscale = 1.0/dx;
99: VecSetValues(x2,1,&i,&dx,ADD_VALUES);
100: } else {
101: wscale = 0.0;
102: }
103: VecAssemblyBegin(x2);
104: VecAssemblyEnd(x2);
105: SNESComputeFunction(snes,x2,j2a);
106: VecAXPY(j2a,-1.0,j1a);
107: /* Communicate scale=1/dx_i to all processors */
108: VecGetOwnershipRanges(x1,&ranges);
109: root = size;
110: for (j=size-1; j>-1; j--) {
111: root--;
112: if (i>=ranges[j]) break;
113: }
114: MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm);
116: VecScale(j2a,wscale);
117: VecNorm(j2a,NORM_INFINITY,&amax); amax *= 1.e-14;
118: VecGetArray(j2a,&y);
119: for (j=start; j<end; j++) {
120: if (PetscAbsScalar(y[j-start]) > amax || j == i) {
121: MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES);
122: }
123: }
124: VecRestoreArray(j2a,&y);
125: }
126: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
128: if (B != J) {
129: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
131: }
132: return(0);
133: }