Actual source code: fdiff.c

  1: #include <petsctao.h>
  2: #include <petsc/private/taoimpl.h>
  3: #include <petscsnes.h>
  4: #include <petscdmshell.h>

  6: /*
  7:    For finited difference computations of the Hessian, we use PETSc's SNESComputeJacobianDefault
  8: */
  9: static PetscErrorCode Fsnes(SNES snes, Vec X, Vec G, void *ctx)
 10: {
 11:   Tao tao = (Tao)ctx;

 14:   TaoComputeGradient(tao, X, G);
 15:   return 0;
 16: }

 18: /*@C
 19:   TaoDefaultComputeGradient - computes the gradient using finite differences.

 21:   Collective

 23:   Input Parameters:
 24: + tao   - the Tao context
 25: . X     - compute gradient at this point
 26: - dummy - not used

 28:   Output Parameter:
 29: . G - Gradient Vector

 31:   Options Database Key:
 32: + -tao_fd_gradient      - activates TaoDefaultComputeGradient()
 33: - -tao_fd_delta <delta> - change in X used to calculate finite differences

 35:   Level: advanced

 37:   Notes:
 38:   This routine is slow and expensive, and is not optimized
 39:   to take advantage of sparsity in the problem.  Although
 40:   not recommended for general use
 41:   in large-scale applications, it can be useful in checking the
 42:   correctness of a user-provided gradient.  Use the tao method TAOTEST
 43:   to get an indication of whether your gradient is correct.
 44:   This finite difference gradient evaluation can be set using the routine `TaoSetGradient()` or by using the command line option -tao_fd_gradient

 46: .seealso: `Tao`, `TaoSetGradient()`
 47: @*/
 48: PetscErrorCode TaoDefaultComputeGradient(Tao tao, Vec Xin, Vec G, void *dummy)
 49: {
 50:   Vec          X;
 51:   PetscScalar *g;
 52:   PetscReal    f, f2;
 53:   PetscInt     low, high, N, i;
 54:   PetscBool    flg;
 55:   PetscReal    h = .5 * PETSC_SQRT_MACHINE_EPSILON;

 57:   PetscOptionsGetReal(((PetscObject)tao)->options, ((PetscObject)tao)->prefix, "-tao_fd_delta", &h, &flg);
 58:   VecDuplicate(Xin, &X);
 59:   VecCopy(Xin, X);
 60:   VecGetSize(X, &N);
 61:   VecGetOwnershipRange(X, &low, &high);
 62:   VecSetOption(X, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
 63:   VecGetArray(G, &g);
 64:   for (i = 0; i < N; i++) {
 65:     VecSetValue(X, i, -h, ADD_VALUES);
 66:     VecAssemblyBegin(X);
 67:     VecAssemblyEnd(X);
 68:     TaoComputeObjective(tao, X, &f);
 69:     VecSetValue(X, i, 2.0 * h, ADD_VALUES);
 70:     VecAssemblyBegin(X);
 71:     VecAssemblyEnd(X);
 72:     TaoComputeObjective(tao, X, &f2);
 73:     VecSetValue(X, i, -h, ADD_VALUES);
 74:     VecAssemblyBegin(X);
 75:     VecAssemblyEnd(X);
 76:     if (i >= low && i < high) g[i - low] = (f2 - f) / (2.0 * h);
 77:   }
 78:   VecRestoreArray(G, &g);
 79:   VecDestroy(&X);
 80:   return 0;
 81: }

 83: /*@C
 84:    TaoDefaultComputeHessian - Computes the Hessian using finite differences.

 86:    Collective

 88:    Input Parameters:
 89: +  tao   - the Tao context
 90: .  V     - compute Hessian at this point
 91: -  dummy - not used

 93:    Output Parameters:
 94: +  H - Hessian matrix (not altered in this routine)
 95: -  B - newly computed Hessian matrix to use with preconditioner (generally the same as H)

 97:    Options Database Key:
 98: .  -tao_fd_hessian - activates TaoDefaultComputeHessian()

100:    Level: advanced

102:    Notes:
103:    This routine is slow and expensive, and is not optimized
104:    to take advantage of sparsity in the problem.  Although
105:    it is not recommended for general use
106:    in large-scale applications, It can be useful in checking the
107:    correctness of a user-provided Hessian.

109: .seealso: `Tao`, `TaoSetHessian()`, `TaoDefaultComputeHessianColor()`, `SNESComputeJacobianDefault()`, `TaoSetGradient()`, `TaoDefaultComputeGradient()`
110: @*/
111: PetscErrorCode TaoDefaultComputeHessian(Tao tao, Vec V, Mat H, Mat B, void *dummy)
112: {
113:   SNES snes;
114:   DM   dm;

116:   PetscInfo(tao, "TAO Using finite differences w/o coloring to compute Hessian matrix\n");
117:   SNESCreate(PetscObjectComm((PetscObject)H), &snes);
118:   SNESSetFunction(snes, NULL, Fsnes, tao);
119:   SNESGetDM(snes, &dm);
120:   DMShellSetGlobalVector(dm, V);
121:   SNESSetUp(snes);
122:   if (H) {
123:     PetscInt n, N;

125:     VecGetSize(V, &N);
126:     VecGetLocalSize(V, &n);
127:     MatSetSizes(H, n, n, N, N);
128:     MatSetUp(H);
129:   }
130:   if (B && B != H) {
131:     PetscInt n, N;

133:     VecGetSize(V, &N);
134:     VecGetLocalSize(V, &n);
135:     MatSetSizes(B, n, n, N, N);
136:     MatSetUp(B);
137:   }
138:   SNESComputeJacobianDefault(snes, V, H, B, NULL);
139:   SNESDestroy(&snes);
140:   return 0;
141: }

143: /*@C
144:    TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.

146:    Collective

148:    Input Parameters:
149: +  tao - the Tao context
150: .  V   - compute Hessian at this point
151: -  ctx - the color object of type `MatFDColoring`

153:    Output Parameters:
154: +  H - Hessian matrix (not altered in this routine)
155: -  B - newly computed Hessian matrix to use with preconditioner (generally the same as H)

157:    Level: advanced

159: .seealso: `Tao`, `MatColoring`, `TaoSetHessian()`, `TaoDefaultComputeHessian()`, `SNESComputeJacobianDefaultColor()`, `TaoSetGradient()`
160: @*/
161: PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H, Mat B, void *ctx)
162: {
163:   MatFDColoring coloring = (MatFDColoring)ctx;

166:   PetscInfo(tao, "TAO computing matrix using finite differences Hessian and coloring\n");
167:   MatFDColoringApply(B, coloring, V, ctx);
168:   if (H != B) {
169:     MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
170:     MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
171:   }
172:   return 0;
173: }

175: PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao, Vec X, Mat H, Mat B, void *ctx)
176: {
177:   PetscInt  n, N;
178:   PetscBool assembled;

181:   MatAssembled(H, &assembled);
182:   if (!assembled) {
183:     VecGetSize(X, &N);
184:     VecGetLocalSize(X, &n);
185:     MatSetSizes(H, n, n, N, N);
186:     MatSetType(H, MATMFFD);
187:     MatSetUp(H);
188:     MatMFFDSetFunction(H, (PetscErrorCode(*)(void *, Vec, Vec))TaoComputeGradient, tao);
189:   }
190:   MatMFFDSetBase(H, X, NULL);
191:   MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
192:   MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
193:   return 0;
194: }