Actual source code: taoimpl.h
1: #pragma once
3: #include <petsctao.h>
4: #include <petsctaolinesearch.h>
5: #include <petsc/private/petscimpl.h>
7: PETSC_EXTERN PetscBool TaoRegisterAllCalled;
8: PETSC_EXTERN PetscErrorCode TaoRegisterAll(void);
10: typedef struct _TaoOps *TaoOps;
12: struct _TaoOps {
13: /* Methods set by application */
14: PetscErrorCode (*computeresidual)(Tao, Vec, Vec, void *);
15: PetscErrorCode (*computeresidualjacobian)(Tao, Vec, Mat, Mat, void *);
16: PetscErrorCode (*computeconstraints)(Tao, Vec, Vec, void *);
17: PetscErrorCode (*computeinequalityconstraints)(Tao, Vec, Vec, void *);
18: PetscErrorCode (*computeequalityconstraints)(Tao, Vec, Vec, void *);
19: PetscErrorCode (*computejacobian)(Tao, Vec, Mat, Mat, void *);
20: PetscErrorCode (*computejacobianstate)(Tao, Vec, Mat, Mat, Mat, void *);
21: PetscErrorCode (*computejacobiandesign)(Tao, Vec, Mat, void *);
22: PetscErrorCode (*computejacobianinequality)(Tao, Vec, Mat, Mat, void *);
23: PetscErrorCode (*computejacobianequality)(Tao, Vec, Mat, Mat, void *);
24: PetscErrorCode (*computebounds)(Tao, Vec, Vec, void *);
25: PetscErrorCode (*update)(Tao, PetscInt, void *);
26: PetscErrorCode (*convergencetest)(Tao, void *);
27: PetscErrorCode (*convergencedestroy)(void *);
29: /* Methods set by solver */
30: PetscErrorCode (*computedual)(Tao, Vec, Vec);
31: PetscErrorCode (*setup)(Tao);
32: PetscErrorCode (*solve)(Tao);
33: PetscErrorCode (*view)(Tao, PetscViewer);
34: PetscErrorCode (*setfromoptions)(Tao, PetscOptionItems);
35: PetscErrorCode (*destroy)(Tao);
36: };
38: #define MAXTAOMONITORS 10
40: typedef struct _n_TaoTermMapping TaoTermMapping;
42: /*S
43: TaoTermMapping - Object held by either `Tao` or `TaoTerm` with `TAOTERMSUM` type
44: that contain necessary information regarding mapping matrix
46: Level: developer
48: Notes:
49: This object is internally handled in `TaoAddTerm()`, which couples `TaoTerm` with
50: appropriately mapped gradients, mapped Hessians, and other necessary information.
52: This object is also used in `TAOTERMSUM` to handle mapping of summands of `TAOTERMSUM`.
54: Users would not need to work directly with the `TaoTermMapping`, rather, they work with
55: `Tao` with added terms via `TaoAddTerm()`, or with `TaoTerm` with type `TAOTERMSUM`.
57: Developer Note:
58: Currently, as `MatPtAP` does not support diagonal matrices, internal work matrices was added
59: for a workaround.
61: .seealso: [](ch_tao), `Tao`, `TaoAddTerm()`, `TAOTERMSUM`,
62: S*/
63: struct _n_TaoTermMapping {
64: char *prefix;
65: TaoTerm term;
66: PetscReal scale;
67: Mat map;
68: Vec _map_output;
69: Vec _unmapped_gradient;
70: Vec _mapped_gradient;
71: Mat _unmapped_H;
72: Mat _unmapped_Hpre;
73: Mat _mapped_H;
74: Mat _mapped_Hpre;
75: Mat _mapped_H_work; /* Temporary work matrices for PtAP for diagonal A */
76: Mat _mapped_Hpre_work;
77: TaoTermMask mask;
78: };
80: #define TaoTermObjectiveMasked(a) ((a) & TAOTERM_MASK_OBJECTIVE)
81: #define TaoTermGradientMasked(a) ((a) & TAOTERM_MASK_GRADIENT)
82: #define TaoTermHessianMasked(a) ((a) & TAOTERM_MASK_HESSIAN)
84: struct _p_Tao {
85: PETSCHEADER(struct _TaoOps);
86: PetscCtx ctx; /* user provided context */
87: void *user_lsresP;
88: void *user_lsjacP;
89: void *user_conP;
90: void *user_con_equalityP;
91: void *user_con_inequalityP;
92: void *user_jacP;
93: void *user_jac_equalityP;
94: void *user_jac_inequalityP;
95: void *user_jac_stateP;
96: void *user_jac_designP;
97: void *user_boundsP;
98: void *user_update;
100: PetscErrorCode (*monitor[MAXTAOMONITORS])(Tao, void *);
101: PetscCtxDestroyFn *monitordestroy[MAXTAOMONITORS];
102: void *monitorcontext[MAXTAOMONITORS];
103: PetscInt numbermonitors;
104: void *cnvP;
105: TaoConvergedReason reason;
107: PetscBool setupcalled;
108: void *data;
110: DM dm; // DM used to create linear algebra objects
111: Vec solution;
112: Vec gradient;
113: Vec stepdirection;
114: Vec XL;
115: Vec XU;
116: Vec IL;
117: Vec IU;
118: Vec DI;
119: Vec DE;
120: Mat hessian;
121: Mat hessian_pre;
122: Mat gradient_norm;
123: Vec gradient_norm_tmp;
124: Vec ls_res;
125: Mat ls_jac;
126: Mat ls_jac_pre;
127: Vec res_weights_v;
128: PetscInt res_weights_n;
129: PetscInt *res_weights_rows;
130: PetscInt *res_weights_cols;
131: PetscReal *res_weights_w;
132: Vec constraints;
133: Vec constraints_equality;
134: Vec constraints_inequality;
135: Mat jacobian;
136: Mat jacobian_pre;
137: Mat jacobian_inequality;
138: Mat jacobian_inequality_pre;
139: Mat jacobian_equality;
140: Mat jacobian_equality_pre;
141: Mat jacobian_state;
142: Mat jacobian_state_inv;
143: Mat jacobian_design;
144: Mat jacobian_state_pre;
145: Mat jacobian_design_pre;
146: IS state_is;
147: IS design_is;
148: PetscReal step;
149: PetscReal residual;
150: PetscReal gnorm0;
151: PetscReal cnorm;
152: PetscReal cnorm0;
153: PetscReal fc;
155: PetscInt max_constraints;
156: PetscInt niter;
157: PetscInt ntotalits;
158: PetscInt nconstraints;
159: PetscInt niconstraints;
160: PetscInt neconstraints;
161: PetscInt nres;
162: PetscInt njac;
163: PetscInt njac_equality;
164: PetscInt njac_inequality;
165: PetscInt njac_state;
166: PetscInt njac_design;
168: PetscInt ksp_its; /* KSP iterations for this solver iteration */
169: PetscInt ksp_tot_its; /* Total (cumulative) KSP iterations */
171: TaoLineSearch linesearch;
172: PetscBool lsflag; /* goes up when line search fails */
173: KSP ksp;
174: PetscReal trust; /* Current trust region */
176: /* EW type forcing term */
177: PetscBool ksp_ewconv;
178: SNES snes_ewdummy;
180: PetscObjectParameterDeclare(PetscReal, gatol);
181: PetscObjectParameterDeclare(PetscReal, grtol);
182: PetscObjectParameterDeclare(PetscReal, gttol);
183: PetscObjectParameterDeclare(PetscReal, catol);
184: PetscObjectParameterDeclare(PetscReal, crtol);
185: PetscObjectParameterDeclare(PetscReal, steptol);
186: PetscObjectParameterDeclare(PetscReal, fmin);
187: PetscObjectParameterDeclare(PetscInt, max_it);
188: PetscObjectParameterDeclare(PetscInt, max_funcs);
189: PetscObjectParameterDeclare(PetscReal, trust0); /* initial trust region radius */
191: PetscBool printreason;
192: PetscBool viewsolution;
193: PetscBool viewgradient;
194: PetscBool viewconstraints;
195: PetscBool viewhessian;
196: PetscBool viewjacobian;
197: PetscBool bounded;
198: PetscBool constrained;
199: PetscBool eq_constrained;
200: PetscBool ineq_constrained;
201: PetscBool ineq_doublesided;
202: PetscBool header_printed;
203: PetscBool recycle;
205: TaoSubsetType subset_type;
206: PetscInt hist_max; /* Number of iteration histories to keep */
207: PetscReal *hist_obj; /* obj value at each iteration */
208: PetscReal *hist_resid; /* residual at each iteration */
209: PetscReal *hist_cnorm; /* constraint norm at each iteration */
210: PetscInt *hist_lits; /* number of ksp its at each TAO iteration */
211: PetscInt hist_len;
212: PetscBool hist_reset;
213: PetscBool hist_malloc;
215: TaoTermMapping objective_term; /* TaoTerm in use */
216: Vec objective_parameters;
217: PetscInt num_terms;
218: PetscBool term_set;
220: TaoTerm callbacks; /* TAOTERMCALLBACKS for the original callbacks */
221: PetscBool uses_hessian_matrices;
222: PetscBool uses_gradient;
223: };
225: PETSC_EXTERN PetscLogEvent TAO_Solve;
226: PETSC_EXTERN PetscLogEvent TAO_ConstraintsEval;
227: PETSC_EXTERN PetscLogEvent TAO_JacobianEval;
228: PETSC_INTERN PetscLogEvent TAO_ResidualEval;
230: PETSC_INTERN PetscLogEvent TAOTERM_ObjectiveEval;
231: PETSC_INTERN PetscLogEvent TAOTERM_GradientEval;
232: PETSC_INTERN PetscLogEvent TAOTERM_ObjGradEval;
233: PETSC_INTERN PetscLogEvent TAOTERM_HessianEval;
235: static inline PetscErrorCode TaoLogConvergenceHistory(Tao tao, PetscReal obj, PetscReal resid, PetscReal cnorm, PetscInt totits)
236: {
237: PetscFunctionBegin;
238: if (tao->hist_max > tao->hist_len) {
239: if (tao->hist_obj) tao->hist_obj[tao->hist_len] = obj;
240: if (tao->hist_resid) tao->hist_resid[tao->hist_len] = resid;
241: if (tao->hist_cnorm) tao->hist_cnorm[tao->hist_len] = cnorm;
242: if (tao->hist_lits) {
243: PetscInt sits = totits;
244: PetscCheck(tao->hist_len >= 0, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "History length cannot be negative");
245: for (PetscInt i = 0; i < tao->hist_len; i++) sits -= tao->hist_lits[i];
246: tao->hist_lits[tao->hist_len] = sits;
247: }
248: tao->hist_len++;
249: }
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: PETSC_INTERN PetscErrorCode TaoTestGradient_Internal(Tao, Vec, Vec, PetscViewer, PetscViewer);
255: typedef struct _TaoTermOps *TaoTermOps;
257: struct _TaoTermOps {
258: PetscErrorCode (*setfromoptions)(TaoTerm, PetscOptionItems);
259: PetscErrorCode (*setup)(TaoTerm);
260: PetscErrorCode (*view)(TaoTerm, PetscViewer);
261: PetscErrorCode (*destroy)(TaoTerm);
263: TaoTermObjectiveFn *objective;
264: TaoTermObjectiveAndGradientFn *objectiveandgradient;
265: TaoTermGradientFn *gradient;
266: TaoTermHessianFn *hessian;
268: PetscErrorCode (*isobjectivedefined)(TaoTerm, PetscBool *);
269: PetscErrorCode (*isgradientdefined)(TaoTerm, PetscBool *);
270: PetscErrorCode (*isobjectiveandgradientdefined)(TaoTerm, PetscBool *);
271: PetscErrorCode (*ishessiandefined)(TaoTerm, PetscBool *);
272: PetscErrorCode (*iscreatehessianmatricesdefined)(TaoTerm, PetscBool *);
273: PetscErrorCode (*iscomputehessianfdpossible)(TaoTerm, PetscBool3 *);
275: PetscErrorCode (*createsolutionvec)(TaoTerm, Vec *);
276: PetscErrorCode (*createparametersvec)(TaoTerm, Vec *);
277: PetscErrorCode (*createhessianmatrices)(TaoTerm, Mat *, Mat *);
278: };
280: struct _p_TaoTerm {
281: PETSCHEADER(struct _TaoTermOps);
282: void *data;
283: PetscBool setup_called;
284: Mat solution_factory; // dummies used to create vectors
285: Mat parameters_factory;
286: Mat parameters_factory_orig; // copy so that parameters_factory can be made a reference of solution_factory if parameter space == vector space
287: TaoTermParametersMode parameters_mode;
288: PetscBool Hpre_is_H; // Hessian mode data
289: MatType H_mattype;
290: MatType Hpre_mattype;
292: PetscInt ngrad_mffd;
293: PetscInt nobj; // actual objective callback invocations
294: PetscInt ngrad; // actual gradient callback invocations
295: PetscInt nobjgrad; // actual objective+gradient callback invocations
296: PetscInt nhess; // actual Hessian callback invocations
298: PetscReal fd_delta; // for TaoTermComputeGradientFD()
299: PetscBool fd_gradient; // use finite differences for the gradient
300: PetscBool fd_hessian; // use finite differences for the Hessian
301: };
303: PETSC_INTERN PetscErrorCode TaoTermRegisterAll(void);
305: PETSC_INTERN PetscErrorCode TaoTermCreateCallbacks(Tao, TaoTerm *);
307: PETSC_INTERN PetscErrorCode TaoTermCreate_ElementwiseDivergence_Internal(TaoTerm);
308: PETSC_INTERN PetscErrorCode TaoTermDestroy_ElementwiseDivergence_Internal(TaoTerm);
310: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetObjective(TaoTerm, PetscErrorCode (*)(Tao, Vec, PetscReal *, PetscCtx), PetscCtx);
311: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetGradient(TaoTerm, PetscErrorCode (*)(Tao, Vec, Vec, PetscCtx), PetscCtx);
312: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetObjectiveAndGradient(TaoTerm, PetscErrorCode (*)(Tao, Vec, PetscReal *, Vec, PetscCtx), PetscCtx);
313: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetHessian(TaoTerm, PetscErrorCode (*)(Tao, Vec, Mat, Mat, PetscCtx), PetscCtx);
315: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetObjective(TaoTerm, PetscErrorCode (**)(Tao, Vec, PetscReal *, PetscCtx), PetscCtxRt);
316: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetGradient(TaoTerm, PetscErrorCode (**)(Tao, Vec, Vec, PetscCtx), PetscCtxRt);
317: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetObjectiveAndGradient(TaoTerm, PetscErrorCode (**)(Tao, Vec, PetscReal *, Vec, PetscCtx), PetscCtxRt);
318: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetHessian(TaoTerm, PetscErrorCode (**)(Tao, Vec, Mat, Mat, PetscCtx), PetscCtxRt);
320: PETSC_INTERN PetscErrorCode TaoTermMappingSetData(TaoTermMapping *, const char *, PetscReal, TaoTerm, Mat);
321: PETSC_INTERN PetscErrorCode TaoTermMappingGetData(TaoTermMapping *, const char **, PetscReal *, TaoTerm *, Mat *);
322: PETSC_INTERN PetscErrorCode TaoTermMappingReset(TaoTermMapping *);
323: PETSC_INTERN PetscErrorCode TaoTermMappingComputeObjective(TaoTermMapping *, Vec, Vec, InsertMode, PetscReal *);
324: PETSC_INTERN PetscErrorCode TaoTermMappingComputeGradient(TaoTermMapping *, Vec, Vec, InsertMode, Vec);
325: PETSC_INTERN PetscErrorCode TaoTermMappingComputeObjectiveAndGradient(TaoTermMapping *, Vec, Vec, InsertMode, PetscReal *, Vec);
326: PETSC_INTERN PetscErrorCode TaoTermMappingComputeHessian(TaoTermMapping *, Vec, Vec, InsertMode, Mat, Mat);
327: PETSC_INTERN PetscErrorCode TaoTermMappingSetUp(TaoTermMapping *);
328: PETSC_INTERN PetscErrorCode TaoTermMappingCreateSolutionVec(TaoTermMapping *, Vec *);
329: PETSC_INTERN PetscErrorCode TaoTermMappingCreateParametersVec(TaoTermMapping *, Vec *);
330: PETSC_INTERN PetscErrorCode TaoTermMappingCreateHessianMatrices(TaoTermMapping *, Mat *, Mat *);
332: PETSC_INTERN PetscErrorCode VecIfNotCongruentGetSameLayoutVec(Vec, Vec *);
334: PETSC_INTERN PetscErrorCode TaoTermCreateHessianMatricesDefault_H_Internal(TaoTerm, Mat *, Mat *, PetscBool, MatType);
335: PETSC_INTERN PetscErrorCode TaoTermCreateHessianMatricesDefault_Hpre_Internal(TaoTerm, Mat *, Mat *, PetscBool, MatType);