PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
TaoUtil.hh
Go to the documentation of this file.
1 // Copyright (C) 2012, 2013, 2014, 2015, 2017, 2022, 2023 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.
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>
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 std::shared_ptr<TerminationReason> solve() {
118  PetscErrorCode ierr;
119 
120  /* Solve the application */
121  Vec x0;
122  auto reason = m_problem.formInitialGuess(&x0);
123  if (reason->failed()) {
124  auto 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 #if PETSC_VERSION_LT(3,17,0)
130  ierr = TaoSetInitialVector(m_tao, x0);
131 #else
132  ierr = TaoSetSolution(m_tao, x0);
133 #endif
134  PISM_CHK(ierr, "TaoSetInitialVector");
135 
136  ierr = TaoSolve(m_tao);
137  PISM_CHK(ierr, "TaoSolve");
138 
139  TaoConvergedReason tao_reason;
140  ierr = TaoGetConvergedReason(m_tao, &tao_reason);
141  PISM_CHK(ierr, "TaoGetConvergedReason");
142 
143  reason.reset(new TAOTerminationReason(tao_reason));
144  return reason;
145  }
146 
147  virtual void setMaximumIterations(int max_it) {
148  PetscErrorCode ierr = TaoSetMaximumIterations(m_tao, max_it);
149  PISM_CHK(ierr, "TaoSetMaximumIterations");
150  }
151 
152  virtual Problem &problem() {
153  return m_problem;
154  }
155 
156 protected:
157  MPI_Comm m_comm;
159  Problem &m_problem;
160 };
161 
162 
163 //! \brief Adaptor to connect a TAO objective function callback to a C++ object method.
164 /*! The TAO library interfaces with user code via C-style callback functions.
165  This class makes it convenient to associate a TAO Objective callback
166  with a C++ object method. To assign
167  \code
168  void MyObject::evaluateObjective(TaoSolver tao,Vec x, double *value);
169  \endcode
170 
171  as the objective function to a `TaoSolver` `tao`,
172 
173  \code
174  MyObject obj;
175  TaoObjectiveCallback<MyObject>::connect(tao,obj);
176  \endcode
177 
178  The method name `evaluateObjective` for the callback is hard-coded.
179  See TaoObjGradCallback for a technique to allow
180  the method name to be specified (at the expense of a little more cumbersome code).
181 */
182 template<class Problem>
184 public:
185 
186  static void connect(Tao tao, Problem &p) {
187  PetscErrorCode ierr;
188  ierr = TaoSetObjectiveRoutine(tao,
190  &p);
191  PISM_CHK(ierr, "TaoSetObjectiveRoutine");
192  }
193 
194 protected:
195  static PetscErrorCode callback(Tao tao, Vec x, double *value, void *ctx) {
196  try {
197  Problem *p = reinterpret_cast<Problem *>(ctx);
198  PetscErrorCode ierr = p->evaluateObjective(tao,x,value); CHKERRQ(ierr);
199  } catch (...) {
200  MPI_Comm com = MPI_COMM_SELF;
201  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
202  handle_fatal_errors(com);
203  SETERRQ(com, 1, "A PISM callback failed");
204  }
205  return 0;
206  }
207 };
208 
209 
210 //! \brief Adaptor to connect a TAO monitoring callback to a C++ object method.
211 /*! The TAO library interfaces with user code via C-style callback functions.
212  This class makes it convenient to associate a TAO Monitor callback
213  with a C++ object method. To assign
214  \code
215  void MyObject::monitorTao(Tao tao)
216  \endcode
217 
218  as the objective function to a `Tao` `tao`,
219 
220  \code
221  MyObject obj;
222  TaoMonitorCallback<MyObject>::connect(tao,obj);
223  \endcode
224 
225  The method name `monitorTao` for the callback is hard-coded.
226  See TaoObjGradCallback for a technique to allow
227  the method name to be specified (at the expense of a little more cumbersome code).
228 */
229 template<class Problem>
231 public:
232  static void connect(Tao tao, Problem &p) {
233  PetscErrorCode ierr = TaoSetMonitor(tao,
235  &p, NULL);
236  PISM_CHK(ierr, "TaoSetMonitor");
237  }
238 protected:
239  static PetscErrorCode callback(Tao tao, void *ctx) {
240  try {
241  Problem *p = reinterpret_cast<Problem *>(ctx);
242  p->monitorTao(tao);
243  } catch (...) {
244  MPI_Comm com = MPI_COMM_SELF;
245  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
246  handle_fatal_errors(com);
247  SETERRQ(com, 1, "A PISM callback failed");
248  }
249  return 0;
250  }
251 };
252 
253 //! \brief Adaptor to connect a TAO objective function callback to a C++ object method.
254 /*! The TAO library interfaces with user code via C-style callback functions.
255  This class makes it convenient to associate a TAO VariableBounds callback
256  with a C++ object method. To assign
257  \code
258  void MyObject::getVariableBounds(Tao tao,Vec lo, Vec hi);
259  \endcode
260 
261  as the objective function to a `Tao` `tao`,
262 
263  \code
264  MyObject obj;
265  TaoGetVariableBoundsCallback<MyObject>::connect(tao,obj);
266  \endcode
267 
268  The method name `getVariableBounds` for the callback is hard-coded.
269  See TaoObjGradCallback for a technique to allow
270  the method name to be specified (at the expense of a little more cumbersome code).
271 */
272 template<class Problem>
274 public:
275  static void connect(Tao tao, Problem &p) {
276  PetscErrorCode ierr;
277  ierr = TaoSetVariableBoundsRoutine(tao,
279  &p);
280  PISM_CHK(ierr, "TaoSetVariableBoundsRoutine");
281  }
282 
283 protected:
284  static PetscErrorCode callback(Tao tao, Vec lo, Vec hi, void *ctx) {
285  try {
286  Problem *p = reinterpret_cast<Problem *>(ctx);
287  p->getVariableBounds(tao,lo,hi);
288  } catch (...) {
289  MPI_Comm com = MPI_COMM_SELF;
290  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
291  handle_fatal_errors(com);
292  SETERRQ(com, 1, "A PISM callback failed");
293  }
294  return 0;
295  }
296 };
297 
298 //! \brief Adaptor to connect a TAO objective gradient callback to a C++ object method.
299 /*! The TAO library interfaces with user code via C-style callback functions.
300  This class makes it convenient to associate a TAO Objective Gradient callback
301  with a C++ object method. To assign
302  \code
303  void MyObject::evaluateGradient(Tao tao,Vec x, Vec gradient);
304  \endcode
305 
306  as the objective function to a `Tao` `tao`,
307 
308  \code
309  MyObject obj;
310  TaoGradientCallback<MyObject>::connect(tao,obj);
311  \endcode
312 
313  The method name `evaluateGradient` for the callback is hard-coded.
314  See TaoObjGradCallback for a technique to allow
315  the method name to be specified (at the expense of a little more cumbersome code).
316 */
317 template<class Problem>
319 public:
320  static void connect(Tao tao, Problem &p) {
321  PetscErrorCode ierr;
322  ierr = TaoSetGradientRoutine(tao,
324  &p);
325  PISM_CHK(ierr, "TaoSetGradientRoutine");
326  }
327 
328 protected:
329  static PetscErrorCode callback(Tao tao, Vec x, Vec gradient, void *ctx) {
330  try {
331  Problem *p = reinterpret_cast<Problem *>(ctx);
332  p->evaluateGradient(tao,x,gradient);
333  } catch (...) {
334  MPI_Comm com = MPI_COMM_SELF;
335  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
336  handle_fatal_errors(com);
337  SETERRQ(com, 1, "A PISM callback failed");
338  }
339  return 0;
340  }
341 };
342 
343 //! \brief Adaptor to connect a TAO objective function callback to a C++ object method.
344 /*! The TAO library interfaces with user code via C-style callback functions.
345  This class makes it convenient to associate a TAO convergence monitoring callback
346  with a C++ object method. To assign
347  \code
348  void MyObject::convergenceTest(Tao tao);
349  \endcode
350 
351  as the convergence test function to a `Tao` `tao`,
352 
353  \code
354  MyObject obj;
355  TaoConvergenceCallback<MyObject>::connect(tao,obj);
356  \endcode
357 
358  The method name `convergenceTest` for the callback is hard-coded.
359  See TaoObjGradCallback for a technique to allow
360  the method name to be specified (at the expense of a little more cumbersome code).
361 */
362 template<class Problem>
364 public:
365  static void connect(Tao tao, Problem &p) {
366  PetscErrorCode ierr;
367  ierr = TaoSetConvergenceTest(tao,
369  &p);
370  PISM_CHK(ierr, "TaoSetConvergenceTest");
371  }
372 protected:
373  static PetscErrorCode callback(Tao tao, void *ctx) {
374  try {
375  Problem *p = reinterpret_cast<Problem *>(ctx);
376  p->convergenceTest(tao);
377  } catch (...) {
378  MPI_Comm com = MPI_COMM_SELF;
379  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
380  handle_fatal_errors(com);
381  SETERRQ(com, 1, "A PISM callback failed");
382  }
383  return 0;
384  }
385 };
386 
387 
388 //! \brief Adaptor to connect a TAO objective and gradient function callback to a C++ object method.
389 /*! The TAO library interfaces with user code via C-style callback functions.
390  This class makes it convenient to associate a TAO combined objective value and gradient
391  callback with a C++ object method. To assign
392  \code
393  void MyObject::someObjectiveFunction(Tao tao,Vec x, double *value, Vec gradient);
394  \endcode
395 
396  as the convergence test function to a `Tao` `tao`,
397 
398  \code
399  MyObject obj;
400  typedef TaoObjGradCallback<MyObject,&MyObject::someObjectiveFunction> ObjGradCallback;
401  ObjGradCallback::connect(tao,obj);
402  \endcode
403 
404  Note that the method name for the callback must be specified explicitly via a template argument.
405 */
406 template<class Problem, void (Problem::*Callback)(Tao,Vec,double*,Vec) >
408 public:
409  static void connect(Tao tao, Problem &p) {
410  PetscErrorCode ierr;
411 #if PETSC_VERSION_LT(3,17,0)
412  ierr = TaoSetObjectiveAndGradientRoutine(tao,
414  &p);
415 #else
416  ierr = TaoSetObjectiveAndGradient(tao,
417  NULL,
419  &p);
420 #endif
421  PISM_CHK(ierr, "TaoSetObjectiveAndGradientRoutine");
422  }
423 protected:
424  static PetscErrorCode callback(Tao tao, Vec x, double *value, Vec gradient, void *ctx) {
425  try {
426  Problem *p = reinterpret_cast<Problem *>(ctx);
427  (p->*Callback)(tao,x,value,gradient);
428  } catch (...) {
429  MPI_Comm com = MPI_COMM_SELF;
430  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
431  handle_fatal_errors(com);
432  SETERRQ(com, 1, "A PISM callback failed");
433  }
434  return 0;
435  }
436 };
437 
438 //! \brief Adaptor to connect a TAO objective function callback to a C++ object method.
439 /*! The TAO library interfaces with user code via C-style callback functions.
440  This class makes it convenient to associate a TAO Linearly Constrained Augmented Lagrangian (LCL)
441  callbacks with C++ object methods. To assign
442  \code
443  void MyObject::evaluateConstraints(Tao tao,Vec x,Vec c);
444  void MyObject::evaluateConstraintsJacobianState(Tao tao, Vec x, Mat *J, Mat *Jpc, Mat *Jinv, MatStructure *structure);
445  void MyObject::evaluateConstraintsJacobianDesign(Tao tao, Vec x, Mat *J);
446  \endcode
447  as the LCL callbacks to a `Tao` `tao`,
448 
449  \code
450  MyObject obj;
451  TaoLCLCallback<MyObject>::connect(tao,obj);
452  \endcode
453 
454  The method names for the callback (`evaluateConstraints`, etc.) are hard-coded.
455 */
456 template<class Problem>
458 public:
459  static void connect(Tao tao, Problem &p, Vec c, Mat Jc, Mat Jd, Mat Jcpc=NULL, Mat Jcinv=NULL) {
460  PetscErrorCode ierr;
461 
462  ierr = TaoSetConstraintsRoutine(tao, c,
464  PISM_CHK(ierr, "TaoSetConstraintsRoutine");
465 
466  if (Jcpc == NULL) {
467  Jcpc = Jc;
468  }
469 
470  ierr = TaoSetJacobianStateRoutine(tao, Jc, Jcpc, Jcinv,
472  PISM_CHK(ierr, "TaoSetJacobianStateRoutine");
473 
474  ierr = TaoSetJacobianDesignRoutine(tao, Jd,
476  PISM_CHK(ierr, "TaoSetJacobianDesignRoutine");
477  }
478 protected:
479  static PetscErrorCode constraints_callback(Tao tao, Vec x, Vec c, void*ctx) {
480  try {
481  Problem *p = reinterpret_cast<Problem *>(ctx);
482  p->evaluateConstraints(tao, x, c);
483  } catch (...) {
484  MPI_Comm com = MPI_COMM_SELF;
485  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
486  handle_fatal_errors(com);
487  SETERRQ(com, 1, "A PISM callback failed");
488  }
489  return 0;
490  }
491 
492  static PetscErrorCode jacobian_state_callback(Tao tao, Vec x, Mat J, Mat Jpc,
493  Mat Jinv, void*ctx) {
494  try {
495  Problem *p = reinterpret_cast<Problem *>(ctx);
496  // The MatStructure argument is not used in PETSc 3.5, but I want
497  // to preserve the signature of
498  // evaluateConstraintsJacobianState(...) for now -- (CK)
499  MatStructure structure;
500  p->evaluateConstraintsJacobianState(tao, x, J, Jpc, Jinv, &structure);
501  } catch (...) {
502  MPI_Comm com = MPI_COMM_SELF;
503  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
504  handle_fatal_errors(com);
505  SETERRQ(com, 1, "A PISM callback failed");
506  }
507  return 0;
508  }
509 
510  static PetscErrorCode jacobian_design_callback(Tao tao, Vec x, Mat J, void*ctx) {
511  try {
512  Problem *p = reinterpret_cast<Problem *>(ctx);
513  p->evaluateConstraintsJacobianDesign(tao, x, J);
514  } catch (...) {
515  MPI_Comm com = MPI_COMM_SELF;
516  PetscErrorCode ierr = PetscObjectGetComm((PetscObject)tao, &com); CHKERRQ(ierr);
517  handle_fatal_errors(com);
518  SETERRQ(com, 1, "A PISM callback failed");
519  }
520  return 0;
521  }
522 };
523 
524 } // end of namespace taoutil
525 } // end of namespace pism
526 
527 #endif /* end of include guard: TAOUTIL_HH_W42GJNRO */
T * rawptr()
Definition: Wrapper.hh:39
TAOTerminationReason(TaoConvergedReason r)
Definition: TaoUtil.cc:24
virtual void get_description(std::ostream &desc, int indent_level=0)
Definition: TaoUtil.cc:28
Encapsulate TAO's TaoSolverTerminationReason codes as a PISM TerminationReason.
Definition: TaoUtil.hh:37
virtual std::shared_ptr< TerminationReason > solve()
Solve the minimization problem.
Definition: TaoUtil.hh:117
virtual Problem & problem()
Definition: TaoUtil.hh:152
TaoBasicSolver(MPI_Comm comm, const std::string &tao_type, Problem &prob)
Construct a solver to solve prob using TAO algorithm tao_type.
Definition: TaoUtil.hh:96
virtual void setMaximumIterations(int max_it)
Definition: TaoUtil.hh:147
An interface for solving an optimization problem with TAO where the problem itself is defined by a se...
Definition: TaoUtil.hh:92
static void connect(Tao tao, Problem &p)
Definition: TaoUtil.hh:365
static PetscErrorCode callback(Tao tao, void *ctx)
Definition: TaoUtil.hh:373
Adaptor to connect a TAO objective function callback to a C++ object method.
Definition: TaoUtil.hh:363
static PetscErrorCode callback(Tao tao, Vec lo, Vec hi, void *ctx)
Definition: TaoUtil.hh:284
static void connect(Tao tao, Problem &p)
Definition: TaoUtil.hh:275
Adaptor to connect a TAO objective function callback to a C++ object method.
Definition: TaoUtil.hh:273
static PetscErrorCode callback(Tao tao, Vec x, Vec gradient, void *ctx)
Definition: TaoUtil.hh:329
static void connect(Tao tao, Problem &p)
Definition: TaoUtil.hh:320
Adaptor to connect a TAO objective gradient callback to a C++ object method.
Definition: TaoUtil.hh:318
static PetscErrorCode constraints_callback(Tao tao, Vec x, Vec c, void *ctx)
Definition: TaoUtil.hh:479
static void connect(Tao tao, Problem &p, Vec c, Mat Jc, Mat Jd, Mat Jcpc=NULL, Mat Jcinv=NULL)
Definition: TaoUtil.hh:459
static PetscErrorCode jacobian_design_callback(Tao tao, Vec x, Mat J, void *ctx)
Definition: TaoUtil.hh:510
static PetscErrorCode jacobian_state_callback(Tao tao, Vec x, Mat J, Mat Jpc, Mat Jinv, void *ctx)
Definition: TaoUtil.hh:492
Adaptor to connect a TAO objective function callback to a C++ object method.
Definition: TaoUtil.hh:457
static PetscErrorCode callback(Tao tao, void *ctx)
Definition: TaoUtil.hh:239
static void connect(Tao tao, Problem &p)
Definition: TaoUtil.hh:232
Adaptor to connect a TAO monitoring callback to a C++ object method.
Definition: TaoUtil.hh:230
static PetscErrorCode callback(Tao tao, Vec x, double *value, Vec gradient, void *ctx)
Definition: TaoUtil.hh:424
static void connect(Tao tao, Problem &p)
Definition: TaoUtil.hh:409
Adaptor to connect a TAO objective and gradient function callback to a C++ object method.
Definition: TaoUtil.hh:407
static void connect(Tao tao, Problem &p)
Definition: TaoUtil.hh:186
static PetscErrorCode callback(Tao tao, Vec x, double *value, void *ctx)
Definition: TaoUtil.hh:195
Adaptor to connect a TAO objective function callback to a C++ object method.
Definition: TaoUtil.hh:183
#define PISM_CHK(errcode, name)
void handle_fatal_errors(MPI_Comm com)
const int J[]
Definition: ssafd_code.cc:34