Actual source code: taotermfdiff.c
1: #include <petsc/private/taoimpl.h>
2: #include <petscsnes.h>
3: #include <petscdmshell.h>
5: typedef struct _n_TaoTermWithParameters {
6: TaoTerm term; // weak-reference
7: Vec params;
8: } TaoTermWithParameters;
10: PETSC_INTERN PetscErrorCode TaoTermWithParametersDestroy(PetscCtxRt ctx)
11: {
12: TaoTermWithParameters *t = (TaoTermWithParameters *)*(void **)ctx;
14: PetscFunctionBegin;
15: PetscCall(VecDestroy(&t->params));
16: PetscCall(PetscFree(t));
17: PetscFunctionReturn(PETSC_SUCCESS);
18: }
20: static PetscErrorCode SNESFunction_TaoTerm(SNES snes, Vec X, Vec G, void *ctx)
21: {
22: TaoTermWithParameters *t = (TaoTermWithParameters *)ctx;
24: PetscFunctionBegin;
26: PetscCall(TaoTermComputeGradient(t->term, X, t->params, G));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: /*@
31: TaoTermComputeGradientFD - Approximate the gradient of a `TaoTerm` using finite differences
33: Collective
35: Input Parameters:
36: + term - a `TaoTerm`
37: . x - a solution vector
38: - params - parameters vector (may be `NULL`, see `TaoTermParametersMode`)
40: Output Parameter:
41: . g - the computed finite difference approximation to the gradient
43: Options Database Keys:
44: + -tao_term_fd_delta <delta> - change in `x` used to calculate finite differences
45: - -tao_term_gradient_use_fd <bool> - Use `TaoTermComputeGradientFD()` in `TaoTermComputeGradient()`
47: Level: advanced
49: Notes:
50: This routine is slow and expensive, and is not optimized to take advantage of
51: sparsity in the problem. Although not recommended for general use in
52: large-scale applications, it can be useful in checking the correctness of a
53: user-provided gradient. Call `TaoTermComputeGradientSetUseFD()` to start using
54: this routine in `TaoTermComputeGradient()`.
56: .seealso: [](sec_tao_term),
57: `TaoTerm`,
58: `TaoTermGetFDDelta()`,
59: `TaoTermSetFDDelta()`,
60: `TaoTermComputeGradientSetUseFD()`,
61: `TaoTermComputeGradientGetUseFD()`,
62: `TaoTermComputeHessianFD()`
63: @*/
64: PetscErrorCode TaoTermComputeGradientFD(TaoTerm term, Vec x, Vec params, Vec g)
65: {
66: Vec x_perturbed;
67: PetscScalar *_g;
68: PetscReal f, f2;
69: PetscInt low, high, N, i;
70: PetscReal h;
72: PetscFunctionBegin;
77: h = term->fd_delta;
79: PetscCall(VecDuplicate(x, &x_perturbed));
80: PetscCall(VecCopy(x, x_perturbed));
81: PetscCall(VecGetSize(x_perturbed, &N));
82: PetscCall(VecGetOwnershipRange(x_perturbed, &low, &high));
83: PetscCall(VecSetOption(x_perturbed, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
84: PetscCall(VecGetArray(g, &_g));
85: for (i = 0; i < N; i++) {
86: PetscCall(VecSetValue(x_perturbed, i, -h, ADD_VALUES));
87: PetscCall(VecAssemblyBegin(x_perturbed));
88: PetscCall(VecAssemblyEnd(x_perturbed));
89: PetscCall(TaoTermComputeObjective(term, x_perturbed, params, &f));
90: PetscCall(VecSetValue(x_perturbed, i, 2.0 * h, ADD_VALUES));
91: PetscCall(VecAssemblyBegin(x_perturbed));
92: PetscCall(VecAssemblyEnd(x_perturbed));
93: PetscCall(TaoTermComputeObjective(term, x_perturbed, params, &f2));
94: PetscCall(VecSetValue(x_perturbed, i, -h, ADD_VALUES));
95: PetscCall(VecAssemblyBegin(x_perturbed));
96: PetscCall(VecAssemblyEnd(x_perturbed));
97: if (i >= low && i < high) _g[i - low] = (f2 - f) / (2.0 * h);
98: }
99: PetscCall(VecRestoreArray(g, &_g));
100: PetscCall(VecDestroy(&x_perturbed));
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /*@
105: TaoTermComputeHessianFD - Use finite difference to compute Hessian matrix.
107: Collective
109: Input Parameters:
110: + term - a `TaoTerm`
111: . x - a solution vector
112: - params - parameters vector (may be `NULL`, see `TaoTermParametersMode`)
114: Output Parameters:
115: + H - (optional) Hessian matrix
116: - Hpre - (optional) Hessian preconditioning matrix
118: Options Database Keys:
119: + -tao_term_fd_delta <delta> - change in X used to calculate finite differences
120: - -tao_term_hessian_use_fd <bool> - Use `TaoTermComputeHessianFD()` in `TaoTermComputeHessian()`
122: Level: advanced
124: Notes:
125: This routine is slow and expensive, and is not optimized to take advantage of
126: sparsity in the problem. Although not recommended for general use in
127: large-scale applications, it can be useful in checking the correctness of a
128: user-provided Hessian. Call `TaoTermComputeHessianSetUseFD()` to start using
129: this routine in `TaoTermComputeHessian()`.
131: .seealso: [](sec_tao_term),
132: `TaoTerm`,
133: `TaoTermComputeHessian()`,
134: `TaoTermGetFDDelta()`,
135: `TaoTermSetFDDelta()`,
136: `TaoTermComputeHessianSetUseFD()`,
137: `TaoTermComputeHessianGetUseFD()`
138: @*/
139: PetscErrorCode TaoTermComputeHessianFD(TaoTerm term, Vec x, Vec params, Mat H, Mat Hpre)
140: {
141: SNES snes;
142: DM dm;
143: TaoTermWithParameters t;
145: PetscFunctionBegin;
151: PetscCheck(H || Hpre, PetscObjectComm((PetscObject)term), PETSC_ERR_ARG_NULL, "At least one of H or Hpre must be non-NULL");
152: // Note: Request for FD takes higher precedence over MATMFFD. Ignore MFFD
153: PetscCall(PetscInfo(term, "%s: TaoTerm using finite differences w/o coloring to compute Hessian matrix.\n", ((PetscObject)term)->prefix));
154: // Note: same routine as in fdiff.c
155: PetscCall(SNESCreate(PetscObjectComm((PetscObject)term), &snes));
157: t.term = term;
158: t.params = params;
159: PetscCall(SNESSetFunction(snes, NULL, SNESFunction_TaoTerm, (void *)&t));
160: PetscCall(SNESGetDM(snes, &dm));
161: PetscCall(DMShellSetGlobalVector(dm, x));
162: PetscCall(SNESSetUp(snes));
163: if (H) {
164: PetscInt n, N;
166: PetscCall(VecGetSize(x, &N));
167: PetscCall(VecGetLocalSize(x, &n));
168: PetscCall(MatSetSizes(H, n, n, N, N));
169: PetscCall(MatSetUp(H));
170: }
171: if (Hpre && Hpre != H) {
172: PetscInt n, N;
174: PetscCall(VecGetSize(x, &N));
175: PetscCall(VecGetLocalSize(x, &n));
176: PetscCall(MatSetSizes(Hpre, n, n, N, N));
177: PetscCall(MatSetUp(Hpre));
178: }
179: PetscCall(SNESComputeJacobianDefault(snes, x, H ? H : Hpre, Hpre ? Hpre : H, NULL));
180: PetscCall(SNESDestroy(&snes));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode MatMFFDFunction_TaoTermHessianShell(void *ctx, Vec x, Vec g)
185: {
186: TaoTermWithParameters *tp = (TaoTermWithParameters *)ctx;
188: PetscFunctionBegin;
189: // we expect the solution to move around in a finite difference method, but not the parameters
190: // TODO but not checking for it now
191: PetscCall(TaoTermComputeGradient(tp->term, x, tp->params, g));
192: tp->term->ngrad_mffd++;
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: static PetscErrorCode TaoTermInitializeHessianMFFD(TaoTerm term, Mat mffd)
197: {
198: TaoTermWithParameters *tp;
199: PetscLayout sol_layout;
200: VecType sol_vec_type;
201: PetscContainer container;
203: PetscFunctionBegin;
204: PetscCall(TaoTermGetSolutionLayout(term, &sol_layout));
205: PetscCall(MatSetLayouts(mffd, sol_layout, sol_layout));
206: PetscCall(TaoTermGetSolutionVecType(term, &sol_vec_type));
207: PetscCall(MatSetVecType(mffd, sol_vec_type));
208: PetscCall(MatSetType(mffd, MATMFFD));
209: PetscCall(PetscNew(&tp));
210: tp->term = term;
211: PetscCall(MatMFFDSetFunction(mffd, MatMFFDFunction_TaoTermHessianShell, tp));
212: PetscCall(MatSetOption(mffd, MAT_SYMMETRIC, PETSC_TRUE));
213: PetscCall(MatSetOption(mffd, MAT_SYMMETRY_ETERNAL, PETSC_TRUE));
214: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)term), &container));
215: PetscCall(PetscContainerSetPointer(container, (void *)tp));
216: PetscCall(PetscContainerSetCtxDestroy(container, TaoTermWithParametersDestroy));
217: PetscCall(PetscObjectCompose((PetscObject)mffd, "__TaoTermWithParameters", (PetscObject)container));
218: PetscCall(PetscContainerDestroy(&container));
219: PetscCall(PetscFree(term->Hpre_mattype));
220: PetscCall(PetscFree(term->H_mattype));
221: PetscCall(PetscStrallocpy(MATMFFD, (char **)&term->H_mattype));
222: PetscCall(PetscStrallocpy(MATMFFD, (char **)&term->Hpre_mattype));
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: /*@
227: TaoTermCreateHessianMFFD - Create a `MATMFFD` for a matrix-free finite-difference approximation of the Hessian of a `TaoTerm`
229: Collective
231: Input Parameter:
232: . term - a `TaoTerm`
234: Output Parameter:
235: . mffd - a `Mat` of type `MATMFFD`
237: Level: advanced
239: .seealso: [](sec_tao_term), `TaoTerm`, `TaoTermComputeHessianFD()`
240: @*/
241: PetscErrorCode TaoTermCreateHessianMFFD(TaoTerm term, Mat *mffd)
242: {
243: PetscFunctionBegin;
245: PetscAssertPointer(mffd, 2);
246: PetscCall(MatCreate(PetscObjectComm((PetscObject)term), mffd));
247: PetscCall(TaoTermInitializeHessianMFFD(term, *mffd));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: PetscErrorCode TaoTermComputeHessianMFFD(TaoTerm term, Vec x, Vec params, Mat H, Mat B)
252: {
253: PetscContainer container;
254: TaoTermWithParameters *tp;
256: PetscFunctionBegin;
262: PetscCall(PetscObjectQuery((PetscObject)H, "__TaoTermWithParameters", (PetscObject *)&container));
263: if (!container) {
264: PetscCall(TaoTermInitializeHessianMFFD(term, H));
265: PetscCall(PetscObjectQuery((PetscObject)H, "__TaoTermWithParameters", (PetscObject *)&container));
266: PetscCheck(container, PetscObjectComm((PetscObject)term), PETSC_ERR_PLIB, "failed to initialize mffd matrix");
267: }
268: PetscCall(PetscContainerGetPointer(container, (void **)&tp));
269: PetscCheck(tp->term == term, PetscObjectComm((PetscObject)term), PETSC_ERR_ARG_INCOMP, "Hessian shell matrix does not come from this TaoTerm");
270: if (params) PetscCall(PetscObjectReference((PetscObject)params));
271: PetscCall(VecDestroy(&tp->params));
272: tp->params = params;
273: PetscCall(MatMFFDSetBase(H, x, NULL));
274: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
275: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }