1 #include <petsctao.h> /*I "petsctao.h" I*/
2 #include <petsc/private/taoimpl.h>
3 #include <petscsnes.h>
4 #include <petscdmshell.h>
5
6 /*
7 For finited difference computations of the Hessian, we use PETSc's SNESComputeJacobianDefault
8 */
Fsnes(SNES snes,Vec X,Vec G,PetscCtx ctx)9 static PetscErrorCode Fsnes(SNES snes, Vec X, Vec G, PetscCtx ctx)
10 {
11 Tao tao = (Tao)ctx;
12
13 PetscFunctionBegin;
14 PetscValidHeaderSpecific(tao, TAO_CLASSID, 4);
15 PetscCall(TaoComputeGradient(tao, X, G));
16 PetscFunctionReturn(PETSC_SUCCESS);
17 }
18
19 /*@C
20 TaoDefaultComputeGradient - computes the gradient using finite differences.
21
22 Collective
23
24 Input Parameters:
25 + tao - the Tao context
26 . Xin - compute gradient at this point
27 - dummy - not used
28
29 Output Parameter:
30 . G - Gradient Vector
31
32 Options Database Key:
33 + -tao_fd_gradient - activates TaoDefaultComputeGradient()
34 - -tao_fd_delta <delta> - change in X used to calculate finite differences
35
36 Level: advanced
37
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
46
47 .seealso: `Tao`, `TaoSetGradient()`
48 @*/
TaoDefaultComputeGradient(Tao tao,Vec Xin,Vec G,void * dummy)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;
57
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 }
84
85 /*@C
86 TaoDefaultComputeHessian - Computes the Hessian using finite differences.
87
88 Collective
89
90 Input Parameters:
91 + tao - the Tao context
92 . V - compute Hessian at this point
93 - dummy - not used
94
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)
98
99 Options Database Key:
100 . -tao_fd_hessian - activates TaoDefaultComputeHessian()
101
102 Level: advanced
103
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.
110
111 .seealso: `Tao`, `TaoSetHessian()`, `TaoDefaultComputeHessianColor()`, `SNESComputeJacobianDefault()`, `TaoSetGradient()`, `TaoDefaultComputeGradient()`
112 @*/
TaoDefaultComputeHessian(Tao tao,Vec V,Mat H,Mat B,void * dummy)113 PetscErrorCode TaoDefaultComputeHessian(Tao tao, Vec V, Mat H, Mat B, void *dummy)
114 {
115 SNES snes;
116 DM dm;
117
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;
127
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;
135
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 }
145
146 /*@C
147 TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences.
148
149 Collective
150
151 Input Parameters:
152 + tao - the Tao context
153 . V - compute Hessian at this point
154 - ctx - the color object of type `MatFDColoring`
155
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)
159
160 Level: advanced
161
162 .seealso: `Tao`, `MatColoring`, `TaoSetHessian()`, `TaoDefaultComputeHessian()`, `SNESComputeJacobianDefaultColor()`, `TaoSetGradient()`
163 @*/
TaoDefaultComputeHessianColor(Tao tao,Vec V,Mat H,Mat B,PetscCtx ctx)164 PetscErrorCode TaoDefaultComputeHessianColor(Tao tao, Vec V, Mat H, Mat B, PetscCtx ctx)
165 {
166 MatFDColoring coloring = (MatFDColoring)ctx;
167
168 PetscFunctionBegin;
169 PetscValidHeaderSpecific(coloring, MAT_FDCOLORING_CLASSID, 5);
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 }
178
TaoDefaultComputeHessianMFFD(Tao tao,Vec X,Mat H,Mat B,PetscCtx ctx)179 PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao, Vec X, Mat H, Mat B, PetscCtx ctx)
180 {
181 PetscInt n, N;
182 PetscBool assembled;
183
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 }
200