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 finite difference computations of the Hessian, we use PETSc's SNESComputeJacobianDefault
8: */
9: static PetscErrorCode Fsnes(SNES snes, Vec X, Vec G, PetscCtx ctx)
10: {
11: Tao tao = (Tao)ctx;
13: PetscFunctionBegin;
15: PetscCall(TaoComputeGradient(tao, X, G));
16: PetscFunctionReturn(PETSC_SUCCESS);
17: }
19: /*@C
20: TaoDefaultComputeGradient - computes the gradient using finite differences.
22: Collective
24: Input Parameters:
25: + tao - the Tao context
26: . Xin - compute gradient at this point
27: - dummy - not used
29: Output Parameter:
30: . G - Gradient Vector
32: Options Database Key:
33: + -tao_fd_gradient - activates TaoDefaultComputeGradient()
34: - -tao_fd_delta delta - change in X used to calculate finite differences
36: Level: advanced
38: Notes:
39: This routine is slow and expensive, and is not optimized
40: to take advantage of sparsity in the problem. Although
41: not recommended for general use
42: in large-scale applications, it can be useful in checking the
43: correctness of a user-provided gradient using the command-line option `-tao_test_gradient`
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()`, `TaoTermComputeGradientFD()`
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: PetscFunctionBegin;
58: PetscCall(PetscOptionsGetReal(((PetscObject)tao)->options, ((PetscObject)tao)->prefix, "-tao_fd_delta", &h, &flg));
59: PetscCall(VecDuplicate(Xin, &X));
60: PetscCall(VecCopy(Xin, X));
61: PetscCall(VecGetSize(X, &N));
62: PetscCall(VecGetOwnershipRange(X, &low, &high));
63: PetscCall(VecSetOption(X, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
64: PetscCall(VecGetArray(G, &g));
65: for (i = 0; i < N; i++) {
66: PetscCall(VecSetValue(X, i, -h, ADD_VALUES));
67: PetscCall(VecAssemblyBegin(X));
68: PetscCall(VecAssemblyEnd(X));
69: PetscCall(TaoComputeObjective(tao, X, &f));
70: PetscCall(VecSetValue(X, i, 2.0 * h, ADD_VALUES));
71: PetscCall(VecAssemblyBegin(X));
72: PetscCall(VecAssemblyEnd(X));
73: PetscCall(TaoComputeObjective(tao, X, &f2));
74: PetscCall(VecSetValue(X, i, -h, ADD_VALUES));
75: PetscCall(VecAssemblyBegin(X));
76: PetscCall(VecAssemblyEnd(X));
77: if (i >= low && i < high) g[i - low] = (f2 - f) / (2.0 * h);
78: }
79: PetscCall(VecRestoreArray(G, &g));
80: PetscCall(VecDestroy(&X));
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: /*@C
85: TaoDefaultComputeHessian - Computes the Hessian using finite differences.
87: Collective
89: Input Parameters:
90: + tao - the Tao context
91: . V - compute Hessian at this point
92: - dummy - not used
94: Output Parameters:
95: + H - Hessian matrix (not altered in this routine)
96: - B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
98: Options Database Key:
99: . -tao_fd_hessian - activates TaoDefaultComputeHessian()
101: Level: advanced
103: Notes:
104: This routine is slow and expensive, and is not optimized
105: to take advantage of sparsity in the problem. Although
106: it is not recommended for general use
107: in large-scale applications, It can be useful in checking the
108: correctness of a user-provided Hessian.
110: .seealso: `Tao`, `TaoSetHessian()`, `TaoDefaultComputeHessianColor()`, `SNESComputeJacobianDefault()`, `TaoSetGradient()`, `TaoDefaultComputeGradient()`
111: @*/
112: PetscErrorCode TaoDefaultComputeHessian(Tao tao, Vec V, Mat H, Mat B, void *dummy)
113: {
114: SNES snes;
115: DM dm;
117: PetscFunctionBegin;
118: PetscCall(PetscInfo(tao, "TAO Using finite differences w/o coloring to compute Hessian matrix\n"));
119: PetscCall(SNESCreate(PetscObjectComm((PetscObject)H), &snes));
120: PetscCall(SNESSetFunction(snes, NULL, Fsnes, tao));
121: PetscCall(SNESGetDM(snes, &dm));
122: PetscCall(DMShellSetGlobalVector(dm, V));
123: PetscCall(SNESSetUp(snes));
124: if (H) {
125: PetscInt n, N;
127: PetscCall(VecGetSize(V, &N));
128: PetscCall(VecGetLocalSize(V, &n));
129: PetscCall(MatSetSizes(H, n, n, N, N));
130: PetscCall(MatSetUp(H));
131: }
132: if (B && B != H) {
133: PetscInt n, N;
135: PetscCall(VecGetSize(V, &N));
136: PetscCall(VecGetLocalSize(V, &n));
137: PetscCall(MatSetSizes(B, n, n, N, N));
138: PetscCall(MatSetUp(B));
139: }
140: PetscCall(SNESComputeJacobianDefault(snes, V, H, B, NULL));
141: PetscCall(SNESDestroy(&snes));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: /*@C
146: TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.
148: Collective
150: Input Parameters:
151: + tao - the Tao context
152: . V - compute Hessian at this point
153: - ctx - the color object of type `MatFDColoring`
155: Output Parameters:
156: + H - Hessian matrix (not altered in this routine)
157: - B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
159: Level: advanced
161: .seealso: `Tao`, `MatColoring`, `TaoSetHessian()`, `TaoDefaultComputeHessian()`, `SNESComputeJacobianDefaultColor()`, `TaoSetGradient()`
162: @*/
163: PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H, Mat B, PetscCtx ctx)
164: {
165: MatFDColoring coloring = (MatFDColoring)ctx;
167: PetscFunctionBegin;
169: PetscCall(PetscInfo(tao, "TAO computing matrix using finite differences Hessian and coloring\n"));
170: PetscCall(MatFDColoringApply(B, coloring, V, ctx));
171: if (H != B) {
172: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
173: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
174: }
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: /*@
179: TaoDefaultComputeHessianMFFD - Computes the Hessian using finite differences with `MATMFFD`.
181: Collective
183: Input Parameters:
184: + tao - the `Tao` context
185: . X - compute Hessian at this point
186: - ctx - ignored
188: Output Parameters:
189: + H - Hessian matrix of type `MATMFFD`
190: - B - should be `NULL` or equal to `H`
192: Level: advanced
194: Note:
195: This can be passed to `TaoSetHessian()` to use `MATMFFD` for approximate Hessian-vector products. The matrix `H` can originate from
196: `MatCreateMFFD()` or from `TaoTermCreateHessianMFFD()`.
198: .seealso: `Tao`, `MATMFFD`, `MatCreateMFFD()`, `TaoTermCreateHessianMFFD()`
199: @*/
200: PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao, Vec X, Mat H, Mat B, PetscCtx ctx)
201: {
202: PetscInt n, N;
203: PetscBool assembled;
205: PetscFunctionBegin;
206: PetscCheck(!B || B == H, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "Preconditioning Hessian matrix");
207: PetscCall(MatAssembled(H, &assembled));
208: if (!assembled) {
209: PetscCall(VecGetSize(X, &N));
210: PetscCall(VecGetLocalSize(X, &n));
211: PetscCall(MatSetSizes(H, n, n, N, N));
212: PetscCall(MatSetType(H, MATMFFD));
213: PetscCall(MatSetUp(H));
214: PetscCall(MatMFFDSetFunction(H, (PetscErrorCode (*)(void *, Vec, Vec))TaoComputeGradient, tao));
215: }
216: PetscCall(MatMFFDSetBase(H, X, NULL));
217: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
218: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }