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: Vec solution;
111: Vec gradient;
112: Vec stepdirection;
113: Vec XL;
114: Vec XU;
115: Vec IL;
116: Vec IU;
117: Vec DI;
118: Vec DE;
119: Mat hessian;
120: Mat hessian_pre;
121: Mat gradient_norm;
122: Vec gradient_norm_tmp;
123: Vec ls_res;
124: Mat ls_jac;
125: Mat ls_jac_pre;
126: Vec res_weights_v;
127: PetscInt res_weights_n;
128: PetscInt *res_weights_rows;
129: PetscInt *res_weights_cols;
130: PetscReal *res_weights_w;
131: Vec constraints;
132: Vec constraints_equality;
133: Vec constraints_inequality;
134: Mat jacobian;
135: Mat jacobian_pre;
136: Mat jacobian_inequality;
137: Mat jacobian_inequality_pre;
138: Mat jacobian_equality;
139: Mat jacobian_equality_pre;
140: Mat jacobian_state;
141: Mat jacobian_state_inv;
142: Mat jacobian_design;
143: Mat jacobian_state_pre;
144: Mat jacobian_design_pre;
145: IS state_is;
146: IS design_is;
147: PetscReal step;
148: PetscReal residual;
149: PetscReal gnorm0;
150: PetscReal cnorm;
151: PetscReal cnorm0;
152: PetscReal fc;
154: PetscInt max_constraints;
155: PetscInt niter;
156: PetscInt ntotalits;
157: PetscInt nconstraints;
158: PetscInt niconstraints;
159: PetscInt neconstraints;
160: PetscInt nres;
161: PetscInt njac;
162: PetscInt njac_equality;
163: PetscInt njac_inequality;
164: PetscInt njac_state;
165: PetscInt njac_design;
167: PetscInt ksp_its; /* KSP iterations for this solver iteration */
168: PetscInt ksp_tot_its; /* Total (cumulative) KSP iterations */
170: TaoLineSearch linesearch;
171: PetscBool lsflag; /* goes up when line search fails */
172: KSP ksp;
173: PetscReal trust; /* Current trust region */
175: /* EW type forcing term */
176: PetscBool ksp_ewconv;
177: SNES snes_ewdummy;
179: PetscObjectParameterDeclare(PetscReal, gatol);
180: PetscObjectParameterDeclare(PetscReal, grtol);
181: PetscObjectParameterDeclare(PetscReal, gttol);
182: PetscObjectParameterDeclare(PetscReal, catol);
183: PetscObjectParameterDeclare(PetscReal, crtol);
184: PetscObjectParameterDeclare(PetscReal, steptol);
185: PetscObjectParameterDeclare(PetscReal, fmin);
186: PetscObjectParameterDeclare(PetscInt, max_it);
187: PetscObjectParameterDeclare(PetscInt, max_funcs);
188: PetscObjectParameterDeclare(PetscReal, trust0); /* initial trust region radius */
190: PetscBool printreason;
191: PetscBool viewsolution;
192: PetscBool viewgradient;
193: PetscBool viewconstraints;
194: PetscBool viewhessian;
195: PetscBool viewjacobian;
196: PetscBool bounded;
197: PetscBool constrained;
198: PetscBool eq_constrained;
199: PetscBool ineq_constrained;
200: PetscBool ineq_doublesided;
201: PetscBool header_printed;
202: PetscBool recycle;
204: TaoSubsetType subset_type;
205: PetscInt hist_max; /* Number of iteration histories to keep */
206: PetscReal *hist_obj; /* obj value at each iteration */
207: PetscReal *hist_resid; /* residual at each iteration */
208: PetscReal *hist_cnorm; /* constraint norm at each iteration */
209: PetscInt *hist_lits; /* number of ksp its at each TAO iteration */
210: PetscInt hist_len;
211: PetscBool hist_reset;
212: PetscBool hist_malloc;
214: TaoTermMapping objective_term; /* TaoTerm in use */
215: Vec objective_parameters;
216: PetscInt num_terms;
217: PetscBool term_set;
219: TaoTerm callbacks; /* TAOTERMCALLBACKS for the original callbacks */
220: PetscBool uses_hessian_matrices;
221: PetscBool uses_gradient;
222: };
224: PETSC_EXTERN PetscLogEvent TAO_Solve;
225: PETSC_EXTERN PetscLogEvent TAO_ConstraintsEval;
226: PETSC_EXTERN PetscLogEvent TAO_JacobianEval;
227: PETSC_INTERN PetscLogEvent TAO_ResidualEval;
229: PETSC_INTERN PetscLogEvent TAOTERM_ObjectiveEval;
230: PETSC_INTERN PetscLogEvent TAOTERM_GradientEval;
231: PETSC_INTERN PetscLogEvent TAOTERM_ObjGradEval;
232: PETSC_INTERN PetscLogEvent TAOTERM_HessianEval;
234: static inline PetscErrorCode TaoLogConvergenceHistory(Tao tao, PetscReal obj, PetscReal resid, PetscReal cnorm, PetscInt totits)
235: {
236: PetscFunctionBegin;
237: if (tao->hist_max > tao->hist_len) {
238: if (tao->hist_obj) tao->hist_obj[tao->hist_len] = obj;
239: if (tao->hist_resid) tao->hist_resid[tao->hist_len] = resid;
240: if (tao->hist_cnorm) tao->hist_cnorm[tao->hist_len] = cnorm;
241: if (tao->hist_lits) {
242: PetscInt sits = totits;
243: PetscCheck(tao->hist_len >= 0, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "History length cannot be negative");
244: for (PetscInt i = 0; i < tao->hist_len; i++) sits -= tao->hist_lits[i];
245: tao->hist_lits[tao->hist_len] = sits;
246: }
247: tao->hist_len++;
248: }
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: PETSC_INTERN PetscErrorCode TaoTestGradient_Internal(Tao, Vec, Vec, PetscViewer, PetscViewer);
254: typedef struct _TaoTermOps *TaoTermOps;
256: struct _TaoTermOps {
257: PetscErrorCode (*setfromoptions)(TaoTerm, PetscOptionItems);
258: PetscErrorCode (*setup)(TaoTerm);
259: PetscErrorCode (*view)(TaoTerm, PetscViewer);
260: PetscErrorCode (*destroy)(TaoTerm);
262: TaoTermObjectiveFn *objective;
263: TaoTermObjectiveAndGradientFn *objectiveandgradient;
264: TaoTermGradientFn *gradient;
265: TaoTermHessianFn *hessian;
267: PetscErrorCode (*isobjectivedefined)(TaoTerm, PetscBool *);
268: PetscErrorCode (*isgradientdefined)(TaoTerm, PetscBool *);
269: PetscErrorCode (*isobjectiveandgradientdefined)(TaoTerm, PetscBool *);
270: PetscErrorCode (*ishessiandefined)(TaoTerm, PetscBool *);
271: PetscErrorCode (*iscreatehessianmatricesdefined)(TaoTerm, PetscBool *);
272: PetscErrorCode (*iscomputehessianfdpossible)(TaoTerm, PetscBool3 *);
274: PetscErrorCode (*createsolutionvec)(TaoTerm, Vec *);
275: PetscErrorCode (*createparametersvec)(TaoTerm, Vec *);
276: PetscErrorCode (*createhessianmatrices)(TaoTerm, Mat *, Mat *);
277: };
279: struct _p_TaoTerm {
280: PETSCHEADER(struct _TaoTermOps);
281: void *data;
282: PetscBool setup_called;
283: Mat solution_factory; // dummies used to create vectors
284: Mat parameters_factory;
285: Mat parameters_factory_orig; // copy so that parameters_factory can be made a reference of solution_factory if parameter space == vector space
286: TaoTermParametersMode parameters_mode;
287: PetscBool Hpre_is_H; // Hessian mode data
288: MatType H_mattype;
289: MatType Hpre_mattype;
291: PetscInt ngrad_mffd;
292: PetscInt nobj; // actual objective callback invocations
293: PetscInt ngrad; // actual gradient callback invocations
294: PetscInt nobjgrad; // actual objective+gradient callback invocations
295: PetscInt nhess; // actual Hessian callback invocations
297: PetscReal fd_delta; // for TaoTermComputeGradientFD()
298: PetscBool fd_gradient; // use finite differences for the gradient
299: PetscBool fd_hessian; // use finite differences for the Hessian
300: };
302: PETSC_INTERN PetscErrorCode TaoTermRegisterAll(void);
304: PETSC_INTERN PetscErrorCode TaoTermCreateCallbacks(Tao, TaoTerm *);
306: PETSC_INTERN PetscErrorCode TaoTermCreate_ElementwiseDivergence_Internal(TaoTerm);
307: PETSC_INTERN PetscErrorCode TaoTermDestroy_ElementwiseDivergence_Internal(TaoTerm);
309: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetObjective(TaoTerm, PetscErrorCode (*)(Tao, Vec, PetscReal *, PetscCtx), PetscCtx);
310: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetGradient(TaoTerm, PetscErrorCode (*)(Tao, Vec, Vec, PetscCtx), PetscCtx);
311: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetObjectiveAndGradient(TaoTerm, PetscErrorCode (*)(Tao, Vec, PetscReal *, Vec, PetscCtx), PetscCtx);
312: PETSC_INTERN PetscErrorCode TaoTermCallbacksSetHessian(TaoTerm, PetscErrorCode (*)(Tao, Vec, Mat, Mat, PetscCtx), PetscCtx);
314: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetObjective(TaoTerm, PetscErrorCode (**)(Tao, Vec, PetscReal *, PetscCtx), PetscCtxRt);
315: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetGradient(TaoTerm, PetscErrorCode (**)(Tao, Vec, Vec, PetscCtx), PetscCtxRt);
316: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetObjectiveAndGradient(TaoTerm, PetscErrorCode (**)(Tao, Vec, PetscReal *, Vec, PetscCtx), PetscCtxRt);
317: PETSC_INTERN PetscErrorCode TaoTermCallbacksGetHessian(TaoTerm, PetscErrorCode (**)(Tao, Vec, Mat, Mat, PetscCtx), PetscCtxRt);
319: PETSC_INTERN PetscErrorCode TaoTermMappingSetData(TaoTermMapping *, const char *, PetscReal, TaoTerm, Mat);
320: PETSC_INTERN PetscErrorCode TaoTermMappingGetData(TaoTermMapping *, const char **, PetscReal *, TaoTerm *, Mat *);
321: PETSC_INTERN PetscErrorCode TaoTermMappingReset(TaoTermMapping *);
322: PETSC_INTERN PetscErrorCode TaoTermMappingComputeObjective(TaoTermMapping *, Vec, Vec, InsertMode, PetscReal *);
323: PETSC_INTERN PetscErrorCode TaoTermMappingComputeGradient(TaoTermMapping *, Vec, Vec, InsertMode, Vec);
324: PETSC_INTERN PetscErrorCode TaoTermMappingComputeObjectiveAndGradient(TaoTermMapping *, Vec, Vec, InsertMode, PetscReal *, Vec);
325: PETSC_INTERN PetscErrorCode TaoTermMappingComputeHessian(TaoTermMapping *, Vec, Vec, InsertMode, Mat, Mat);
326: PETSC_INTERN PetscErrorCode TaoTermMappingSetUp(TaoTermMapping *);
327: PETSC_INTERN PetscErrorCode TaoTermMappingCreateSolutionVec(TaoTermMapping *, Vec *);
328: PETSC_INTERN PetscErrorCode TaoTermMappingCreateParametersVec(TaoTermMapping *, Vec *);
329: PETSC_INTERN PetscErrorCode TaoTermMappingCreateHessianMatrices(TaoTermMapping *, Mat *, Mat *);
331: PETSC_INTERN PetscErrorCode VecIfNotCongruentGetSameLayoutVec(Vec, Vec *);
333: PETSC_INTERN PetscErrorCode TaoTermCreateHessianMatricesDefault_H_Internal(TaoTerm, Mat *, Mat *, PetscBool, MatType);
334: PETSC_INTERN PetscErrorCode TaoTermCreateHessianMatricesDefault_Hpre_Internal(TaoTerm, Mat *, Mat *, PetscBool, MatType);