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;
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. Use the tao method TAOTEST
44: to get an indication of whether your gradient is correct.
45: This finite difference gradient evaluation can be set using the routine `TaoSetGradient()` or by using the command line option -tao_fd_gradient
47: .seealso: `Tao`, `TaoSetGradient()`
48: @*/
49: PetscErrorCode TaoDefaultComputeGradient(Tao tao, Vec Xin, Vec G, void *dummy)
50: {
51: Vec X;
52: PetscScalar *g;
53: PetscReal f, f2;
54: PetscInt low, high, N, i;
55: PetscBool flg;
56: PetscReal h = .5 * PETSC_SQRT_MACHINE_EPSILON;
58: PetscFunctionBegin;
59: PetscCall(PetscOptionsGetReal(((PetscObject)tao)->options, ((PetscObject)tao)->prefix, "-tao_fd_delta", &h, &flg));
60: PetscCall(VecDuplicate(Xin, &X));
61: PetscCall(VecCopy(Xin, X));
62: PetscCall(VecGetSize(X, &N));
63: PetscCall(VecGetOwnershipRange(X, &low, &high));
64: PetscCall(VecSetOption(X, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
65: PetscCall(VecGetArray(G, &g));
66: for (i = 0; i < N; i++) {
67: PetscCall(VecSetValue(X, i, -h, ADD_VALUES));
68: PetscCall(VecAssemblyBegin(X));
69: PetscCall(VecAssemblyEnd(X));
70: PetscCall(TaoComputeObjective(tao, X, &f));
71: PetscCall(VecSetValue(X, i, 2.0 * h, ADD_VALUES));
72: PetscCall(VecAssemblyBegin(X));
73: PetscCall(VecAssemblyEnd(X));
74: PetscCall(TaoComputeObjective(tao, X, &f2));
75: PetscCall(VecSetValue(X, i, -h, ADD_VALUES));
76: PetscCall(VecAssemblyBegin(X));
77: PetscCall(VecAssemblyEnd(X));
78: if (i >= low && i < high) g[i - low] = (f2 - f) / (2.0 * h);
79: }
80: PetscCall(VecRestoreArray(G, &g));
81: PetscCall(VecDestroy(&X));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@C
86: TaoDefaultComputeHessian - Computes the Hessian using finite differences.
88: Collective
90: Input Parameters:
91: + tao - the Tao context
92: . V - compute Hessian at this point
93: - dummy - not used
95: Output Parameters:
96: + H - Hessian matrix (not altered in this routine)
97: - B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
99: Options Database Key:
100: . -tao_fd_hessian - activates TaoDefaultComputeHessian()
102: Level: advanced
104: Notes:
105: This routine is slow and expensive, and is not optimized
106: to take advantage of sparsity in the problem. Although
107: it is not recommended for general use
108: in large-scale applications, It can be useful in checking the
109: correctness of a user-provided Hessian.
111: .seealso: `Tao`, `TaoSetHessian()`, `TaoDefaultComputeHessianColor()`, `SNESComputeJacobianDefault()`, `TaoSetGradient()`, `TaoDefaultComputeGradient()`
112: @*/
113: PetscErrorCode TaoDefaultComputeHessian(Tao tao, Vec V, Mat H, Mat B, void *dummy)
114: {
115: SNES snes;
116: DM dm;
118: PetscFunctionBegin;
119: PetscCall(PetscInfo(tao, "TAO Using finite differences w/o coloring to compute Hessian matrix\n"));
120: PetscCall(SNESCreate(PetscObjectComm((PetscObject)H), &snes));
121: PetscCall(SNESSetFunction(snes, NULL, Fsnes, tao));
122: PetscCall(SNESGetDM(snes, &dm));
123: PetscCall(DMShellSetGlobalVector(dm, V));
124: PetscCall(SNESSetUp(snes));
125: if (H) {
126: PetscInt n, N;
128: PetscCall(VecGetSize(V, &N));
129: PetscCall(VecGetLocalSize(V, &n));
130: PetscCall(MatSetSizes(H, n, n, N, N));
131: PetscCall(MatSetUp(H));
132: }
133: if (B && B != H) {
134: PetscInt n, N;
136: PetscCall(VecGetSize(V, &N));
137: PetscCall(VecGetLocalSize(V, &n));
138: PetscCall(MatSetSizes(B, n, n, N, N));
139: PetscCall(MatSetUp(B));
140: }
141: PetscCall(SNESComputeJacobianDefault(snes, V, H, B, NULL));
142: PetscCall(SNESDestroy(&snes));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: /*@C
147: TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.
149: Collective
151: Input Parameters:
152: + tao - the Tao context
153: . V - compute Hessian at this point
154: - ctx - the color object of type `MatFDColoring`
156: Output Parameters:
157: + H - Hessian matrix (not altered in this routine)
158: - B - newly computed Hessian matrix to use with preconditioner (generally the same as H)
160: Level: advanced
162: .seealso: `Tao`, `MatColoring`, `TaoSetHessian()`, `TaoDefaultComputeHessian()`, `SNESComputeJacobianDefaultColor()`, `TaoSetGradient()`
163: @*/
164: PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H, Mat B, void *ctx)
165: {
166: MatFDColoring coloring = (MatFDColoring)ctx;
168: PetscFunctionBegin;
170: PetscCall(PetscInfo(tao, "TAO computing matrix using finite differences Hessian and coloring\n"));
171: PetscCall(MatFDColoringApply(B, coloring, V, ctx));
172: if (H != B) {
173: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
174: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
175: }
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao, Vec X, Mat H, Mat B, void *ctx)
180: {
181: PetscInt n, N;
182: PetscBool assembled;
184: PetscFunctionBegin;
185: PetscCheck(!B || B == H, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "Preconditioning Hessian matrix");
186: PetscCall(MatAssembled(H, &assembled));
187: if (!assembled) {
188: PetscCall(VecGetSize(X, &N));
189: PetscCall(VecGetLocalSize(X, &n));
190: PetscCall(MatSetSizes(H, n, n, N, N));
191: PetscCall(MatSetType(H, MATMFFD));
192: PetscCall(MatSetUp(H));
193: PetscCall(MatMFFDSetFunction(H, (PetscErrorCode (*)(void *, Vec, Vec))TaoComputeGradient, tao));
194: }
195: PetscCall(MatMFFDSetBase(H, X, NULL));
196: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
197: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }