1 #pragma once
2
3 #include <petsctao.h>
4 #include <petsctaolinesearch.h>
5 #include <petsc/private/petscimpl.h>
6
7 PETSC_EXTERN PetscBool TaoRegisterAllCalled;
8 PETSC_EXTERN PetscErrorCode TaoRegisterAll(void);
9
10 typedef struct _TaoOps *TaoOps;
11
12 struct _TaoOps {
13 /* Methods set by application */
14 PetscErrorCode (*computeobjective)(Tao, Vec, PetscReal *, void *);
15 PetscErrorCode (*computeobjectiveandgradient)(Tao, Vec, PetscReal *, Vec, void *);
16 PetscErrorCode (*computegradient)(Tao, Vec, Vec, void *);
17 PetscErrorCode (*computehessian)(Tao, Vec, Mat, Mat, void *);
18 PetscErrorCode (*computeresidual)(Tao, Vec, Vec, void *);
19 PetscErrorCode (*computeresidualjacobian)(Tao, Vec, Mat, Mat, void *);
20 PetscErrorCode (*computeconstraints)(Tao, Vec, Vec, void *);
21 PetscErrorCode (*computeinequalityconstraints)(Tao, Vec, Vec, void *);
22 PetscErrorCode (*computeequalityconstraints)(Tao, Vec, Vec, void *);
23 PetscErrorCode (*computejacobian)(Tao, Vec, Mat, Mat, void *);
24 PetscErrorCode (*computejacobianstate)(Tao, Vec, Mat, Mat, Mat, void *);
25 PetscErrorCode (*computejacobiandesign)(Tao, Vec, Mat, void *);
26 PetscErrorCode (*computejacobianinequality)(Tao, Vec, Mat, Mat, void *);
27 PetscErrorCode (*computejacobianequality)(Tao, Vec, Mat, Mat, void *);
28 PetscErrorCode (*computebounds)(Tao, Vec, Vec, void *);
29 PetscErrorCode (*update)(Tao, PetscInt, void *);
30 PetscErrorCode (*convergencetest)(Tao, void *);
31 PetscErrorCode (*convergencedestroy)(void *);
32
33 /* Methods set by solver */
34 PetscErrorCode (*computedual)(Tao, Vec, Vec);
35 PetscErrorCode (*setup)(Tao);
36 PetscErrorCode (*solve)(Tao);
37 PetscErrorCode (*view)(Tao, PetscViewer);
38 PetscErrorCode (*setfromoptions)(Tao, PetscOptionItems);
39 PetscErrorCode (*destroy)(Tao);
40 };
41
42 #define MAXTAOMONITORS 10
43
44 struct _p_Tao {
45 PETSCHEADER(struct _TaoOps);
46 PetscCtx ctx; /* user provided context */
47 void *user_objP;
48 void *user_objgradP;
49 void *user_gradP;
50 void *user_hessP;
51 void *user_lsresP;
52 void *user_lsjacP;
53 void *user_conP;
54 void *user_con_equalityP;
55 void *user_con_inequalityP;
56 void *user_jacP;
57 void *user_jac_equalityP;
58 void *user_jac_inequalityP;
59 void *user_jac_stateP;
60 void *user_jac_designP;
61 void *user_boundsP;
62 void *user_update;
63
64 PetscErrorCode (*monitor[MAXTAOMONITORS])(Tao, void *);
65 PetscCtxDestroyFn *monitordestroy[MAXTAOMONITORS];
66 void *monitorcontext[MAXTAOMONITORS];
67 PetscInt numbermonitors;
68 void *cnvP;
69 TaoConvergedReason reason;
70
71 PetscBool setupcalled;
72 void *data;
73
74 Vec solution;
75 Vec gradient;
76 Vec stepdirection;
77 Vec XL;
78 Vec XU;
79 Vec IL;
80 Vec IU;
81 Vec DI;
82 Vec DE;
83 Mat hessian;
84 Mat hessian_pre;
85 Mat gradient_norm;
86 Vec gradient_norm_tmp;
87 Vec ls_res;
88 Mat ls_jac;
89 Mat ls_jac_pre;
90 Vec res_weights_v;
91 PetscInt res_weights_n;
92 PetscInt *res_weights_rows;
93 PetscInt *res_weights_cols;
94 PetscReal *res_weights_w;
95 Vec constraints;
96 Vec constraints_equality;
97 Vec constraints_inequality;
98 Mat jacobian;
99 Mat jacobian_pre;
100 Mat jacobian_inequality;
101 Mat jacobian_inequality_pre;
102 Mat jacobian_equality;
103 Mat jacobian_equality_pre;
104 Mat jacobian_state;
105 Mat jacobian_state_inv;
106 Mat jacobian_design;
107 Mat jacobian_state_pre;
108 Mat jacobian_design_pre;
109 IS state_is;
110 IS design_is;
111 PetscReal step;
112 PetscReal residual;
113 PetscReal gnorm0;
114 PetscReal cnorm;
115 PetscReal cnorm0;
116 PetscReal fc;
117
118 PetscInt max_constraints;
119 PetscInt nfuncs;
120 PetscInt ngrads;
121 PetscInt nfuncgrads;
122 PetscInt nhess;
123 PetscInt niter;
124 PetscInt ntotalits;
125 PetscInt nconstraints;
126 PetscInt niconstraints;
127 PetscInt neconstraints;
128 PetscInt njac;
129 PetscInt njac_equality;
130 PetscInt njac_inequality;
131 PetscInt njac_state;
132 PetscInt njac_design;
133
134 PetscInt ksp_its; /* KSP iterations for this solver iteration */
135 PetscInt ksp_tot_its; /* Total (cumulative) KSP iterations */
136
137 TaoLineSearch linesearch;
138 PetscBool lsflag; /* goes up when line search fails */
139 KSP ksp;
140 PetscReal trust; /* Current trust region */
141
142 /* EW type forcing term */
143 PetscBool ksp_ewconv;
144 SNES snes_ewdummy;
145
146 PetscObjectParameterDeclare(PetscReal, gatol);
147 PetscObjectParameterDeclare(PetscReal, grtol);
148 PetscObjectParameterDeclare(PetscReal, gttol);
149 PetscObjectParameterDeclare(PetscReal, catol);
150 PetscObjectParameterDeclare(PetscReal, crtol);
151 PetscObjectParameterDeclare(PetscReal, steptol);
152 PetscObjectParameterDeclare(PetscReal, fmin);
153 PetscObjectParameterDeclare(PetscInt, max_it);
154 PetscObjectParameterDeclare(PetscInt, max_funcs);
155 PetscObjectParameterDeclare(PetscReal, trust0); /* initial trust region radius */
156
157 PetscBool printreason;
158 PetscBool viewsolution;
159 PetscBool viewgradient;
160 PetscBool viewconstraints;
161 PetscBool viewhessian;
162 PetscBool viewjacobian;
163 PetscBool bounded;
164 PetscBool constrained;
165 PetscBool eq_constrained;
166 PetscBool ineq_constrained;
167 PetscBool ineq_doublesided;
168 PetscBool header_printed;
169 PetscBool recycle;
170
171 TaoSubsetType subset_type;
172 PetscInt hist_max; /* Number of iteration histories to keep */
173 PetscReal *hist_obj; /* obj value at each iteration */
174 PetscReal *hist_resid; /* residual at each iteration */
175 PetscReal *hist_cnorm; /* constraint norm at each iteration */
176 PetscInt *hist_lits; /* number of ksp its at each TAO iteration */
177 PetscInt hist_len;
178 PetscBool hist_reset;
179 PetscBool hist_malloc;
180 };
181
182 PETSC_EXTERN PetscLogEvent TAO_Solve;
183 PETSC_EXTERN PetscLogEvent TAO_ObjectiveEval;
184 PETSC_EXTERN PetscLogEvent TAO_GradientEval;
185 PETSC_EXTERN PetscLogEvent TAO_ObjGradEval;
186 PETSC_EXTERN PetscLogEvent TAO_HessianEval;
187 PETSC_EXTERN PetscLogEvent TAO_ConstraintsEval;
188 PETSC_EXTERN PetscLogEvent TAO_JacobianEval;
189
TaoLogConvergenceHistory(Tao tao,PetscReal obj,PetscReal resid,PetscReal cnorm,PetscInt totits)190 static inline PetscErrorCode TaoLogConvergenceHistory(Tao tao, PetscReal obj, PetscReal resid, PetscReal cnorm, PetscInt totits)
191 {
192 PetscFunctionBegin;
193 if (tao->hist_max > tao->hist_len) {
194 if (tao->hist_obj) tao->hist_obj[tao->hist_len] = obj;
195 if (tao->hist_resid) tao->hist_resid[tao->hist_len] = resid;
196 if (tao->hist_cnorm) tao->hist_cnorm[tao->hist_len] = cnorm;
197 if (tao->hist_lits) {
198 PetscInt sits = totits;
199 PetscCheck(tao->hist_len >= 0, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "History length cannot be negative");
200 for (PetscInt i = 0; i < tao->hist_len; i++) sits -= tao->hist_lits[i];
201 tao->hist_lits[tao->hist_len] = sits;
202 }
203 tao->hist_len++;
204 }
205 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207