tTaoUtil.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tTaoUtil.hh (17530B)
---
1 // Copyright (C) 2012, 2013, 2014, 2015, 2017 David Maxwell and Constantine Khroulev
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #ifndef TAOUTIL_HH_W42GJNRO
20 #define TAOUTIL_HH_W42GJNRO
21
22 #include <petsc.h>
23 #include <stdexcept>
24 #include <string>
25
26 #include "pism/util/pism_utilities.hh"
27 #include "pism/util/TerminationReason.hh"
28 #include "pism/util/petscwrappers/Tao.hh"
29 #include "pism/util/error_handling.hh"
30
31 namespace pism {
32
33 //! TAO (inverse modeling) utilities.
34 namespace taoutil {
35
36 //! Encapsulate TAO's TaoSolverTerminationReason codes as a PISM TerminationReason.
37 class TAOTerminationReason: public TerminationReason {
38 public:
39 TAOTerminationReason(TaoConvergedReason r);
40 virtual void get_description(std::ostream &desc,int indent_level=0);
41 };
42
43 //! \brief An interface for solving an optimization problem with TAO where the
44 //! problem itself is defined by a separate Problem class.
45 /*! The primary interface to a TAO optimization problem is mediated by a PETSc-style
46 `TaoSolver` object. The PISM TaoBasicSolver C++ class wraps a `TaoSolver` and some of
47 its initialization boilierplate, and allows a separate class to define the function to be minimized.
48
49 To use a TaoBasicSolver you create a `Problem` class that defines the objective function and initial
50 guess, as well any auxilliary callbacks desired. The Problem class must define a
51
52 \code
53 void Problem::connect(TaoSolver solver);
54 \endcode
55
56 method which gives the `Problem` an opportunity to register its methods as callbacks to the solver,
57 perhaps taking advantage of the various `TaoFooCallback` classes provided in TaoUtil.hh to facilitate this.
58 For example, a problem class MyProblem that did nothing more than register a combined objective/gradient
59 callback could define
60
61 \code
62 void MyProblem::connect(TaoSolver tao) {
63 typedef TaoObjGradCallback<Problem,&MyProblem::evaluateObjectiveAndGradient> ObjGradCallback;
64 ObjGradCallback::connect(tao,*this);
65 }
66 \endcode
67
68 In addition to the `connect` method, a `Problem` must define
69 \code
70 TerminationReason::Ptr MyProblem::formInitialGuess(Vec *v)
71 \endcode
72 which allows the problem to set the initial guess for optimization. If the minimization
73 is successful, the solution will be found in the same vector that was returned by this method.
74
75 Assuming a `MyProblem` called `problem` has been constructed, solution
76 of the minimization is done using, for example, the TAO algorithm
77 tao_cg:
78
79 \code
80 TaoBasicSolver<MyProblem> solver(com,"tao_cg",problem);
81
82 TerminationReason::Ptr reason = solver.solve();
83
84 if (reason->succeeded()) {
85 printf("Success: %s\n",reason->description().c_str());
86 } else {
87 printf("Failure: %s\n",reason->description().c_str());
88 }
89 \endcode
90 */
91 template<class Problem>
92 class TaoBasicSolver {
93 public:
94
95 //! Construct a solver to solve `prob` using TAO algorithm `tao_type`.
96 TaoBasicSolver(MPI_Comm comm, const std::string & tao_type, Problem &prob)
97 : m_comm(comm), m_problem(prob) {
98 PetscErrorCode ierr;
99
100 ierr = TaoCreate(m_comm, m_tao.rawptr());
101 PISM_CHK(ierr, "TaoCreate");
102
103 ierr = TaoSetType(m_tao,tao_type.c_str());
104 PISM_CHK(ierr, "TaoSetType");
105
106 m_problem.connect(m_tao);
107
108 ierr = TaoSetFromOptions(m_tao);
109 PISM_CHK(ierr, "TaoSetFromOptions");
110 }
111
112 virtual ~TaoBasicSolver() {
113 // empty
114 }
115
116 //! Solve the minimization problem.
117 virtual TerminationReason::Ptr solve() {
118 PetscErrorCode ierr;
119
120 /* Solve the application */
121 Vec x0;
122 TerminationReason::Ptr reason = m_problem.formInitialGuess(&x0);
123 if (reason->failed()) {
124 TerminationReason::Ptr root_cause = reason;
125 reason.reset(new GenericTerminationReason(-1, "Unable to form initial guess"));
126 reason->set_root_cause(root_cause);
127 return reason;
128 }
129 ierr = TaoSetInitialVector(m_tao, x0);
130 PISM_CHK(ierr, "TaoSetInitialVector");
131
132 ierr = TaoSolve(m_tao);
133 PISM_CHK(ierr, "TaoSolve");
134
135 TaoConvergedReason tao_reason;
136 ierr = TaoGetConvergedReason(m_tao, &tao_reason);
137 PISM_CHK(ierr, "TaoGetConvergedReason");
138
139 reason.reset(new TAOTerminationReason(tao_reason));
140 return reason;
141 }
142
143 virtual void setMaximumIterations(int max_it) {
144 PetscErrorCode ierr = TaoSetMaximumIterations(m_tao, max_it);
145 PISM_CHK(ierr, "TaoSetMaximumIterations");
146 }
147
148 virtual Problem &problem() {
149 return m_problem;
150 }
151
152 protected:
153 MPI_Comm m_comm;
154 petsc::Tao m_tao;
155 Problem &m_problem;
156 };
157
158
159 //! \brief Adaptor to connect a TAO objective function callback to a C++ object method.
160 /*! The TAO library interfaces with user code via C-style callback functions.
161 This class makes it convenient to associate a TAO Objective callback
162 with a C++ object method. To assign
163 \code
164 void MyObject::evaluateObjective(TaoSolver tao,Vec x, double *value);
165 \endcode
166
167 as the objective function to a `TaoSolver` `tao`,
168
169 \code
170 MyObject obj;
171 TaoObjectiveCallback<MyObject>::connect(tao,obj);
172 \endcode
173
174 The method name `evaluateObjective` for the callback is hard-coded.
175 See TaoObjGradCallback for a technique to allow
176 the method name to be specified (at the expense of a little more cumbersome code).
177 */
178 template<class Problem>
179 class TaoObjectiveCallback {
180 public:
181
182 static void connect(Tao tao, Problem &p) {
183 PetscErrorCode ierr;
184 ierr = TaoSetObjectiveRoutine(tao,
185 TaoObjectiveCallback<Problem>::callback,
186 &p);
187 PISM_CHK(ierr, "TaoSetObjectiveRoutine");
188 }
189
190 protected:
191 static PetscErrorCode callback(Tao tao, Vec x, double *value, void *ctx) {
192 try {
193 Problem *p = reinterpret_cast<Problem *>(ctx);
194 PetscErrorCode ierr = p->evaluateObjective(tao,x,value); CHKERRQ(ierr);
195 } catch (...) {
196 MPI_Comm com = MPI_COMM_SELF;
197 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
198 handle_fatal_errors(com);
199 SETERRQ(com, 1, "A PISM callback failed");
200 }
201 return 0;
202 }
203 };
204
205
206 //! \brief Adaptor to connect a TAO monitoring callback to a C++ object method.
207 /*! The TAO library interfaces with user code via C-style callback functions.
208 This class makes it convenient to associate a TAO Monitor callback
209 with a C++ object method. To assign
210 \code
211 void MyObject::monitorTao(Tao tao)
212 \endcode
213
214 as the objective function to a `Tao` `tao`,
215
216 \code
217 MyObject obj;
218 TaoMonitorCallback<MyObject>::connect(tao,obj);
219 \endcode
220
221 The method name `monitorTao` for the callback is hard-coded.
222 See TaoObjGradCallback for a technique to allow
223 the method name to be specified (at the expense of a little more cumbersome code).
224 */
225 template<class Problem>
226 class TaoMonitorCallback {
227 public:
228 static void connect(Tao tao, Problem &p) {
229 PetscErrorCode ierr = TaoSetMonitor(tao,
230 TaoMonitorCallback<Problem>::callback,
231 &p, NULL);
232 PISM_CHK(ierr, "TaoSetMonitor");
233 }
234 protected:
235 static PetscErrorCode callback(Tao tao, void *ctx) {
236 try {
237 Problem *p = reinterpret_cast<Problem *>(ctx);
238 p->monitorTao(tao);
239 } catch (...) {
240 MPI_Comm com = MPI_COMM_SELF;
241 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
242 handle_fatal_errors(com);
243 SETERRQ(com, 1, "A PISM callback failed");
244 }
245 return 0;
246 }
247 };
248
249 //! \brief Adaptor to connect a TAO objective function callback to a C++ object method.
250 /*! The TAO library interfaces with user code via C-style callback functions.
251 This class makes it convenient to associate a TAO VariableBounds callback
252 with a C++ object method. To assign
253 \code
254 void MyObject::getVariableBounds(Tao tao,Vec lo, Vec hi);
255 \endcode
256
257 as the objective function to a `Tao` `tao`,
258
259 \code
260 MyObject obj;
261 TaoGetVariableBoundsCallback<MyObject>::connect(tao,obj);
262 \endcode
263
264 The method name `getVariableBounds` for the callback is hard-coded.
265 See TaoObjGradCallback for a technique to allow
266 the method name to be specified (at the expense of a little more cumbersome code).
267 */
268 template<class Problem>
269 class TaoGetVariableBoundsCallback {
270 public:
271 static void connect(Tao tao, Problem &p) {
272 PetscErrorCode ierr;
273 ierr = TaoSetVariableBoundsRoutine(tao,
274 TaoGetVariableBoundsCallback<Problem>::callback,
275 &p);
276 PISM_CHK(ierr, "TaoSetVariableBoundsRoutine");
277 }
278
279 protected:
280 static PetscErrorCode callback(Tao tao, Vec lo, Vec hi, void *ctx) {
281 try {
282 Problem *p = reinterpret_cast<Problem *>(ctx);
283 p->getVariableBounds(tao,lo,hi);
284 } catch (...) {
285 MPI_Comm com = MPI_COMM_SELF;
286 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
287 handle_fatal_errors(com);
288 SETERRQ(com, 1, "A PISM callback failed");
289 }
290 return 0;
291 }
292 };
293
294 //! \brief Adaptor to connect a TAO objective gradient callback to a C++ object method.
295 /*! The TAO library interfaces with user code via C-style callback functions.
296 This class makes it convenient to associate a TAO Objective Gradient callback
297 with a C++ object method. To assign
298 \code
299 void MyObject::evaluateGradient(Tao tao,Vec x, Vec gradient);
300 \endcode
301
302 as the objective function to a `Tao` `tao`,
303
304 \code
305 MyObject obj;
306 TaoGradientCallback<MyObject>::connect(tao,obj);
307 \endcode
308
309 The method name `evaluateGradient` for the callback is hard-coded.
310 See TaoObjGradCallback for a technique to allow
311 the method name to be specified (at the expense of a little more cumbersome code).
312 */
313 template<class Problem>
314 class TaoGradientCallback {
315 public:
316 static void connect(Tao tao, Problem &p) {
317 PetscErrorCode ierr;
318 ierr = TaoSetGradientRoutine(tao,
319 TaoGradientCallback<Problem>::callback,
320 &p);
321 PISM_CHK(ierr, "TaoSetGradientRoutine");
322 }
323
324 protected:
325 static PetscErrorCode callback(Tao tao, Vec x, Vec gradient, void *ctx) {
326 try {
327 Problem *p = reinterpret_cast<Problem *>(ctx);
328 p->evaluateGradient(tao,x,gradient);
329 } catch (...) {
330 MPI_Comm com = MPI_COMM_SELF;
331 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
332 handle_fatal_errors(com);
333 SETERRQ(com, 1, "A PISM callback failed");
334 }
335 return 0;
336 }
337 };
338
339 //! \brief Adaptor to connect a TAO objective function callback to a C++ object method.
340 /*! The TAO library interfaces with user code via C-style callback functions.
341 This class makes it convenient to associate a TAO convergence monitoring callback
342 with a C++ object method. To assign
343 \code
344 void MyObject::convergenceTest(Tao tao);
345 \endcode
346
347 as the convergence test function to a `Tao` `tao`,
348
349 \code
350 MyObject obj;
351 TaoConvergenceCallback<MyObject>::connect(tao,obj);
352 \endcode
353
354 The method name `convergenceTest` for the callback is hard-coded.
355 See TaoObjGradCallback for a technique to allow
356 the method name to be specified (at the expense of a little more cumbersome code).
357 */
358 template<class Problem>
359 class TaoConvergenceCallback {
360 public:
361 static void connect(Tao tao, Problem &p) {
362 PetscErrorCode ierr;
363 ierr = TaoSetConvergenceTest(tao,
364 TaoConvergenceCallback<Problem>::callback,
365 &p);
366 PISM_CHK(ierr, "TaoSetConvergenceTest");
367 }
368 protected:
369 static PetscErrorCode callback(Tao tao, void *ctx) {
370 try {
371 Problem *p = reinterpret_cast<Problem *>(ctx);
372 p->convergenceTest(tao);
373 } catch (...) {
374 MPI_Comm com = MPI_COMM_SELF;
375 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
376 handle_fatal_errors(com);
377 SETERRQ(com, 1, "A PISM callback failed");
378 }
379 return 0;
380 }
381 };
382
383
384 //! \brief Adaptor to connect a TAO objective and gradient function callback to a C++ object method.
385 /*! The TAO library interfaces with user code via C-style callback functions.
386 This class makes it convenient to associate a TAO combined objective value and gradient
387 callback with a C++ object method. To assign
388 \code
389 void MyObject::someObjectiveFunction(Tao tao,Vec x, double *value, Vec gradient);
390 \endcode
391
392 as the convergence test function to a `Tao` `tao`,
393
394 \code
395 MyObject obj;
396 typedef TaoObjGradCallback<MyObject,&MyObject::someObjectiveFunction> ObjGradCallback;
397 ObjGradCallback::connect(tao,obj);
398 \endcode
399
400 Note that the method name for the callback must be specified explicitly via a template argument.
401 */
402 template<class Problem, void (Problem::*Callback)(Tao,Vec,double*,Vec) >
403 class TaoObjGradCallback {
404 public:
405 static void connect(Tao tao, Problem &p) {
406 PetscErrorCode ierr;
407 ierr = TaoSetObjectiveAndGradientRoutine(tao,
408 TaoObjGradCallback<Problem,Callback>::callback,
409 &p);
410 PISM_CHK(ierr, "TaoSetObjectiveAndGradientRoutine");
411 }
412 protected:
413 static PetscErrorCode callback(Tao tao, Vec x, double *value, Vec gradient, void *ctx) {
414 try {
415 Problem *p = reinterpret_cast<Problem *>(ctx);
416 (p->*Callback)(tao,x,value,gradient);
417 } catch (...) {
418 MPI_Comm com = MPI_COMM_SELF;
419 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
420 handle_fatal_errors(com);
421 SETERRQ(com, 1, "A PISM callback failed");
422 }
423 return 0;
424 }
425 };
426
427 //! \brief Adaptor to connect a TAO objective function callback to a C++ object method.
428 /*! The TAO library interfaces with user code via C-style callback functions.
429 This class makes it convenient to associate a TAO Linearly Constrained Augmented Lagrangian (LCL)
430 callbacks with C++ object methods. To assign
431 \code
432 void MyObject::evaluateConstraints(Tao tao,Vec x,Vec c);
433 void MyObject::evaluateConstraintsJacobianState(Tao tao, Vec x, Mat *J, Mat *Jpc, Mat *Jinv, MatStructure *structure);
434 void MyObject::evaluateConstraintsJacobianDesign(Tao tao, Vec x, Mat *J);
435 \endcode
436 as the LCL callbacks to a `Tao` `tao`,
437
438 \code
439 MyObject obj;
440 TaoLCLCallback<MyObject>::connect(tao,obj);
441 \endcode
442
443 The method names for the callback (`evaluateConstraints`, etc.) are hard-coded.
444 */
445 template<class Problem>
446 class TaoLCLCallbacks {
447 public:
448 static void connect(Tao tao, Problem &p, Vec c, Mat Jc, Mat Jd, Mat Jcpc=NULL, Mat Jcinv=NULL) {
449 PetscErrorCode ierr;
450
451 ierr = TaoSetConstraintsRoutine(tao, c,
452 TaoLCLCallbacks<Problem>::constraints_callback, &p);
453 PISM_CHK(ierr, "TaoSetConstraintsRoutine");
454
455 if (Jcpc == NULL) {
456 Jcpc = Jc;
457 }
458
459 ierr = TaoSetJacobianStateRoutine(tao, Jc, Jcpc, Jcinv,
460 TaoLCLCallbacks<Problem>::jacobian_state_callback, &p);
461 PISM_CHK(ierr, "TaoSetJacobianStateRoutine");
462
463 ierr = TaoSetJacobianDesignRoutine(tao, Jd,
464 TaoLCLCallbacks<Problem>::jacobian_design_callback, &p);
465 PISM_CHK(ierr, "TaoSetJacobianDesignRoutine");
466 }
467 protected:
468 static PetscErrorCode constraints_callback(Tao tao, Vec x, Vec c, void*ctx) {
469 try {
470 Problem *p = reinterpret_cast<Problem *>(ctx);
471 p->evaluateConstraints(tao, x, c);
472 } catch (...) {
473 MPI_Comm com = MPI_COMM_SELF;
474 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
475 handle_fatal_errors(com);
476 SETERRQ(com, 1, "A PISM callback failed");
477 }
478 return 0;
479 }
480
481 static PetscErrorCode jacobian_state_callback(Tao tao, Vec x, Mat J, Mat Jpc,
482 Mat Jinv, void*ctx) {
483 try {
484 Problem *p = reinterpret_cast<Problem *>(ctx);
485 // The MatStructure argument is not used in PETSc 3.5, but I want
486 // to preserve the signature of
487 // evaluateConstraintsJacobianState(...) for now -- (CK)
488 MatStructure structure;
489 p->evaluateConstraintsJacobianState(tao, x, J, Jpc, Jinv, &structure);
490 } catch (...) {
491 MPI_Comm com = MPI_COMM_SELF;
492 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
493 handle_fatal_errors(com);
494 SETERRQ(com, 1, "A PISM callback failed");
495 }
496 return 0;
497 }
498
499 static PetscErrorCode jacobian_design_callback(Tao tao, Vec x, Mat J, void*ctx) {
500 try {
501 Problem *p = reinterpret_cast<Problem *>(ctx);
502 p->evaluateConstraintsJacobianDesign(tao, x, J);
503 } catch (...) {
504 MPI_Comm com = MPI_COMM_SELF;
505 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
506 handle_fatal_errors(com);
507 SETERRQ(com, 1, "A PISM callback failed");
508 }
509 return 0;
510 }
511 };
512
513 } // end of namespace taoutil
514 } // end of namespace pism
515
516 #endif /* end of include guard: TAOUTIL_HH_W42GJNRO */