v0.16.3
Loading...
Searching...
No Matches
TaoNullSpace.cpp
Go to the documentation of this file.
1/** \file TaoNullSpace.cpp
2 * \brief Registration and lifecycle of the null-space TAO shell.
3 */
4
5// TAOSHELL exposes no setup/options/view or context-destruction setters.
6// Private access also safeguards rejected More-Thuente trials: PETSc exposes
7// no public callback-domain rejection or line-search apply-hook setter.
8#include <petsc/private/taoimpl.h>
9#include <petsc/private/taolinesearchimpl.h>
10#include <TaoNullSpace.hpp>
11
12#include <memory>
13
14namespace MoFEM {
15namespace {
16
17constexpr char nullSpaceContextKey[] = "_MoFEM_TaoNullSpace";
18constexpr char nullSpaceLineSearch[] = "nullspace";
19
20/** Algorithm 1: halve a fixed flow step, accepting the final finite trial. */
21MoFEMErrorCode createNullSpaceLineSearch(TaoLineSearch ls) {
22 MoFEMFunctionBegin;
23 ls->ops->setup = nullptr;
24 ls->ops->reset = nullptr;
25 ls->ops->monitor = nullptr;
26 ls->ops->view = nullptr;
27 ls->ops->destroy = nullptr;
28 ls->ops->setfromoptions = nullptr;
29 ls->ops->apply = [](TaoLineSearch ls, Vec x, PetscReal *f, Vec g,
30 Vec direction) {
31 MoFEMFunctionBegin;
32 if (!std::isfinite(ls->initstep) || !(ls->initstep > 0) ||
33 ls->max_funcs < 1)
34 SETERRQ(PetscObjectComm(getPetscObject(ls)), PETSC_ERR_ARG_OUTOFRANGE,
35 "Null-space backtracking requires positive stepinit and "
36 "max_funcs >= 1");
37 auto start = vectorDuplicate(x);
38 auto start_gradient = vectorDuplicate(g);
39 CHKERR VecCopy(x, start);
40 CHKERR VecCopy(g, start_gradient);
41 const auto initial_value = *f;
42 CHKERR TaoLineSearchMonitor(ls, 0, *f, 0);
43 ls->step = std::min(ls->initstep, ls->stepmax);
44 ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
45 for (PetscInt trial = 1;
46 trial <= ls->max_funcs && ls->step >= ls->stepmin && ls->step > 0;
47 ++trial, ls->step *= 0.5) {
48 CHKERR VecWAXPY(x, ls->step, direction, start);
49 CHKERR TaoLineSearchComputeObjective(ls, x, f);
50 if (trial == 1)
51 ls->f_fullstep = *f;
52 CHKERR TaoLineSearchMonitor(ls, trial, *f, ls->step);
53 if (std::isfinite(*f) && std::abs(*f) < PETSC_INFINITY &&
54 (*f < initial_value || trial == ls->max_funcs)) {
55 // Only a merit-acceptable trial needs a derivative. The combined
56 // callback also preserves the rejection sentinel if that fails.
57 CHKERR TaoLineSearchComputeObjectiveAndGradient(ls, x, f, g);
58 PetscReal norm;
59 CHKERR VecNorm(g, NORM_2, &norm);
60 if (std::isfinite(*f) && std::abs(*f) < PETSC_INFINITY &&
61 std::isfinite(norm)) {
62 ls->reason = *f < initial_value ? TAOLINESEARCH_SUCCESS
63 : TAOLINESEARCH_SUCCESS_USER;
64 if (ls->reason == TAOLINESEARCH_SUCCESS_USER)
65 CHKERR PetscInfo(ls, "Null-space backtracking accepted the final "
66 "finite trial without merit decrease\n");
67 MoFEMFunctionReturnHot(0);
68 }
69 }
70 if (trial == ls->max_funcs) {
71 ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
72 break;
73 }
74 }
75 // A failed physical evaluation is never accepted by the budget fallback.
76 CHKERR VecCopy(start, x);
77 CHKERR VecCopy(start_gradient, g);
78 *f = initial_value;
79 ls->step = 0;
80 MoFEMFunctionReturn(0);
81 };
82 MoFEMFunctionReturn(0);
83}
84
85/** One Euler step; returning to the parent refreshes all active-set data. */
86MoFEMErrorCode solveNullSpaceFlowStep(Tao tao) {
87 MoFEMFunctionBegin;
88 Vec x, g;
89 CHKERR TaoGetSolution(tao, &x);
90 CHKERR TaoGetGradient(tao, &g, nullptr, nullptr);
91 PetscReal objective, norm, step;
92 CHKERR TaoComputeObjectiveAndGradient(tao, x, &objective, g);
93 CHKERR VecNorm(g, NORM_2, &norm);
94 CHKERR TaoSetConvergedReason(tao, TAO_CONTINUE_ITERATING);
95 CHKERR TaoSetIterationNumber(tao, 0);
96 CHKERR TaoLogConvergenceHistory(tao, objective, norm, 0, 0);
97 CHKERR TaoMonitor(tao, 0, objective, norm, 0, 0);
98 auto direction = vectorDuplicate(x);
99 CHKERR VecCopy(g, direction);
100 CHKERR VecScale(direction, -1);
101 TaoLineSearch ls;
102 CHKERR TaoGetLineSearch(tao, &ls);
103 CHKERR TaoLineSearchUseTaoRoutines(ls, tao);
104 TaoLineSearchConvergedReason ls_reason;
105 CHKERR TaoLineSearchApply(ls, x, &objective, g, direction, &step, &ls_reason);
106 tao->lsflag = ls_reason > 0 ? PETSC_FALSE : PETSC_TRUE;
107 CHKERR TaoAddLineSearchCounts(tao);
108 CHKERR VecNorm(g, NORM_2, &norm);
109 CHKERR TaoSetIterationNumber(tao, 1);
110 CHKERR TaoSetConvergedReason(tao, ls_reason > 0 ? TAO_CONVERGED_USER
111 : TAO_DIVERGED_LS_FAILURE);
112 CHKERR TaoLogConvergenceHistory(tao, objective, norm, 0, 0);
113 CHKERR TaoMonitor(tao, 1, objective, norm, 0, step);
114 MoFEMFunctionReturn(0);
115}
116
117/** Ordered constraint rows and an independent accepted-iterate Jacobian. */
118struct NullSpaceRows {
119 SmartPetscObj<IS> rows;
120 SmartPetscObj<Mat> jacobian;
121 SmartPetscObj<Mat> inverse;
122 SmartPetscObj<Vec> values;
123 SmartPetscObj<Vec> work;
124 boost::shared_ptr<ConstrainMatrixCtx> projection;
125 PetscInt size = 0;
126
127 MoFEMErrorCode select(Mat full_jacobian, IS columns,
128 const std::vector<PetscInt> &indices) {
129 MoFEMFunctionBegin;
130 *this = NullSpaceRows();
131 const auto comm = PetscObjectComm(getPetscObject(full_jacobian));
132 IS is;
133 CHKERR ISCreateGeneral(comm, indices.size(), indices.data(),
134 PETSC_COPY_VALUES, &is);
135 rows = SmartPetscObj<IS>(is);
136 CHKERR ISGetSize(rows, &size);
137 if (size) {
138 Mat matrix;
139 CHKERR MatCreateSubMatrix(full_jacobian, rows, columns,
140 MAT_INITIAL_MATRIX, &matrix);
141 jacobian = SmartPetscObj<Mat>(matrix);
142 Vec vector;
143 CHKERR MatCreateVecs(jacobian, nullptr, &vector);
144 values = SmartPetscObj<Vec>(vector);
145 work = vectorDuplicate(values);
146 }
147 MoFEMFunctionReturn(0);
148 }
149
150 /** Copy a selected subvector before applying potentially failing solvers. */
151 MoFEMErrorCode extract(Vec full, Vec selected) {
152 MoFEMFunctionBegin;
153 Vec sub;
154 CHKERR VecGetSubVector(full, rows, &sub);
155 CHKERR VecCopy(sub, selected);
156 CHKERR VecRestoreSubVector(full, rows, &sub);
157 MoFEMFunctionReturn(0);
158 }
159
160 MoFEMErrorCode buildProjection(Vec control, IS columns,
161 const std::string &prefix,
162 PetscOptions options) {
163 MoFEMFunctionBegin;
164 if (!size)
165 MoFEMFunctionReturnHot(0);
166 auto ksp = createKSP(PetscObjectComm(getPetscObject(control)));
167 CHKERR KSPSetType(ksp, KSPPREONLY);
168 PC pc;
169 CHKERR KSPGetPC(ksp, &pc);
170 CHKERR PCSetType(pc, PCLU);
171 CHKERR PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);
172 CHKERR KSPSetTolerances(ksp, 1e-12, 1e-14, PETSC_DEFAULT, 1000);
173 CHKERR KSPSetOptionsPrefix(ksp, prefix.c_str());
174 CHKERR PetscObjectSetOptions(getPetscObject(ksp.get()), options);
175 CHKERR KSPSetFromOptions(ksp);
176 inverse = createCCTInverseKSP(jacobian, ksp, prefix);
177 Vec vector;
178 CHKERR MatCreateVecs(jacobian, &vector, nullptr);
179 SmartPetscObj<Vec> column_vector(vector);
180 VecScatter scatter;
181 CHKERR VecScatterCreate(control, columns, column_vector, columns, &scatter);
182 SmartPetscObj<VecScatter> identity(scatter);
183 projection = createConstrainMatrixCtx(identity, jacobian, inverse);
184 CHKERR projection->initializeQorP(control);
185 MoFEMFunctionReturn(0);
186 }
187};
188
189struct TaoNullSpace {
190 // Borrowed to evaluate the original problem without an ownership cycle.
191 Tao parentTao = nullptr;
192 SmartPetscObj<DM> constraintDM;
193 SmartPetscObj<Vec> constraints;
194 SmartPetscObj<Mat> jacobian;
195 SmartPetscObj<Mat> preconditioner;
196 SmartPetscObj<Vec> gradient;
197 SmartPetscObj<Vec> multipliers;
198 SmartPetscObj<Tao> subsolver;
199 SmartPetscObj<Vec> subSolution;
200 SmartPetscObj<Vec> subGradient;
201 SmartPetscObj<Vec> baseSolution;
202 SmartPetscObj<Vec> baseGradient;
203 SmartPetscObj<Vec> projectedGradient;
204 SmartPetscObj<Vec> stationarity;
205 SmartPetscObj<Vec> stepDirection;
206 SmartPetscObj<Vec> constraintWork;
207 SmartPetscObj<Vec> constraintWeights;
208 SmartPetscObj<Vec> bindingMultipliers;
209 SmartPetscObj<IS> allColumns;
210 NullSpaceRows binding;
211 NullSpaceRows restoration;
212 PetscReal alphaJ = 1;
213 PetscReal alphaC = 1;
214 PetscReal activeTolerance = 1e-6;
215 PetscReal multiplierTolerance = 1e-10;
216 PetscReal minimumStep = 1e-14;
217 PetscBool monitorKkt = PETSC_FALSE;
218 bool epochValid = false;
219 bool subproblemValid = false;
220 TaoLineSearch activeLineSearch = nullptr;
221 PetscErrorCode (*applyLineSearch)(TaoLineSearch, Vec, PetscReal *, Vec,
222 Vec) = nullptr;
223 PetscErrorCode (*destroyShell)(Tao) = nullptr;
224
225 /** Allocate vectors with the application's distributed layouts. */
226 MoFEMErrorCode setUp(Tao tao) {
227 MoFEMFunctionBegin;
228 if (tao->eq_constrained || tao->bounded || tao->ineq_doublesided ||
229 tao->ops->computebounds || tao->ops->computeconstraints)
230 SETERRQ(PetscObjectComm(getPetscObject(tao)), PETSC_ERR_SUP,
231 "TAONULLSPACE currently requires constraints in c_i >= 0 form; "
232 "separate equality constraints and bounds are unsupported");
233 if (!constraints || !jacobian)
234 SETERRQ(PetscObjectComm(getPetscObject(tao)), PETSC_ERR_ARG_WRONGSTATE,
235 "Call TaoNullSpaceSetConstraints before TaoSetUp");
236 Vec solution = nullptr;
237 CHKERR TaoGetSolution(tao, &solution);
238 gradient = vectorDuplicate(solution);
239 multipliers = vectorDuplicate(constraints);
240 CHKERR VecZeroEntries(multipliers);
241 subSolution = vectorDuplicate(solution);
242 subGradient = vectorDuplicate(solution);
243 baseSolution = vectorDuplicate(solution);
244 baseGradient = vectorDuplicate(solution);
245 projectedGradient = vectorDuplicate(solution);
246 stationarity = vectorDuplicate(solution);
247 stepDirection = vectorDuplicate(solution);
248 constraintWork = vectorDuplicate(constraints);
249 constraintWeights = vectorDuplicate(constraints);
250 CHKERR ::TaoSetGradient(tao, stationarity, nullptr, nullptr);
251 CHKERR VecDestroy(&tao->stepdirection);
252 tao->stepdirection = stepDirection;
253 CHKERR PetscObjectReference(getPetscObject(tao->stepdirection));
254 CHKERR VecZeroEntries(stepDirection);
255 PetscInt first, last;
256 CHKERR MatGetOwnershipRangeColumn(jacobian, &first, &last);
257 IS columns;
258 CHKERR ISCreateStride(PetscObjectComm(getPetscObject(tao)), last - first,
259 first, 1, &columns);
260 allColumns = SmartPetscObj<IS>(columns);
261 CHKERR VecCopy(solution, subSolution);
262 CHKERR TaoSetSolution(subsolver, subSolution);
263 CHKERR ::TaoSetObjective(
264 subsolver,
265 [](Tao, Vec x, PetscReal *objective, void *ctx) {
266 return static_cast<TaoNullSpace *>(ctx)->evaluateSubObjective(
267 x, objective);
268 },
269 this);
270 CHKERR ::TaoSetGradient(
271 subsolver, subGradient,
272 [](Tao, Vec x, Vec g, void *ctx) {
273 return static_cast<TaoNullSpace *>(ctx)->evaluateSubGradient(x, g);
274 },
275 this);
276 CHKERR ::TaoSetObjectiveAndGradient(
277 subsolver, subGradient,
278 [](Tao, Vec x, PetscReal *objective, Vec g, void *ctx) {
279 return static_cast<TaoNullSpace *>(ctx)->evaluateSubProblem(
280 x, objective, g);
281 },
282 this);
283 MoFEMFunctionReturn(0);
284 }
285
286 /** Evaluate values, and derivatives only when a gradient is requested. */
287 MoFEMErrorCode evaluate(Vec solution, PetscReal &objective, Vec g) {
288 MoFEMFunctionBegin;
289 if (g)
290 CHKERR TaoComputeObjectiveAndGradient(parentTao, solution, &objective, g);
291 else
292 CHKERR TaoComputeObjective(parentTao, solution, &objective);
293 // PETSC_INFINITY is a finite sentinel, also used by application callbacks
294 // when the equilibrium solve cannot evaluate a trial.
295 if (!std::isfinite(objective) || std::abs(objective) >= PETSC_INFINITY)
296 MoFEMFunctionReturnHot(0);
297 CHKERR TaoComputeInequalityConstraints(parentTao, solution, constraints);
298 PetscReal constraint_norm;
299 CHKERR VecNorm(constraints, NORM_2, &constraint_norm);
300 if (!std::isfinite(constraint_norm)) {
301 objective = PETSC_INFINITY;
302 MoFEMFunctionReturnHot(0);
303 }
304 if (g) {
305 PetscReal gradient_norm;
306 CHKERR VecNorm(g, NORM_2, &gradient_norm);
307 if (!std::isfinite(gradient_norm)) {
308 objective = PETSC_INFINITY;
309 MoFEMFunctionReturnHot(0);
310 }
311 CHKERR TaoComputeJacobianInequality(parentTao, solution, jacobian,
312 preconditioner);
313 }
314 MoFEMFunctionReturn(0);
315 }
316
317 /** NNLS selects binding rows without assuming every near-active row binds. */
318 MoFEMErrorCode selectMultipliers(NullSpaceRows &candidates) {
319 MoFEMFunctionBegin;
320 if (candidates.size) {
321 struct DualProblem {
322 Mat jacobian;
323 Vec gradient;
324 SmartPetscObj<Vec> residual;
325 } dual_problem{candidates.jacobian, baseGradient,
326 vectorDuplicate(baseGradient)};
327 auto dual = createTao(PetscObjectComm(getPetscObject(parentTao)));
328 auto lambda = vectorDuplicate(candidates.values);
329 auto dual_gradient = vectorDuplicate(lambda);
330 auto lower = vectorDuplicate(lambda);
331 auto upper = vectorDuplicate(lambda);
332 CHKERR candidates.extract(multipliers, lambda);
333 CHKERR VecZeroEntries(lower);
334 CHKERR VecSet(upper, PETSC_INFINITY);
335 CHKERR TaoSetType(dual, TAOBLMVM);
336 CHKERR TaoSetSolution(dual, lambda);
337 CHKERR TaoSetVariableBounds(dual, lower, upper);
338 CHKERR ::TaoSetObjectiveAndGradient(
339 dual, dual_gradient,
340 [](Tao, Vec mu, PetscReal *f, Vec g, void *ctx) {
341 MoFEMFunctionBegin;
342 auto *problem = static_cast<DualProblem *>(ctx);
343 CHKERR MatMultTranspose(problem->jacobian, mu, problem->residual);
344 CHKERR VecAXPY(problem->residual, -1, problem->gradient);
345 PetscScalar squared;
346 CHKERR VecDot(problem->residual, problem->residual, &squared);
347 *f = 0.5 * PetscRealPart(squared);
348 CHKERR MatMult(problem->jacobian, problem->residual, g);
349 MoFEMFunctionReturn(0);
350 },
351 &dual_problem);
352 const char *prefix = nullptr;
353 PetscOptions options = nullptr;
354 CHKERR TaoGetOptionsPrefix(parentTao, &prefix);
355 CHKERR PetscObjectGetOptions(getPetscObject(parentTao), &options);
356 CHKERR PetscObjectSetOptions(getPetscObject(dual.get()), options);
357 CHKERR TaoSetOptionsPrefix(dual, prefix);
358 CHKERR TaoAppendOptionsPrefix(dual, "tao_nullspace_dual_");
359 CHKERR TaoSetTolerances(dual, 1e-12, 0, 0);
360 CHKERR TaoSetMaximumIterations(dual, 500);
361 CHKERR TaoSetFromOptions(dual);
362 PetscBool diagonal_nnls;
363 CHKERR PetscObjectTypeCompare(getPetscObject(dual.get()), TAODIAGONALNNLS,
364 &diagonal_nnls);
365 if (diagonal_nnls)
366 CHKERR TaoDiagonalNNLSSetProblem(dual, candidates.jacobian,
367 baseGradient);
368 CHKERR TaoSolve(dual);
369 TaoConvergedReason reason;
370 PetscInt iterations;
371 PetscReal residual_norm;
372 CHKERR TaoGetSolutionStatus(dual, &iterations, nullptr, &residual_norm,
373 nullptr, nullptr, &reason);
374 if (reason == TAO_DIVERGED_MAXITS) {
375 MOFEM_LOG("WORLD", Sev::warning)
376 << "Null-space multiplier selection reached its iteration limit ("
377 << iterations << "), residual " << residual_norm
378 << "; continuing with the current multipliers";
379 } else if (reason <= 0)
380 SETERRQ(PetscObjectComm(getPetscObject(parentTao)),
381 PETSC_ERR_NOT_CONVERGED,
382 "Null-space multiplier selection failed: %s",
383 TaoConvergedReasons[reason]);
384 CHKERR VecZeroEntries(multipliers);
385 Vec sub;
386 CHKERR VecGetSubVector(multipliers, candidates.rows, &sub);
387 CHKERR VecCopy(lambda, sub);
388 CHKERR VecRestoreSubVector(multipliers, candidates.rows, &sub);
389 } else {
390 CHKERR VecZeroEntries(multipliers);
391 }
392 MoFEMFunctionReturn(0);
393 }
394
395 /** Freeze ordered active sets, copied Jacobians and inverse actions. */
396 MoFEMErrorCode buildEpoch(Vec solution) {
397 MoFEMFunctionBegin;
398 epochValid = false;
399 CHKERR VecCopy(solution, baseSolution);
400 CHKERR VecCopy(gradient, baseGradient);
401 PetscInt first, last;
402 CHKERR VecGetOwnershipRange(constraints, &first, &last);
403 const PetscScalar *c;
404 CHKERR VecGetArrayRead(constraints, &c);
405 std::vector<PetscInt> candidates_indices;
406 for (PetscInt i = 0; i != last - first; ++i) {
407 PetscInt count;
408 const PetscScalar *values;
409 CHKERR MatGetRow(jacobian, first + i, &count, nullptr, &values);
410 PetscReal row_norm = 0;
411 for (PetscInt j = 0; j != count; ++j)
412 row_norm = std::hypot(row_norm, PetscAbsScalar(values[j]));
413 CHKERR MatRestoreRow(jacobian, first + i, &count, nullptr, &values);
414 // Feppon et al., Section 4.1: epsilon_i = h ||grad c_i|| makes
415 // candidate selection invariant to positive row scaling.
416 if (PetscRealPart(c[i]) <= activeTolerance * row_norm)
417 candidates_indices.push_back(first + i);
418 }
419 CHKERR VecRestoreArrayRead(constraints, &c);
420 NullSpaceRows candidates;
421 CHKERR candidates.select(jacobian, allColumns, candidates_indices);
422 CHKERR selectMultipliers(candidates);
423 const PetscScalar *mu;
424 CHKERR VecGetArrayRead(constraints, &c);
425 CHKERR VecGetArrayRead(multipliers, &mu);
426 std::vector<PetscInt> binding_indices;
427 std::vector<PetscInt> restoration_indices;
428 for (PetscInt i = 0; i != last - first; ++i) {
429 const bool binds = PetscRealPart(mu[i]) > multiplierTolerance;
430 if (binds)
431 binding_indices.push_back(first + i);
432 if (binds || PetscRealPart(c[i]) <= 0)
433 restoration_indices.push_back(first + i);
434 }
435 CHKERR VecRestoreArrayRead(multipliers, &mu);
436 CHKERR VecRestoreArrayRead(constraints, &c);
437 CHKERR binding.select(jacobian, allColumns, binding_indices);
438 CHKERR restoration.select(jacobian, allColumns, restoration_indices);
439 const char *prefix = nullptr;
440 PetscOptions options = nullptr;
441 CHKERR TaoGetOptionsPrefix(parentTao, &prefix);
442 CHKERR PetscObjectGetOptions(getPetscObject(parentTao), &options);
443 const std::string base_prefix = prefix ? prefix : "";
444 CHKERR binding.buildProjection(baseGradient, allColumns,
445 base_prefix + "tao_nullspace_binding_",
446 options);
447 CHKERR restoration.buildProjection(
448 baseGradient, allColumns, base_prefix + "tao_nullspace_restoration_",
449 options);
450 if (binding.size) {
451 bindingMultipliers = vectorDuplicate(binding.values);
452 CHKERR MatMult(binding.jacobian, baseGradient, binding.work);
453 CHKERR MatMult(binding.inverse, binding.work, bindingMultipliers);
454 CHKERR binding.projection->applyQ(baseGradient, projectedGradient);
455 } else {
456 bindingMultipliers.reset();
457 CHKERR VecCopy(baseGradient, projectedGradient);
458 }
459 epochValid = true;
460 MoFEMFunctionReturn(0);
461 }
462
463 /** Scalar merit and its exact derivative for one frozen active-set epoch.
464 * Phi = alphaJ*(f - lambda_B^T c_B) + alphaC/2*c_R^T S_R*c_R.
465 * Live trial Jacobians never overwrite the accepted-iterate inverse actions.
466 */
467 MoFEMErrorCode evaluateSubProblem(Vec x, PetscReal *objective, Vec g) {
468 MoFEMFunctionBegin;
469 subproblemValid = false;
470 if (!epochValid)
471 SETERRQ(PetscObjectComm(getPetscObject(parentTao)),
472 PETSC_ERR_ARG_WRONGSTATE,
473 "The null-space subproblem needs an accepted-iterate active set");
474 PetscReal f;
475 CHKERR evaluate(x, f, g ? gradient.get() : nullptr);
476 if (!std::isfinite(f) || std::abs(f) >= PETSC_INFINITY) {
477 CHKERR rejectSubProblem(objective, g);
478 MoFEMFunctionReturnHot(0);
479 }
480 *objective = alphaJ * f;
481 if (g)
482 CHKERR VecZeroEntries(constraintWeights);
483 if (binding.size) {
484 CHKERR binding.extract(constraints, binding.values);
485 PetscScalar value;
486 CHKERR VecDot(bindingMultipliers, binding.values, &value);
487 *objective -= alphaJ * PetscRealPart(value);
488 if (g) {
489 Vec sub;
490 CHKERR VecGetSubVector(constraintWeights, binding.rows, &sub);
491 CHKERR VecCopy(bindingMultipliers, sub);
492 CHKERR VecScale(sub, -alphaJ);
493 CHKERR VecRestoreSubVector(constraintWeights, binding.rows, &sub);
494 }
495 }
496 if (restoration.size) {
497 CHKERR restoration.extract(constraints, restoration.values);
498 CHKERR MatMult(restoration.inverse, restoration.values, restoration.work);
499 PetscScalar value;
500 CHKERR VecDot(restoration.values, restoration.work, &value);
501 *objective += 0.5 * alphaC * PetscRealPart(value);
502 if (g) {
503 Vec sub;
504 CHKERR VecGetSubVector(constraintWeights, restoration.rows, &sub);
505 CHKERR VecAXPY(sub, alphaC, restoration.work);
506 CHKERR VecRestoreSubVector(constraintWeights, restoration.rows, &sub);
507 }
508 }
509 if (g) {
510 CHKERR MatMultTranspose(jacobian, constraintWeights, g);
511 CHKERR VecAXPY(g, alphaJ, gradient);
512 PetscReal norm;
513 CHKERR VecNorm(g, NORM_2, &norm);
514 if (!std::isfinite(norm)) {
515 CHKERR rejectSubProblem(objective, g);
516 MoFEMFunctionReturnHot(0);
517 }
518 }
519 if (!std::isfinite(*objective) || std::abs(*objective) >= PETSC_INFINITY) {
520 CHKERR rejectSubProblem(objective, g);
521 MoFEMFunctionReturnHot(0);
522 }
523 subproblemValid = true;
524 MoFEMFunctionReturn(0);
525 }
526
527 /** A rejected trial has no usable objective/derivative pair. */
528 MoFEMErrorCode rejectSubProblem(PetscReal *objective, Vec g) {
529 MoFEMFunctionBegin;
530 *objective = PETSC_INFINITY;
531 if (g)
532 CHKERR VecZeroEntries(g);
533 if (activeLineSearch) {
534 // PETSc More-Thuente checks the current step bounds before calling its
535 // cubic interpolation. With rejected f and zero g, collapsing the bounds
536 // halts there; the apply adapter then backtracks using actual evaluations.
537 activeLineSearch->stepmin = activeLineSearch->step;
538 activeLineSearch->stepmax = activeLineSearch->step;
539 }
540 MoFEMFunctionReturn(0);
541 }
542
543 /** Preserve native More-Thuente on valid trials; backtrack after rejection.
544 * A failed physical evaluation can prevent a strong-Wolfe step.
545 * A valid Armijo step lets the outer iteration continue.
546 */
547 MoFEMErrorCode applyMoreThuente(TaoLineSearch ls, Vec x, PetscReal *f, Vec g,
548 Vec direction) {
549 MoFEMFunctionBegin;
550 auto start = vectorDuplicate(x);
551 auto start_gradient = vectorDuplicate(g);
552 CHKERR VecCopy(x, start);
553 CHKERR VecCopy(g, start_gradient);
554 const auto initial_value = *f;
555 const auto step_min = ls->stepmin;
556 const auto step_max = ls->stepmax;
557 PetscScalar slope;
558 CHKERR VecDot(g, direction, &slope);
559 subproblemValid = true;
560 activeLineSearch = ls;
561 const auto error = applyLineSearch(ls, x, f, g, direction);
562 activeLineSearch = nullptr;
563 ls->stepmin = step_min;
564 ls->stepmax = step_max;
565 CHKERR error;
566 if (subproblemValid)
567 MoFEMFunctionReturnHot(0);
568
569 CHKERR PetscInfo(ls, "Null-space trial rejected; backtracking without "
570 "interpolating rejected objective/gradient data\n");
571 ls->step = std::min(0.5 * ls->step, step_max);
572 while (ls->nfeval + ls->nfgeval < ls->max_funcs && ls->step >= step_min &&
573 ls->step > 0) {
574 CHKERR VecWAXPY(x, ls->step, direction, start);
575 CHKERR TaoLineSearchComputeObjective(ls, x, f);
576 CHKERR TaoLineSearchMonitor(ls, ls->nfeval + ls->nfgeval, *f, ls->step);
577 if (subproblemValid &&
578 *f <= initial_value + ls->ftol * ls->step * PetscRealPart(slope)) {
579 CHKERR TaoLineSearchComputeObjectiveAndGradient(ls, x, f, g);
580 if (subproblemValid) {
581 // This is sufficient-decrease acceptance at a domain boundary; it
582 // does not claim the native More-Thuente curvature test succeeded.
583 ls->reason = TAOLINESEARCH_SUCCESS_USER;
584 CHKERR PetscInfo(ls,
585 "Null-space domain safeguard accepted a valid "
586 "sufficient-decrease step %g\n",
587 static_cast<double>(ls->step));
588 MoFEMFunctionReturnHot(0);
589 }
590 }
591 ls->step *= 0.5;
592 }
593 ls->reason = ls->nfeval + ls->nfgeval >= ls->max_funcs
594 ? TAOLINESEARCH_HALTED_MAXFCN
595 : TAOLINESEARCH_HALTED_LOWERBOUND;
596 CHKERR VecCopy(start, x);
597 CHKERR VecCopy(start_gradient, g);
598 *f = initial_value;
599 ls->step = 0;
600 MoFEMFunctionReturn(0);
601 }
602
603 /** Install the apply adapter only while the owned inner TAO is solving. */
604 MoFEMErrorCode solveSubProblem() {
605 MoFEMFunctionBegin;
606 TaoLineSearch ls = nullptr;
607 CHKERR TaoGetLineSearch(subsolver, &ls);
608 PetscBool more_thuente = PETSC_FALSE;
609 if (ls)
610 CHKERR PetscObjectTypeCompare(getPetscObject(ls), TAOLINESEARCHMT,
611 &more_thuente);
612 if (!more_thuente) {
613 CHKERR TaoSolve(subsolver);
614 MoFEMFunctionReturnHot(0);
615 }
616 CHKERR TaoLineSearchUseTaoRoutines(ls, subsolver);
617 applyLineSearch = ls->ops->apply;
618 ls->ops->apply = [](TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s) {
619 return static_cast<TaoNullSpace *>(ls->tao->user_objgradP)
620 ->applyMoreThuente(ls, x, f, g, s);
621 };
622 const auto error = TaoSolve(subsolver);
623 // Also detach on errors, so a retained subsolver has no adapter borrowing
624 // the parent context once this solve returns or the parent is destroyed.
625 ls->ops->apply = applyLineSearch;
626 applyLineSearch = nullptr;
627 activeLineSearch = nullptr;
628 CHKERR error;
629 MoFEMFunctionReturn(0);
630 }
631
632 MoFEMErrorCode evaluateSubObjective(Vec x, PetscReal *objective) {
633 MoFEMFunctionBegin;
634 CHKERR evaluateSubProblem(x, objective, nullptr);
635 MoFEMFunctionReturn(0);
636 }
637
638 MoFEMErrorCode evaluateSubGradient(Vec x, Vec g) {
639 MoFEMFunctionBegin;
640 PetscReal objective;
641 CHKERR evaluateSubProblem(x, &objective, g);
642 MoFEMFunctionReturn(0);
643 }
644
645 /** Solve the frozen inner problem, then re-evaluate the full KKT system. */
646 MoFEMErrorCode solve(Tao tao) {
647 MoFEMFunctionBegin;
648 CHKERR TaoSetConvergedReason(tao, TAO_CONTINUE_ITERATING);
649 Vec solution;
650 CHKERR TaoGetSolution(tao, &solution);
651 PetscReal step = 1;
652 for (PetscInt iteration = 0;; ++iteration) {
653 PetscReal objective;
654 CHKERR evaluate(solution, objective, gradient);
655 if (!std::isfinite(objective) || std::abs(objective) >= PETSC_INFINITY)
656 SETERRQ(PetscObjectComm(getPetscObject(tao)), PETSC_ERR_FP,
657 "Null-space accepted iterate has non-finite objective, "
658 "gradient or constraint data");
659 CHKERR buildEpoch(solution);
660 CHKERR MatMultTranspose(jacobian, multipliers, stationarity);
661 CHKERR VecAYPX(stationarity, -1, gradient);
662 PetscReal stationarity_norm, projected_norm, complementarity_norm,
663 constraint_norm;
664 CHKERR VecNorm(stationarity, NORM_2, &stationarity_norm);
665 CHKERR VecNorm(projectedGradient, NORM_2, &projected_norm);
666 CHKERR VecPointwiseMult(constraintWork, multipliers, constraints);
667 CHKERR VecNorm(constraintWork, NORM_2, &complementarity_norm);
668 CHKERR VecZeroEntries(constraintWork);
669 CHKERR VecPointwiseMin(constraintWork, constraintWork, constraints);
670 CHKERR VecNorm(constraintWork, NORM_2, &constraint_norm);
671 if (!std::isfinite(stationarity_norm) || !std::isfinite(projected_norm) ||
672 !std::isfinite(complementarity_norm) ||
673 !std::isfinite(constraint_norm))
674 SETERRQ(PetscObjectComm(getPetscObject(tao)), PETSC_ERR_FP,
675 "Null-space KKT residual is not finite");
676 const auto residual =
677 std::max({stationarity_norm, projected_norm, complementarity_norm});
678 CHKERR TaoSetIterationNumber(tao, iteration);
679 CHKERR TaoLogConvergenceHistory(tao, objective, residual, constraint_norm,
680 0);
681 CHKERR TaoMonitor(tao, iteration, objective, residual, constraint_norm,
682 step);
683 if (monitorKkt)
684 MOFEM_LOG("WORLD", Sev::inform)
685 << "Null-space KKT iteration " << iteration
686 << ": stationarity " << stationarity_norm
687 << ", projected gradient " << projected_norm
688 << ", complementarity " << complementarity_norm;
689 CHKERR (*tao->ops->convergencetest)(tao, tao->cnvP);
690 if (tao->reason == TAO_CONVERGED_STEPTOL) {
691 PetscReal gatol, catol;
692 CHKERR TaoGetTolerances(tao, &gatol, nullptr, nullptr);
693 CHKERR TaoGetConstraintTolerances(tao, &catol, nullptr);
694 if (residual > gatol || constraint_norm > catol)
695 CHKERR TaoSetConvergedReason(tao, TAO_DIVERGED_LS_FAILURE);
696 }
697 if (tao->reason != TAO_CONTINUE_ITERATING)
698 break;
699 PetscReal merit_before, merit_after;
700 CHKERR evaluateSubProblem(baseSolution, &merit_before, subGradient);
701 CHKERR VecCopy(baseSolution, subSolution);
702 CHKERR TaoSetRecycleHistory(subsolver, PETSC_FALSE);
703 CHKERR solveSubProblem();
704 TaoConvergedReason reason;
705 CHKERR TaoGetConvergedReason(subsolver, &reason);
706 CHKERR evaluateSubProblem(subSolution, &merit_after, subGradient);
707 CHKERR VecWAXPY(stepDirection, -1, baseSolution, subSolution);
708 CHKERR VecNorm(stepDirection, NORM_2, &step);
709 TaoLineSearch ls;
710 CHKERR TaoGetLineSearch(subsolver, &ls);
711 PetscBool algorithm_one = PETSC_FALSE;
712 if (ls)
713 CHKERR PetscObjectTypeCompare(getPetscObject(ls), nullSpaceLineSearch,
714 &algorithm_one);
715 const bool final_trial =
716 algorithm_one && ls->reason == TAOLINESEARCH_SUCCESS_USER;
717 if ((reason <= 0 && reason != TAO_DIVERGED_MAXITS) || !subproblemValid ||
718 !std::isfinite(step) || step <= minimumStep ||
719 (!final_trial &&
720 merit_after > merit_before + 8 * PETSC_MACHINE_EPSILON *
721 (1 + std::abs(merit_before)))) {
722 // Restore the application's live evaluation when rejecting the trial.
723 CHKERR evaluate(solution, objective, gradient);
724 CHKERR TaoSetConvergedReason(tao, TAO_DIVERGED_LS_FAILURE);
725 break;
726 }
727 CHKERR VecCopy(subSolution, solution);
728 }
729 MoFEMFunctionReturn(0);
730 }
731};
732
733/** Get the composed context: TaoShellGetContext only accepts type 'shell'. */
734MoFEMErrorCode getNullSpaceContext(Tao tao, TaoNullSpace *&ctx_ptr) {
735 MoFEMFunctionBegin;
736 PetscBool is_nullspace = PETSC_FALSE;
737 CHKERR PetscObjectTypeCompare(getPetscObject(tao), TAONULLSPACE,
738 &is_nullspace);
739 if (!is_nullspace)
740 SETERRQ(PetscObjectComm(getPetscObject(tao)), PETSC_ERR_ARG_WRONG,
741 "TAO type must be nullspace");
742 PetscContainer container = nullptr;
743 CHKERR PetscObjectQuery(getPetscObject(tao), nullSpaceContextKey,
744 reinterpret_cast<PetscObject *>(&container));
745 if (!container)
746 SETERRQ(PetscObjectComm(getPetscObject(tao)), PETSC_ERR_ARG_WRONGSTATE,
747 "TAONULLSPACE context is not attached");
748 void *ptr = nullptr;
749 CHKERR PetscContainerGetPointer(container, &ptr);
750 ctx_ptr = static_cast<TaoNullSpace *>(ptr);
751 MoFEMFunctionReturn(0);
752}
753
754/** Construct a TAOSHELL alias with an owned context and lifecycle hooks. */
755MoFEMErrorCode createTaoNullSpace(Tao tao) {
756 MoFEMFunctionBegin;
757 // TaoSetType assigns the registered name after this constructor returns.
758 // Clear the previous name so switching from 'shell' also constructs new data.
759 CHKERR PetscObjectChangeTypeName(getPetscObject(tao), nullptr);
760 CHKERR TaoSetType(tao, TAOSHELL);
761 const auto comm = PetscObjectComm(getPetscObject(tao));
762 auto ctx_ptr = std::make_unique<TaoNullSpace>();
763 ctx_ptr->parentTao = tao;
764 ctx_ptr->destroyShell = tao->ops->destroy;
765 ctx_ptr->subsolver = createTao(comm);
766 CHKERR TaoSetType(ctx_ptr->subsolver, TAOSHELL);
767 CHKERR TaoShellSetSolve(ctx_ptr->subsolver, solveNullSpaceFlowStep);
768 CHKERR TaoSetTolerances(ctx_ptr->subsolver, 1e-14, 0, 0);
769 CHKERR TaoSetMaximumIterations(ctx_ptr->subsolver, 1);
770 // TAOSHELL does not create a line search; the owning TAO destroys it.
771 CHKERR TaoLineSearchCreate(comm, &ctx_ptr->subsolver.get()->linesearch);
772 CHKERR TaoLineSearchSetType(ctx_ptr->subsolver.get()->linesearch,
773 nullSpaceLineSearch);
774
775 PetscContainer container = nullptr;
776 CHKERR PetscContainerCreate(comm, &container);
777 CHKERR PetscContainerSetPointer(container, ctx_ptr.get());
778 CHKERR PetscContainerSetUserDestroy(container, [](void *ptr) {
779 MoFEMFunctionBegin;
780 auto *ctx_ptr = static_cast<TaoNullSpace *>(ptr);
781 // A caller may retain the inner TAO beyond the parent. PETSc's public
782 // setters ignore null callbacks, so detach our borrowed context here.
783 auto inner = ctx_ptr->subsolver.get();
784 inner->ops->computeobjective = nullptr;
785 inner->user_objP = nullptr;
786 inner->ops->computegradient = nullptr;
787 inner->user_gradP = nullptr;
788 inner->ops->computeobjectiveandgradient = nullptr;
789 inner->user_objgradP = nullptr;
790 delete ctx_ptr;
791 MoFEMFunctionReturn(0);
792 });
793 ctx_ptr.release();
794 const auto compose_error = PetscObjectCompose(
795 getPetscObject(tao), nullSpaceContextKey, getPetscObject(container));
796 CHKERR PetscContainerDestroy(&container);
797 CHKERR compose_error;
798
799 tao->ops->destroy = [](Tao tao) {
800 MoFEMFunctionBegin;
801 TaoNullSpace *ctx_ptr = nullptr;
802 CHKERR getNullSpaceContext(tao, ctx_ptr);
803 const auto destroy_shell = ctx_ptr->destroyShell;
804 // Release immediately on TaoSetType too, not only on TaoDestroy.
805 CHKERR PetscObjectCompose(getPetscObject(tao), nullSpaceContextKey,
806 nullptr);
807 CHKERR destroy_shell(tao);
808 MoFEMFunctionReturn(0);
809 };
810 tao->ops->setup = [](Tao tao) {
811 MoFEMFunctionBegin;
812 TaoNullSpace *ctx_ptr = nullptr;
813 CHKERR getNullSpaceContext(tao, ctx_ptr);
814 CHKERR ctx_ptr->setUp(tao);
815 MoFEMFunctionReturn(0);
816 };
817 tao->ops->setfromoptions = [](Tao tao, PetscOptionItems *PetscOptionsObject) {
818 MoFEMFunctionBegin;
819 TaoNullSpace *ctx_ptr = nullptr;
820 CHKERR getNullSpaceContext(tao, ctx_ptr);
821 PetscOptionsHeadBegin(PetscOptionsObject, "Null-space TAO shell options");
822 CHKERR PetscOptionsBool(
823 "-tao_nullspace_monitor_kkt", "Print individual KKT residual norms",
824 nullptr, ctx_ptr->monitorKkt, &ctx_ptr->monitorKkt, nullptr);
825 CHKERR PetscOptionsReal("-tao_nullspace_alpha_j", "Objective descent scale",
826 nullptr, ctx_ptr->alphaJ, &ctx_ptr->alphaJ,
827 nullptr);
828 CHKERR PetscOptionsReal("-tao_nullspace_alpha_c",
829 "Constraint restoration scale", nullptr,
830 ctx_ptr->alphaC, &ctx_ptr->alphaC, nullptr);
831 CHKERR PetscOptionsReal(
832 "-tao_nullspace_active_tol", "Activation distance h in control space",
833 nullptr, ctx_ptr->activeTolerance, &ctx_ptr->activeTolerance, nullptr);
834 CHKERR PetscOptionsReal("-tao_nullspace_multiplier_tol",
835 "Binding multiplier tolerance", nullptr,
836 ctx_ptr->multiplierTolerance,
837 &ctx_ptr->multiplierTolerance, nullptr);
838 CHKERR PetscOptionsReal(
839 "-tao_nullspace_min_step", "Minimum accepted step norm", nullptr,
840 ctx_ptr->minimumStep, &ctx_ptr->minimumStep, nullptr);
841 PetscOptionsHeadEnd();
842 if (!(ctx_ptr->alphaJ > 0) || !(ctx_ptr->alphaC > 0) ||
843 !std::isfinite(ctx_ptr->alphaJ) || !std::isfinite(ctx_ptr->alphaC) ||
844 !std::isfinite(ctx_ptr->activeTolerance) ||
845 ctx_ptr->activeTolerance < 0 ||
846 !(ctx_ptr->multiplierTolerance > 0) ||
847 !(ctx_ptr->minimumStep > 0))
848 SETERRQ(PetscObjectComm(getPetscObject(tao)), PETSC_ERR_ARG_OUTOFRANGE,
849 "Require finite positive null-space scales, finite active_tol "
850 ">= 0, and positive multiplier_tol/min_step");
851 const char *prefix = nullptr;
852 PetscOptions options = nullptr;
853 CHKERR TaoGetOptionsPrefix(tao, &prefix);
854 CHKERR PetscObjectGetOptions(getPetscObject(tao), &options);
855 CHKERR PetscObjectSetOptions(getPetscObject(ctx_ptr->subsolver.get()),
856 options);
857 CHKERR TaoSetOptionsPrefix(ctx_ptr->subsolver, prefix);
858 CHKERR TaoAppendOptionsPrefix(ctx_ptr->subsolver, "tao_nullspace_");
859 if (ctx_ptr->subsolver.get()->linesearch)
860 CHKERR PetscObjectSetOptions(
861 getPetscObject(ctx_ptr->subsolver.get()->linesearch), options);
862 CHKERR TaoSetFromOptions(ctx_ptr->subsolver);
863 PetscBool shell;
864 CHKERR PetscObjectTypeCompare(getPetscObject(ctx_ptr->subsolver.get()),
865 TAOSHELL, &shell);
866 if (shell) {
867 CHKERR TaoShellSetSolve(ctx_ptr->subsolver, solveNullSpaceFlowStep);
868 // Switching an iterative inner TAO back to shell destroys its search.
869 auto inner = ctx_ptr->subsolver.get();
870 if (!inner->linesearch) {
871 CHKERR TaoLineSearchCreate(PetscObjectComm(getPetscObject(inner)),
872 &inner->linesearch);
873 CHKERR TaoLineSearchSetType(inner->linesearch, nullSpaceLineSearch);
874 const char *inner_prefix;
875 CHKERR TaoGetOptionsPrefix(inner, &inner_prefix);
876 CHKERR TaoLineSearchSetOptionsPrefix(inner->linesearch, inner_prefix);
877 CHKERR PetscObjectSetOptions(getPetscObject(inner->linesearch), options);
878 CHKERR TaoLineSearchSetFromOptions(inner->linesearch);
879 }
880 }
881 MoFEMFunctionReturn(0);
882 };
883 tao->ops->view = [](Tao tao, PetscViewer viewer) {
884 MoFEMFunctionBegin;
885 TaoNullSpace *ctx_ptr = nullptr;
886 CHKERR getNullSpaceContext(tao, ctx_ptr);
887 PetscBool is_ascii = PETSC_FALSE;
888 CHKERR PetscObjectTypeCompare(getPetscObject(viewer), PETSCVIEWERASCII,
889 &is_ascii);
890 if (is_ascii) {
891 CHKERR PetscViewerASCIIPrintf(
892 viewer,
893 " Null-space gradient flow (default: Algorithm 1 Euler step)\n"
894 " Objective scale %g, constraint scale %g\n"
895 " Activation distance h %g (epsilon_i = h ||J_i||_2)\n",
896 static_cast<double>(ctx_ptr->alphaJ),
897 static_cast<double>(ctx_ptr->alphaC),
898 static_cast<double>(ctx_ptr->activeTolerance));
899 CHKERR TaoView(ctx_ptr->subsolver, viewer);
900 }
901 MoFEMFunctionReturn(0);
902 };
903 CHKERR TaoShellSetSolve(tao, [](Tao tao) {
904 MoFEMFunctionBegin;
905 TaoNullSpace *ctx_ptr = nullptr;
906 CHKERR getNullSpaceContext(tao, ctx_ptr);
907 CHKERR ctx_ptr->solve(tao);
908 MoFEMFunctionReturn(0);
909 });
910 MoFEMFunctionReturn(0);
911}
912
913} // namespace
914
915MoFEMErrorCode TaoRegisterNullSpace() {
916 MoFEMFunctionBegin;
917 CHKERR TaoRegisterDiagonalNNLS();
918 CHKERR TaoLineSearchRegister(nullSpaceLineSearch, createNullSpaceLineSearch);
919 CHKERR TaoRegister(TAONULLSPACE, createTaoNullSpace);
920 MoFEMFunctionReturn(0);
921}
922
923MoFEMErrorCode TaoNullSpaceSetConstraints(Tao tao, DM constraint_dm,
924 Vec constraints, Mat jacobian,
925 Mat preconditioner) {
926 MoFEMFunctionBegin;
927 TaoNullSpace *ctx_ptr = nullptr;
928 CHKERR getNullSpaceContext(tao, ctx_ptr);
929 if (tao->setupcalled)
930 SETERRQ(PetscObjectComm(getPetscObject(tao)), PETSC_ERR_ARG_WRONGSTATE,
931 "Set null-space constraints before TaoSetUp");
932 ctx_ptr->constraintDM = SmartPetscObj<DM>(constraint_dm, true);
933 ctx_ptr->constraints = SmartPetscObj<Vec>(constraints, true);
934 ctx_ptr->jacobian = SmartPetscObj<Mat>(jacobian, true);
935 ctx_ptr->preconditioner =
936 SmartPetscObj<Mat>(preconditioner ? preconditioner : jacobian, true);
937 MoFEMFunctionReturn(0);
938}
939
940MoFEMErrorCode TaoNullSpaceGetSubsolver(Tao tao, Tao *subsolver) {
941 MoFEMFunctionBegin;
942 TaoNullSpace *ctx_ptr = nullptr;
943 CHKERR getNullSpaceContext(tao, ctx_ptr);
944 *subsolver = ctx_ptr->subsolver;
945 MoFEMFunctionReturn(0);
946}
947
948MoFEMErrorCode TaoNullSpaceGetMultipliers(Tao tao, Vec *multipliers) {
949 MoFEMFunctionBegin;
950 TaoNullSpace *ctx_ptr = nullptr;
951 CHKERR getNullSpaceContext(tao, ctx_ptr);
952 if (!ctx_ptr->multipliers)
953 SETERRQ(PetscObjectComm(getPetscObject(tao)), PETSC_ERR_ARG_WRONGSTATE,
954 "Call TaoSetUp before requesting null-space multipliers");
955 *multipliers = ctx_ptr->multipliers;
956 MoFEMFunctionReturn(0);
957}
958
959} // namespace MoFEM
comm
MPI_Comm comm
Definition ConstrainMatrixCtx.cpp:64
ksp
SmartPetscObj< KSP > ksp
Definition ConstrainMatrixCtx.cpp:85
work
SmartPetscObj< Vec > work
Definition TaoDiagonalNNLS.cpp:22
residual
SmartPetscObj< Vec > residual
Definition TaoDiagonalNNLS.cpp:19
parentTao
Tao parentTao
Definition TaoDiagonalNNLS.cpp:16
gradient
SmartPetscObj< Vec > gradient
Definition TaoDiagonalNNLS.cpp:20
jacobian
SmartPetscObj< Mat > jacobian
Definition TaoDiagonalNNLS.cpp:17
alphaJ
PetscReal alphaJ
Definition TaoNullSpace.cpp:212
values
SmartPetscObj< Vec > values
Definition TaoNullSpace.cpp:122
baseSolution
SmartPetscObj< Vec > baseSolution
Definition TaoNullSpace.cpp:201
multiplierTolerance
PetscReal multiplierTolerance
Definition TaoNullSpace.cpp:215
preconditioner
SmartPetscObj< Mat > preconditioner
Definition TaoNullSpace.cpp:195
projectedGradient
SmartPetscObj< Vec > projectedGradient
Definition TaoNullSpace.cpp:203
minimumStep
PetscReal minimumStep
Definition TaoNullSpace.cpp:216
baseGradient
SmartPetscObj< Vec > baseGradient
Definition TaoNullSpace.cpp:202
allColumns
SmartPetscObj< IS > allColumns
Definition TaoNullSpace.cpp:209
alphaC
PetscReal alphaC
Definition TaoNullSpace.cpp:213
constraintWork
SmartPetscObj< Vec > constraintWork
Definition TaoNullSpace.cpp:206
constraintDM
SmartPetscObj< DM > constraintDM
Definition TaoNullSpace.cpp:192
stepDirection
SmartPetscObj< Vec > stepDirection
Definition TaoNullSpace.cpp:205
constraintWeights
SmartPetscObj< Vec > constraintWeights
Definition TaoNullSpace.cpp:207
monitorKkt
PetscBool monitorKkt
Definition TaoNullSpace.cpp:217
projection
boost::shared_ptr< ConstrainMatrixCtx > projection
Definition TaoNullSpace.cpp:124
subsolver
SmartPetscObj< Tao > subsolver
Definition TaoNullSpace.cpp:198
size
PetscInt size
Definition TaoNullSpace.cpp:125
constraints
SmartPetscObj< Vec > constraints
Definition TaoNullSpace.cpp:193
subproblemValid
bool subproblemValid
Definition TaoNullSpace.cpp:219
stationarity
SmartPetscObj< Vec > stationarity
Definition TaoNullSpace.cpp:204
inverse
SmartPetscObj< Mat > inverse
Definition TaoNullSpace.cpp:121
activeLineSearch
TaoLineSearch activeLineSearch
Definition TaoNullSpace.cpp:220
bindingMultipliers
SmartPetscObj< Vec > bindingMultipliers
Definition TaoNullSpace.cpp:208
subGradient
SmartPetscObj< Vec > subGradient
Definition TaoNullSpace.cpp:200
epochValid
bool epochValid
Definition TaoNullSpace.cpp:218
multipliers
SmartPetscObj< Vec > multipliers
Definition TaoNullSpace.cpp:197
activeTolerance
PetscReal activeTolerance
Definition TaoNullSpace.cpp:214
binding
NullSpaceRows binding
Definition TaoNullSpace.cpp:210
applyLineSearch
PetscErrorCode(* applyLineSearch)(TaoLineSearch, Vec, PetscReal *, Vec, Vec)
Definition TaoNullSpace.cpp:221
restoration
NullSpaceRows restoration
Definition TaoNullSpace.cpp:211
subSolution
SmartPetscObj< Vec > subSolution
Definition TaoNullSpace.cpp:199
rows
SmartPetscObj< IS > rows
Definition TaoNullSpace.cpp:119
TaoNullSpace.hpp
TAO shell interface for null-space constrained optimization.
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
mu
double mu
Definition dynamic_first_order_con_law.cpp:97
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MoFEM::TAONULLSPACE
constexpr char TAONULLSPACE[]
Definition TaoNullSpace.hpp:13
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
lambda
static double lambda
Definition incompressible_elasticity.cpp:181
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:40
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition MatrixFunction.hpp:64
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition PetscSmartObj.hpp:269
MoFEM::TaoNullSpaceGetMultipliers
MoFEMErrorCode TaoNullSpaceGetMultipliers(Tao tao, Vec *multipliers)
Definition TaoNullSpace.cpp:948
MoFEM::TaoDiagonalNNLSSetProblem
MoFEMErrorCode TaoDiagonalNNLSSetProblem(Tao tao, Mat jacobian, Vec target)
Definition TaoDiagonalNNLS.cpp:289
MoFEM::TaoNullSpaceSetConstraints
MoFEMErrorCode TaoNullSpaceSetConstraints(Tao tao, DM constraint_dm, Vec constraints, Mat jacobian, Mat preconditioner)
Definition TaoNullSpace.cpp:923
MoFEM::createCCTInverseKSP
SmartPetscObj< Mat > createCCTInverseKSP(Mat constrain_matrix, KSP ksp=nullptr, const std::string &options_prefix="constraint_matrix_")
Definition ConstrainMatrixCtx.cpp:1023
MoFEM::TaoRegisterDiagonalNNLS
MoFEMErrorCode TaoRegisterDiagonalNNLS()
Definition TaoDiagonalNNLS.cpp:283
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
MoFEM::solve
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
Definition Schur.cpp:2951
MoFEM::TAODIAGONALNNLS
constexpr char TAODIAGONALNNLS[]
Definition TaoDiagonalNNLS.hpp:11
MoFEM::createConstrainMatrixCtx
boost::shared_ptr< ConstrainMatrixCtx > createConstrainMatrixCtx(DM x_dm, DM y_dm, Mat constrain_matrix, Mat inverse_cct)
Definition ConstrainMatrixCtx.cpp:977
MoFEM::getPetscObject
PetscObject getPetscObject(T obj)
Definition PetscSmartObj.hpp:10
MoFEM::TaoNullSpaceGetSubsolver
MoFEMErrorCode TaoNullSpaceGetSubsolver(Tao tao, Tao *subsolver)
Definition TaoNullSpace.cpp:940
MoFEM::TaoRegisterNullSpace
MoFEMErrorCode TaoRegisterNullSpace()
Definition TaoNullSpace.cpp:915
MoFEM::createTao
auto createTao(MPI_Comm comm)
Definition PetscSmartObj.hpp:251
g
constexpr double g
Definition shallow_wave.cpp:63
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
Generated by Doxygen 1.9.8 and hosted at University of Glasgow