Actual source code: almmutils.c

  1: #include <../src/tao/constrained/impls/almm/almm.h>
  2: #include <petsctao.h>
  3: #include <petsc/private/petscimpl.h>
  4: #include <petsc/private/vecimpl.h>

  6: /*@
  7:    TaoALMMGetType - Retrieve the augmented Lagrangian formulation type for the subproblem.

  9:    Input Parameter:
 10: .  tao - the `Tao` context for the `TAOALMM` solver

 12:    Output Parameter:
 13: .  type - augmented Lagragrangian type

 15:    Level: advanced

 17: .seealso: `Tao`, `TAOALMM`, `TaoALMMSetType()`, `TaoALMMType`
 18: @*/
 19: PetscErrorCode TaoALMMGetType(Tao tao, TaoALMMType *type)
 20: {
 21:   PetscFunctionBegin;
 24:   PetscUseMethod(tao, "TaoALMMGetType_C", (Tao, TaoALMMType *), (tao, type));
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

 28: PetscErrorCode TaoALMMGetType_Private(Tao tao, TaoALMMType *type)
 29: {
 30:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

 32:   PetscFunctionBegin;
 33:   *type = auglag->type;
 34:   PetscFunctionReturn(PETSC_SUCCESS);
 35: }

 37: /*@
 38:    TaoALMMSetType - Determine the augmented Lagrangian formulation type for the subproblem.

 40:    Input Parameters:
 41: +  tao - the `Tao` context for the `TAOALMM` solver
 42: -  type - augmented Lagragrangian type

 44:    Level: advanced

 46: .seealso: `Tao`, `TAOALMM`, `TaoALMMGetType()`, `TaoALMMType`
 47: @*/
 48: PetscErrorCode TaoALMMSetType(Tao tao, TaoALMMType type)
 49: {
 50:   PetscFunctionBegin;
 52:   PetscTryMethod(tao, "TaoALMMSetType_C", (Tao, TaoALMMType), (tao, type));
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: PetscErrorCode TaoALMMSetType_Private(Tao tao, TaoALMMType type)
 57: {
 58:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

 60:   PetscFunctionBegin;
 61:   PetscCheck(!tao->setupcalled, PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoALMMSetType() must be called before TaoSetUp()");
 62:   auglag->type = type;
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: /*@
 67:    TaoALMMGetSubsolver - Retrieve the subsolver being used by `TAOALMM`.

 69:    Input Parameter:
 70: .  tao - the `Tao` context for the `TAOALMM` solver

 72:    Output Parameter:
 73: .  subsolver - the `Tao` context for the subsolver

 75:    Level: advanced

 77: .seealso: `Tao`, `TAOALMM`, `TaoALMMSetSubsolver()`
 78: @*/
 79: PetscErrorCode TaoALMMGetSubsolver(Tao tao, Tao *subsolver)
 80: {
 81:   PetscFunctionBegin;
 84:   PetscUseMethod(tao, "TaoALMMGetSubsolver_C", (Tao, Tao *), (tao, subsolver));
 85:   PetscFunctionReturn(PETSC_SUCCESS);
 86: }

 88: PetscErrorCode TaoALMMGetSubsolver_Private(Tao tao, Tao *subsolver)
 89: {
 90:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

 92:   PetscFunctionBegin;
 93:   *subsolver = auglag->subsolver;
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: /*@
 98:    TaoALMMSetSubsolver - Changes the subsolver inside `TAOALMM` with the user provided one.

100:    Input Parameters:
101: +  tao - the `Tao` context for the `TAOALMM` solver
102: -  subsolver - the Tao context for the subsolver

104:    Level: advanced

106:    Note:
107:    This is not recommended, instead call `TaoALMMGetSubsolver()` and set the type as desired.

109: .seealso: `Tao`, `TAOALMM`, `TaoALMMGetSubsolver()`
110: @*/
111: PetscErrorCode TaoALMMSetSubsolver(Tao tao, Tao subsolver)
112: {
113:   PetscFunctionBegin;
116:   PetscTryMethod(tao, "TaoALMMSetSubsolver_C", (Tao, Tao), (tao, subsolver));
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: PetscErrorCode TaoALMMSetSubsolver_Private(Tao tao, Tao subsolver)
121: {
122:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
123:   PetscBool compatible;

125:   PetscFunctionBegin;
126:   if (subsolver == auglag->subsolver) PetscFunctionReturn(PETSC_SUCCESS);
127:   if (tao->bounded) {
128:     PetscCall(PetscObjectTypeCompareAny((PetscObject)subsolver, &compatible, TAOSHELL, TAOBNCG, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, ""));
129:     PetscCheck(compatible, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Subsolver must be a bound-constrained first-order method");
130:   } else {
131:     PetscCall(PetscObjectTypeCompareAny((PetscObject)subsolver, &compatible, TAOSHELL, TAOCG, TAOLMVM, TAOBNCG, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, ""));
132:     PetscCheck(compatible, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Subsolver must be a first-order method");
133:   }
134:   PetscCheck(compatible, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Subsolver must be a first-order method");
135:   PetscCall(PetscObjectReference((PetscObject)subsolver));
136:   PetscCall(TaoDestroy(&auglag->subsolver));
137:   auglag->subsolver = subsolver;
138:   if (tao->setupcalled) {
139:     PetscCall(TaoSetSolution(auglag->subsolver, auglag->P));
140:     PetscCall(TaoSetObjective(auglag->subsolver, TaoALMMSubsolverObjective_Private, (void *)auglag));
141:     PetscCall(TaoSetObjectiveAndGradient(auglag->subsolver, NULL, TaoALMMSubsolverObjectiveAndGradient_Private, (void *)auglag));
142:     PetscCall(TaoSetVariableBounds(auglag->subsolver, auglag->PL, auglag->PU));
143:   }
144:   PetscFunctionReturn(PETSC_SUCCESS);
145: }

147: /*@
148:    TaoALMMGetMultipliers - Retrieve a pointer to the Lagrange multipliers.

150:    Input Parameter:
151: .  tao - the `Tao` context for the `TAOALMM` solver

153:    Output Parameter:
154: .  Y - vector of Lagrange multipliers

156:    Level: advanced

158:    Notes:
159:    For problems with both equality and inequality constraints,
160:    the multipliers are combined together as Y = (Ye, Yi). Users
161:    can recover copies of the subcomponents using index sets
162:    provided by `TaoALMMGetDualIS()` and use `VecGetSubVector()`.

164: .seealso: `TAOALMM`, `Tao`, `TaoALMMSetMultipliers()`, `TaoALMMGetDualIS()`
165: @*/
166: PetscErrorCode TaoALMMGetMultipliers(Tao tao, Vec *Y)
167: {
168:   PetscFunctionBegin;
171:   PetscUseMethod(tao, "TaoALMMGetMultipliers_C", (Tao, Vec *), (tao, Y));
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: PetscErrorCode TaoALMMGetMultipliers_Private(Tao tao, Vec *Y)
176: {
177:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

179:   PetscFunctionBegin;
180:   PetscCheck(tao->setupcalled, PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetUp() must be called first for scatters to be constructed");
181:   *Y = auglag->Y;
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: /*@
186:    TaoALMMSetMultipliers - Set user-defined Lagrange multipliers.

188:    Input Parameters:
189: +  tao - the `Tao` context for the `TAOALMM` solver
190: -  Y - vector of Lagrange multipliers

192:    Level: advanced

194:    Notes:
195:    The vector type and parallel layout must match the equality and inequality constraints.

197:    The vector must have a local size equal to the sum of the local sizes for the constraint vectors, and a
198:    global size equal to the sum of the global sizes of the constraint vectors.

200:    This routine is only useful if the user wants to change the
201:    parallel distribution of the combined dual vector in problems that
202:    feature both equality and inequality constraints. For other tasks,
203:    it is strongly recommended that the user retrieve the dual vector
204:    created by the solver using `TaoALMMGetMultipliers()`.

206: .seealso: `TAOALMM`, `Tao`, `TaoALMMGetMultipliers()`
207: @*/
208: PetscErrorCode TaoALMMSetMultipliers(Tao tao, Vec Y)
209: {
210:   PetscFunctionBegin;
213:   PetscTryMethod(tao, "TaoALMMSetMultipliers_C", (Tao, Vec), (tao, Y));
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: PetscErrorCode TaoALMMSetMultipliers_Private(Tao tao, Vec Y)
218: {
219:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
220:   VecType   Ytype;
221:   PetscInt  Nuser, Neq, Nineq, N;
222:   PetscBool same = PETSC_FALSE;

224:   PetscFunctionBegin;
225:   /* no-op if user provides vector from TaoALMMGetMultipliers() */
226:   if (Y == auglag->Y) PetscFunctionReturn(PETSC_SUCCESS);
227:   /* make sure vector type is same as equality and inequality constraints */
228:   if (tao->eq_constrained) {
229:     PetscCall(VecGetType(tao->constraints_equality, &Ytype));
230:   } else {
231:     PetscCall(VecGetType(tao->constraints_inequality, &Ytype));
232:   }
233:   PetscCall(PetscObjectTypeCompare((PetscObject)Y, Ytype, &same));
234:   PetscCheck(same, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given vector for multipliers is not the same type as constraint vectors");
235:   /* make sure global size matches sum of equality and inequality */
236:   if (tao->eq_constrained) {
237:     PetscCall(VecGetSize(tao->constraints_equality, &Neq));
238:   } else {
239:     Neq = 0;
240:   }
241:   if (tao->ineq_constrained) {
242:     PetscCall(VecGetSize(tao->constraints_inequality, &Nineq));
243:   } else {
244:     Nineq = 0;
245:   }
246:   N = Neq + Nineq;
247:   PetscCall(VecGetSize(Y, &Nuser));
248:   PetscCheck(Nuser == N, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given vector has wrong global size");
249:   /* if there is only one type of constraint, then we need the local size to match too */
250:   if (Neq == 0) {
251:     PetscCall(VecGetLocalSize(tao->constraints_inequality, &Nineq));
252:     PetscCall(VecGetLocalSize(Y, &Nuser));
253:     PetscCheck(Nuser == Nineq, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given vector has wrong local size");
254:   }
255:   if (Nineq == 0) {
256:     PetscCall(VecGetLocalSize(tao->constraints_equality, &Neq));
257:     PetscCall(VecGetLocalSize(Y, &Nuser));
258:     PetscCheck(Nuser == Neq, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given vector has wrong local size");
259:   }
260:   /* if we got here, the given vector is compatible so we can replace the current one */
261:   PetscCall(PetscObjectReference((PetscObject)Y));
262:   PetscCall(VecDestroy(&auglag->Y));
263:   auglag->Y = Y;
264:   /* if there are both types of constraints and the solver has already been set up,
265:      then we need to regenerate VecScatter objects for the new combined dual vector */
266:   if (tao->setupcalled && tao->eq_constrained && tao->ineq_constrained) {
267:     PetscCall(VecDestroy(&auglag->C));
268:     PetscCall(VecDuplicate(auglag->Y, &auglag->C));
269:     PetscCall(VecScatterDestroy(&auglag->Yscatter[0]));
270:     PetscCall(VecScatterCreate(auglag->Y, auglag->Yis[0], auglag->Ye, NULL, &auglag->Yscatter[0]));
271:     PetscCall(VecScatterDestroy(&auglag->Yscatter[1]));
272:     PetscCall(VecScatterCreate(auglag->Y, auglag->Yis[1], auglag->Yi, NULL, &auglag->Yscatter[1]));
273:   }
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }

277: /*@
278:    TaoALMMGetPrimalIS - Retrieve the index set that identifies optimization
279:                         and slack variable components of the subsolver's solution vector.

281:    Input Parameter:
282: .  tao - the `Tao` context for the `TAOALMM` solver

284:    Output Parameters:
285: +  opt_is - index set associated with the optimization variables (`NULL` if not needed)
286: -  slack_is - index set associated with the slack variables (`NULL` if not needed)

288:    Level: advanced

290: .seealso: `TAOALMM`, `Tao`, `IS`, `TaoALMMGetPrimalVector()`
291: @*/
292: PetscErrorCode TaoALMMGetPrimalIS(Tao tao, IS *opt_is, IS *slack_is)
293: {
294:   PetscFunctionBegin;
296:   PetscUseMethod(tao, "TaoALMMGetPrimalIS_C", (Tao, IS *, IS *), (tao, opt_is, slack_is));
297:   PetscFunctionReturn(PETSC_SUCCESS);
298: }

300: PetscErrorCode TaoALMMGetPrimalIS_Private(Tao tao, IS *opt_is, IS *slack_is)
301: {
302:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

304:   PetscFunctionBegin;
305:   PetscCheck(tao->ineq_constrained, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "Primal space has index sets only for inequality constrained problems");
306:   PetscCheck(tao->setupcalled, PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetUp() must be called first for index sets to be constructed");
307:   if (opt_is) *opt_is = auglag->Pis[0];
308:   if (slack_is) *slack_is = auglag->Pis[1];
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

312: /*@
313:    TaoALMMGetDualIS - Retrieve the index set that identifies equality
314:                       and inequality constraint components of the dual vector returned
315:                       by `TaoALMMGetMultipliers()`.

317:    Input Parameter:
318: .  tao - the Tao context for the `TAOALMM` solver

320:    Output Parameters:
321: +  eq_is - index set associated with the equality constraints (`NULL` if not needed)
322: -  ineq_is - index set associated with the inequality constraints (`NULL` if not needed)

324:    Level: advanced

326: .seealso: `TAOALMM`, `Tao`, `TaoALMMGetMultipliers()`
327: @*/
328: PetscErrorCode TaoALMMGetDualIS(Tao tao, IS *eq_is, IS *ineq_is)
329: {
330:   PetscFunctionBegin;
332:   PetscUseMethod(tao, "TaoALMMGetDualIS_C", (Tao, IS *, IS *), (tao, eq_is, ineq_is));
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }

336: PetscErrorCode TaoALMMGetDualIS_Private(Tao tao, IS *eq_is, IS *ineq_is)
337: {
338:   TAO_ALMM *auglag = (TAO_ALMM *)tao->data;

340:   PetscFunctionBegin;
342:   PetscCheck(tao->ineq_constrained && tao->ineq_constrained, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "Dual space has index sets only when problem has both equality and inequality constraints");
343:   PetscCheck(tao->setupcalled, PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetUp() must be called first for index sets to be constructed");
344:   if (eq_is) *eq_is = auglag->Yis[0];
345:   if (ineq_is) *ineq_is = auglag->Yis[1];
346:   PetscFunctionReturn(PETSC_SUCCESS);
347: }