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: }