Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #include "libmesh/libmesh_common.h"
21 :
22 : #if defined(LIBMESH_HAVE_PETSC_TAO) && !defined(LIBMESH_USE_COMPLEX_NUMBERS)
23 :
24 :
25 : // C++ includes
26 :
27 : // Local Includes
28 : #include "libmesh/libmesh_logging.h"
29 : #include "libmesh/petsc_vector.h"
30 : #include "libmesh/petsc_matrix.h"
31 : #include "libmesh/dof_map.h"
32 : #include "libmesh/tao_optimization_solver.h"
33 : #include "libmesh/equation_systems.h"
34 :
35 : namespace libMesh
36 : {
37 :
38 : //--------------------------------------------------------------------
39 : // Functions with C linkage to pass to Tao. Tao will call these
40 : // methods as needed.
41 : //
42 : // Since they must have C linkage they have no knowledge of a namespace.
43 : // Give them an obscure name to avoid namespace pollution.
44 : extern "C"
45 : {
46 :
47 : //---------------------------------------------------------------
48 : // This function is called by Tao to evaluate objective function at x
49 : PetscErrorCode
50 990 : __libmesh_tao_objective (Tao /*tao*/, Vec x, PetscReal * objective, void * ctx)
51 : {
52 : PetscFunctionBegin;
53 :
54 44 : LOG_SCOPE("objective()", "TaoOptimizationSolver");
55 :
56 22 : libmesh_assert(x);
57 22 : libmesh_assert(objective);
58 22 : libmesh_assert(ctx);
59 :
60 : // ctx should be a pointer to the solver (it was passed in as void *)
61 22 : TaoOptimizationSolver<Number> * solver =
62 : static_cast<TaoOptimizationSolver<Number> *> (ctx);
63 :
64 44 : OptimizationSystem & sys = solver->system();
65 :
66 : // We'll use current_local_solution below, so let's ensure that it's consistent
67 : // with the vector x that was passed in.
68 22 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
69 990 : PetscVector<Number> X(x, sys.comm());
70 :
71 : // Perform a swap so that sys.solution points to the input vector
72 : // "x", update sys.current_local_solution based on "x", then swap
73 : // back.
74 990 : X.swap(X_sys);
75 990 : sys.update();
76 990 : X.swap(X_sys);
77 :
78 : // Enforce constraints (if any) exactly on the
79 : // current_local_solution. This is the solution vector that is
80 : // actually used in the computation of the objective function
81 : // below, and is not locked by debug-enabled PETSc the way that
82 : // the solution vector is.
83 990 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
84 :
85 990 : if (solver->objective_object != nullptr)
86 1012 : (*objective) = PS(solver->objective_object->objective(*(sys.current_local_solution), sys));
87 : else
88 0 : libmesh_error_msg("Objective function not defined in __libmesh_tao_objective");
89 :
90 1012 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
91 946 : }
92 :
93 :
94 :
95 : //---------------------------------------------------------------
96 : // This function is called by Tao to evaluate the gradient at x
97 : PetscErrorCode
98 640 : __libmesh_tao_gradient(Tao /*tao*/, Vec x, Vec g, void * ctx)
99 : {
100 : PetscFunctionBegin;
101 :
102 24 : LOG_SCOPE("gradient()", "TaoOptimizationSolver");
103 :
104 12 : libmesh_assert(x);
105 12 : libmesh_assert(g);
106 12 : libmesh_assert(ctx);
107 :
108 : // ctx should be a pointer to the solver (it was passed in as void *)
109 12 : TaoOptimizationSolver<Number> * solver =
110 : static_cast<TaoOptimizationSolver<Number> *> (ctx);
111 :
112 24 : OptimizationSystem & sys = solver->system();
113 :
114 : // We'll use current_local_solution below, so let's ensure that it's consistent
115 : // with the vector x that was passed in.
116 12 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
117 664 : PetscVector<Number> X(x, sys.comm());
118 :
119 : // Perform a swap so that sys.solution points to the input vector
120 : // "x", update sys.current_local_solution based on "x", then swap
121 : // back.
122 640 : X.swap(X_sys);
123 640 : sys.update();
124 640 : X.swap(X_sys);
125 :
126 : // We'll also pass the gradient in to the assembly routine
127 : // so let's make a PETSc vector for that too.
128 652 : PetscVector<Number> gradient(g, sys.comm());
129 :
130 : // Clear the gradient prior to assembly
131 640 : gradient.zero();
132 :
133 : // Enforce constraints exactly on the current_local_solution.
134 640 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
135 :
136 640 : if (solver->gradient_object != nullptr)
137 652 : solver->gradient_object->gradient(*(sys.current_local_solution), gradient, sys);
138 : else
139 0 : libmesh_error_msg("Gradient function not defined in __libmesh_tao_gradient");
140 :
141 640 : gradient.close();
142 :
143 652 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
144 616 : }
145 :
146 : //---------------------------------------------------------------
147 : // This function is called by Tao to evaluate the Hessian at x
148 : PetscErrorCode
149 524 : __libmesh_tao_hessian(Tao /*tao*/, Vec x, Mat h, Mat pc, void * ctx)
150 : {
151 : PetscFunctionBegin;
152 :
153 20 : LOG_SCOPE("hessian()", "TaoOptimizationSolver");
154 :
155 10 : libmesh_assert(x);
156 10 : libmesh_assert(h);
157 10 : libmesh_assert(pc);
158 10 : libmesh_assert(ctx);
159 :
160 : // ctx should be a pointer to the solver (it was passed in as void *)
161 10 : TaoOptimizationSolver<Number> * solver =
162 : static_cast<TaoOptimizationSolver<Number> *> (ctx);
163 :
164 20 : OptimizationSystem & sys = solver->system();
165 :
166 : // We'll use current_local_solution below, so let's ensure that it's consistent
167 : // with the vector x that was passed in.
168 10 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
169 544 : PetscVector<Number> X(x, sys.comm());
170 :
171 : // Perform a swap so that sys.solution points to the input vector
172 : // "x", update sys.current_local_solution based on "x", then swap
173 : // back.
174 524 : X.swap(X_sys);
175 524 : sys.update();
176 524 : X.swap(X_sys);
177 :
178 : // Let's also wrap pc and h in PetscMatrix objects for convenience
179 544 : PetscMatrix<Number> PC(pc, sys.comm());
180 534 : PetscMatrix<Number> hessian(h, sys.comm());
181 524 : PC.attach_dof_map(sys.get_dof_map());
182 524 : hessian.attach_dof_map(sys.get_dof_map());
183 :
184 : // Enforce constraints exactly on the current_local_solution.
185 524 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
186 :
187 524 : if (solver->hessian_object != nullptr)
188 : {
189 : // Following PetscNonlinearSolver by passing in PC. It's not clear
190 : // why we pass in PC and not hessian though?
191 1038 : solver->hessian_object->hessian(*(sys.current_local_solution), PC, sys);
192 : }
193 : else
194 0 : libmesh_error_msg("Hessian function not defined in __libmesh_tao_hessian");
195 :
196 524 : PC.close();
197 524 : hessian.close();
198 :
199 534 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
200 504 : }
201 :
202 :
203 : //---------------------------------------------------------------
204 : // This function is called by Tao to evaluate the equality constraints at x
205 : PetscErrorCode
206 129 : __libmesh_tao_equality_constraints(Tao /*tao*/, Vec x, Vec ce, void * ctx)
207 : {
208 : PetscFunctionBegin;
209 :
210 0 : LOG_SCOPE("equality_constraints()", "TaoOptimizationSolver");
211 :
212 0 : libmesh_assert(x);
213 0 : libmesh_assert(ce);
214 0 : libmesh_assert(ctx);
215 :
216 : // ctx should be a pointer to the solver (it was passed in as void *)
217 0 : TaoOptimizationSolver<Number> * solver =
218 : static_cast<TaoOptimizationSolver<Number> *> (ctx);
219 :
220 0 : OptimizationSystem & sys = solver->system();
221 :
222 : // We'll use current_local_solution below, so let's ensure that it's consistent
223 : // with the vector x that was passed in.
224 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
225 129 : PetscVector<Number> X(x, sys.comm());
226 :
227 : // Perform a swap so that sys.solution points to the input vector
228 : // "x", update sys.current_local_solution based on "x", then swap
229 : // back.
230 129 : X.swap(X_sys);
231 129 : sys.update();
232 129 : X.swap(X_sys);
233 :
234 : // We'll also pass the constraints vector ce into the assembly routine
235 : // so let's make a PETSc vector for that too.
236 129 : PetscVector<Number> eq_constraints(ce, sys.comm());
237 :
238 : // Clear the gradient prior to assembly
239 129 : eq_constraints.zero();
240 :
241 : // Enforce constraints exactly on the current_local_solution.
242 129 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
243 :
244 129 : if (solver->equality_constraints_object != nullptr)
245 129 : solver->equality_constraints_object->equality_constraints(*(sys.current_local_solution), eq_constraints, sys);
246 : else
247 0 : libmesh_error_msg("Constraints function not defined in __libmesh_tao_equality_constraints");
248 :
249 129 : eq_constraints.close();
250 :
251 129 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
252 129 : }
253 :
254 : //---------------------------------------------------------------
255 : // This function is called by Tao to evaluate the Jacobian of the
256 : // equality constraints at x
257 : PetscErrorCode
258 129 : __libmesh_tao_equality_constraints_jacobian(Tao /*tao*/, Vec x, Mat J, Mat Jpre, void * ctx)
259 : {
260 : PetscFunctionBegin;
261 :
262 0 : LOG_SCOPE("equality_constraints_jacobian()", "TaoOptimizationSolver");
263 :
264 0 : libmesh_assert(x);
265 0 : libmesh_assert(J);
266 0 : libmesh_assert(Jpre);
267 :
268 : // ctx should be a pointer to the solver (it was passed in as void *)
269 0 : TaoOptimizationSolver<Number> * solver =
270 : static_cast<TaoOptimizationSolver<Number> *> (ctx);
271 :
272 0 : OptimizationSystem & sys = solver->system();
273 :
274 : // We'll use current_local_solution below, so let's ensure that it's consistent
275 : // with the vector x that was passed in.
276 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
277 129 : PetscVector<Number> X(x, sys.comm());
278 :
279 : // Perform a swap so that sys.solution points to the input vector
280 : // "x", update sys.current_local_solution based on "x", then swap
281 : // back.
282 129 : X.swap(X_sys);
283 129 : sys.update();
284 129 : X.swap(X_sys);
285 :
286 : // Let's also wrap J and Jpre in PetscMatrix objects for convenience
287 129 : PetscMatrix<Number> J_petsc(J, sys.comm());
288 129 : PetscMatrix<Number> Jpre_petsc(Jpre, sys.comm());
289 :
290 : // Enforce constraints exactly on the current_local_solution.
291 129 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
292 :
293 129 : if (solver->equality_constraints_jacobian_object != nullptr)
294 129 : solver->equality_constraints_jacobian_object->equality_constraints_jacobian(*(sys.current_local_solution), J_petsc, sys);
295 : else
296 0 : libmesh_error_msg("Constraints function not defined in __libmesh_tao_equality_constraints_jacobian");
297 :
298 129 : J_petsc.close();
299 129 : Jpre_petsc.close();
300 :
301 129 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
302 129 : }
303 :
304 : //---------------------------------------------------------------
305 : // This function is called by Tao to evaluate the inequality constraints at x
306 : PetscErrorCode
307 129 : __libmesh_tao_inequality_constraints(Tao /*tao*/, Vec x, Vec cineq, void * ctx)
308 : {
309 : PetscFunctionBegin;
310 :
311 0 : LOG_SCOPE("inequality_constraints()", "TaoOptimizationSolver");
312 :
313 0 : libmesh_assert(x);
314 0 : libmesh_assert(cineq);
315 0 : libmesh_assert(ctx);
316 :
317 : // ctx should be a pointer to the solver (it was passed in as void *)
318 0 : TaoOptimizationSolver<Number> * solver =
319 : static_cast<TaoOptimizationSolver<Number> *> (ctx);
320 :
321 0 : OptimizationSystem & sys = solver->system();
322 :
323 : // We'll use current_local_solution below, so let's ensure that it's consistent
324 : // with the vector x that was passed in.
325 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
326 129 : PetscVector<Number> X(x, sys.comm());
327 :
328 : // Perform a swap so that sys.solution points to the input vector
329 : // "x", update sys.current_local_solution based on "x", then swap
330 : // back.
331 129 : X.swap(X_sys);
332 129 : sys.update();
333 129 : X.swap(X_sys);
334 :
335 : // We'll also pass the constraints vector ce into the assembly routine
336 : // so let's make a PETSc vector for that too.
337 129 : PetscVector<Number> ineq_constraints(cineq, sys.comm());
338 :
339 : // Clear the gradient prior to assembly
340 129 : ineq_constraints.zero();
341 :
342 : // Enforce constraints exactly on the current_local_solution.
343 129 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
344 :
345 129 : if (solver->inequality_constraints_object != nullptr)
346 129 : solver->inequality_constraints_object->inequality_constraints(*(sys.current_local_solution), ineq_constraints, sys);
347 : else
348 0 : libmesh_error_msg("Constraints function not defined in __libmesh_tao_inequality_constraints");
349 :
350 129 : ineq_constraints.close();
351 :
352 129 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
353 129 : }
354 :
355 : //---------------------------------------------------------------
356 : // This function is called by Tao to evaluate the Jacobian of the
357 : // equality constraints at x
358 : PetscErrorCode
359 129 : __libmesh_tao_inequality_constraints_jacobian(Tao /*tao*/, Vec x, Mat J, Mat Jpre, void * ctx)
360 : {
361 : PetscFunctionBegin;
362 :
363 0 : LOG_SCOPE("inequality_constraints_jacobian()", "TaoOptimizationSolver");
364 :
365 0 : libmesh_assert(x);
366 0 : libmesh_assert(J);
367 0 : libmesh_assert(Jpre);
368 :
369 : // ctx should be a pointer to the solver (it was passed in as void *)
370 0 : TaoOptimizationSolver<Number> * solver =
371 : static_cast<TaoOptimizationSolver<Number> *> (ctx);
372 :
373 0 : OptimizationSystem & sys = solver->system();
374 :
375 : // We'll use current_local_solution below, so let's ensure that it's consistent
376 : // with the vector x that was passed in.
377 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
378 129 : PetscVector<Number> X(x, sys.comm());
379 :
380 : // Perform a swap so that sys.solution points to the input vector
381 : // "x", update sys.current_local_solution based on "x", then swap
382 : // back.
383 129 : X.swap(X_sys);
384 129 : sys.update();
385 129 : X.swap(X_sys);
386 :
387 : // Let's also wrap J and Jpre in PetscMatrix objects for convenience
388 129 : PetscMatrix<Number> J_petsc(J, sys.comm());
389 129 : PetscMatrix<Number> Jpre_petsc(Jpre, sys.comm());
390 :
391 : // Enforce constraints exactly on the current_local_solution.
392 129 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
393 :
394 129 : if (solver->inequality_constraints_jacobian_object != nullptr)
395 129 : solver->inequality_constraints_jacobian_object->inequality_constraints_jacobian(*(sys.current_local_solution), J_petsc, sys);
396 : else
397 0 : libmesh_error_msg("Constraints function not defined in __libmesh_tao_inequality_constraints_jacobian");
398 :
399 129 : J_petsc.close();
400 129 : Jpre_petsc.close();
401 :
402 129 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
403 129 : }
404 :
405 : } // end extern "C"
406 : //---------------------------------------------------------------------
407 :
408 :
409 :
410 : //---------------------------------------------------------------------
411 : // TaoOptimizationSolver<> methods
412 : template <typename T>
413 142 : TaoOptimizationSolver<T>::TaoOptimizationSolver (OptimizationSystem & system_in) :
414 : OptimizationSolver<T>(system_in),
415 142 : _reason(TAO_CONVERGED_USER) // Arbitrary initial value...
416 : {
417 142 : }
418 :
419 :
420 :
421 : template <typename T>
422 284 : TaoOptimizationSolver<T>::~TaoOptimizationSolver ()
423 : {
424 8 : this->TaoOptimizationSolver::clear ();
425 351 : }
426 :
427 :
428 :
429 : template <typename T>
430 79 : void TaoOptimizationSolver<T>::clear () noexcept
431 : {
432 213 : if (this->initialized())
433 : {
434 71 : this->_is_initialized = false;
435 :
436 71 : PetscErrorCode ierr = TaoDestroy(&_tao);
437 : if (ierr)
438 : libmesh_warning("Warning: TaoDestroy returned a non-zero error code which we ignored.");
439 : }
440 79 : }
441 :
442 :
443 :
444 : template <typename T>
445 142 : void TaoOptimizationSolver<T>::init ()
446 : {
447 : // Initialize the data structures if not done so already.
448 142 : if (!this->initialized())
449 : {
450 71 : this->_is_initialized = true;
451 :
452 71 : LibmeshPetscCall(TaoCreate(this->comm().get(),&_tao));
453 : }
454 142 : }
455 :
456 : template <typename T>
457 71 : void TaoOptimizationSolver<T>::solve ()
458 : {
459 4 : LOG_SCOPE("solve()", "TaoOptimizationSolver");
460 :
461 71 : this->init ();
462 :
463 73 : this->system().solution->zero();
464 :
465 71 : PetscMatrixBase<T> * hessian = cast_ptr<PetscMatrixBase<T> *>(this->system().matrix);
466 : // PetscVector<T> * gradient = cast_ptr<PetscVector<T> *>(this->system().rhs);
467 2 : PetscVector<T> * x = cast_ptr<PetscVector<T> *>(this->system().solution.get());
468 2 : PetscVector<T> * ceq = cast_ptr<PetscVector<T> *>(this->system().C_eq.get());
469 2 : PetscMatrixBase<T> * ceq_jac = cast_ptr<PetscMatrixBase<T> *>(this->system().C_eq_jac.get());
470 2 : PetscVector<T> * cineq = cast_ptr<PetscVector<T> *>(this->system().C_ineq.get());
471 2 : PetscMatrixBase<T> * cineq_jac = cast_ptr<PetscMatrixBase<T> *>(this->system().C_ineq_jac.get());
472 71 : PetscVector<T> * lb = cast_ptr<PetscVector<T> *>(&this->system().get_vector("lower_bounds"));
473 71 : PetscVector<T> * ub = cast_ptr<PetscVector<T> *>(&this->system().get_vector("upper_bounds"));
474 :
475 : // Set the starting guess to zero.
476 71 : x->zero();
477 :
478 : // Workaround for bug where TaoSetFromOptions *reset*
479 : // programmatically set tolerance and max. function evaluation
480 : // values when "-tao_type ipm" was specified on the command line: we
481 : // call TaoSetFromOptions twice (both before and after setting
482 : // custom options programmatically)
483 71 : LibmeshPetscCall(TaoSetFromOptions(_tao));
484 :
485 : // Set convergence tolerances
486 : // f(X) - f(X*) (estimated) <= fatol
487 : // |f(X) - f(X*)| (estimated) / |f(X)| <= frtol
488 : // ||g(X)|| <= gatol
489 : // ||g(X)|| / |f(X)| <= grtol
490 : // ||g(X)|| / ||g(X0)|| <= gttol
491 : // Command line equivalents: -tao_fatol, -tao_frtol, -tao_gatol, -tao_grtol, -tao_gttol
492 : // Releases up to 3.7.0 had fatol and frtol, after that they were removed.
493 : #if PETSC_VERSION_LESS_THAN(3,7,0)
494 : LibmeshPetscCall(TaoSetTolerances(_tao,
495 : /*fatol=*/PETSC_DEFAULT,
496 : /*frtol=*/PETSC_DEFAULT,
497 : /*gatol=*/PETSC_DEFAULT,
498 : /*grtol=*/this->objective_function_relative_tolerance,
499 : /*gttol=*/PETSC_DEFAULT));
500 : #else
501 71 : LibmeshPetscCall(TaoSetTolerances(_tao,
502 : /*gatol=*/PETSC_DEFAULT,
503 : /*grtol=*/this->objective_function_relative_tolerance,
504 : /*gttol=*/PETSC_DEFAULT));
505 : #endif
506 :
507 : // Set the max-allowed number of objective function evaluations
508 : // Command line equivalent: -tao_max_funcs
509 71 : LibmeshPetscCall(TaoSetMaximumFunctionEvaluations(_tao, this->max_objective_function_evaluations));
510 :
511 : // Set the max-allowed number of optimization iterations.
512 : // Command line equivalent: -tao_max_it
513 : // Not implemented for now as it seems fairly similar to
514 : // LibmeshPetscCall(TaoSetMaximumIterations(_tao, 4));
515 :
516 : // Set solution vec and an initial guess
517 : #if PETSC_VERSION_LESS_THAN(3,17,0)
518 4 : LibmeshPetscCall(TaoSetInitialVector(_tao, x->vec()));
519 : #else
520 67 : LibmeshPetscCall(TaoSetSolution(_tao, x->vec()));
521 : #endif
522 :
523 : // We have to have an objective function
524 2 : libmesh_assert( this->objective_object );
525 :
526 : // Set routines for objective, gradient, hessian evaluation
527 : #if PETSC_VERSION_LESS_THAN(3,17,0)
528 4 : LibmeshPetscCall(TaoSetObjectiveRoutine(_tao, __libmesh_tao_objective, this));
529 : #else
530 67 : LibmeshPetscCall(TaoSetObjective(_tao, __libmesh_tao_objective, this));
531 : #endif
532 :
533 71 : if (this->gradient_object)
534 : #if PETSC_VERSION_LESS_THAN(3,17,0)
535 4 : LibmeshPetscCall(TaoSetGradientRoutine(_tao, __libmesh_tao_gradient, this));
536 : #else
537 67 : LibmeshPetscCall(TaoSetGradient(_tao, NULL, __libmesh_tao_gradient, this));
538 : #endif
539 :
540 71 : if (this->hessian_object)
541 : #if PETSC_VERSION_LESS_THAN(3,17,0)
542 4 : LibmeshPetscCall(TaoSetHessianRoutine(_tao, hessian->mat(), hessian->mat(), __libmesh_tao_hessian, this));
543 : #else
544 67 : LibmeshPetscCall(TaoSetHessian(_tao, hessian->mat(), hessian->mat(), __libmesh_tao_hessian, this));
545 : #endif
546 :
547 71 : if (this->lower_and_upper_bounds_object)
548 : {
549 : // Need to actually compute the bounds vectors first
550 1 : this->lower_and_upper_bounds_object->lower_and_upper_bounds(this->system());
551 :
552 1 : LibmeshPetscCall(TaoSetVariableBounds(_tao,
553 : lb->vec(),
554 : ub->vec()));
555 : }
556 :
557 71 : if (this->equality_constraints_object)
558 1 : LibmeshPetscCall(TaoSetEqualityConstraintsRoutine(_tao, ceq->vec(), __libmesh_tao_equality_constraints, this));
559 :
560 71 : if (this->equality_constraints_jacobian_object)
561 1 : LibmeshPetscCall(TaoSetJacobianEqualityRoutine(_tao,
562 : ceq_jac->mat(),
563 : ceq_jac->mat(),
564 : __libmesh_tao_equality_constraints_jacobian,
565 : this));
566 :
567 : // Optionally set inequality constraints
568 71 : if (this->inequality_constraints_object)
569 1 : LibmeshPetscCall(TaoSetInequalityConstraintsRoutine(_tao, cineq->vec(), __libmesh_tao_inequality_constraints, this));
570 :
571 : // Optionally set inequality constraints Jacobian
572 71 : if (this->inequality_constraints_jacobian_object)
573 1 : LibmeshPetscCall(TaoSetJacobianInequalityRoutine(_tao,
574 : cineq_jac->mat(),
575 : cineq_jac->mat(),
576 : __libmesh_tao_inequality_constraints_jacobian,
577 : this));
578 :
579 : // Check for Tao command line options
580 71 : LibmeshPetscCall(TaoSetFromOptions(_tao));
581 :
582 : // Perform the optimization
583 71 : LibmeshPetscCall(TaoSolve(_tao));
584 :
585 : // Enforce constraints exactly now that the solve is done. We have
586 : // been enforcing them on the current_local_solution during the
587 : // solve, but now need to be sure they are enforced on the parallel
588 : // solution vector as well.
589 73 : this->system().get_dof_map().enforce_constraints_exactly(this->system());
590 :
591 : // Store the convergence/divergence reason
592 71 : LibmeshPetscCall(TaoGetConvergedReason(_tao, &_reason));
593 71 : }
594 :
595 :
596 : template <typename T>
597 129 : void TaoOptimizationSolver<T>::get_dual_variables()
598 : {
599 0 : LOG_SCOPE("get_dual_variables()", "TaoOptimizationSolver");
600 :
601 : PetscVector<T> * lambda_eq_petsc =
602 0 : cast_ptr<PetscVector<T> *>(this->system().lambda_eq.get());
603 : PetscVector<T> * lambda_ineq_petsc =
604 0 : cast_ptr<PetscVector<T> *>(this->system().lambda_ineq.get());
605 :
606 129 : Vec lambda_eq_petsc_vec = lambda_eq_petsc->vec();
607 129 : Vec lambda_ineq_petsc_vec = lambda_ineq_petsc->vec();
608 :
609 129 : LibmeshPetscCall(TaoGetDualVariables(_tao,
610 : &lambda_eq_petsc_vec,
611 : &lambda_ineq_petsc_vec));
612 129 : }
613 :
614 :
615 : template <typename T>
616 71 : void TaoOptimizationSolver<T>::print_converged_reason()
617 : {
618 2 : libMesh::out << "Tao optimization solver convergence/divergence reason: "
619 71 : << TaoConvergedReasons[this->get_converged_reason()] << std::endl;
620 71 : }
621 :
622 :
623 :
624 : template <typename T>
625 72 : int TaoOptimizationSolver<T>::get_converged_reason()
626 : {
627 72 : if (this->initialized())
628 72 : LibmeshPetscCall(TaoGetConvergedReason(_tao, &_reason));
629 :
630 72 : return static_cast<int>(_reason);
631 : }
632 :
633 :
634 : //------------------------------------------------------------------
635 : // Explicit instantiations
636 : template class LIBMESH_EXPORT TaoOptimizationSolver<Number>;
637 :
638 : } // namespace libMesh
639 :
640 :
641 :
642 : #endif // #if defined(LIBMESH_HAVE_PETSC_TAO) && !defined(LIBMESH_USE_COMPLEX_NUMBERS)
|
