Actual source code: bqpip.c

  1: /*
  2:     This file implements a Mehrotra predictor-corrector method for
  3:     bound-constrained quadratic programs.

  5:  */

  7: #include <../src/tao/quadratic/impls/bqpip/bqpipimpl.h>
  8: #include <petscksp.h>

 10: static PetscErrorCode QPIPComputeResidual(TAO_BQPIP *qp, Tao tao)
 11: {
 12:   PetscReal dtmp = 1.0 - qp->psteplength;

 14:   PetscFunctionBegin;
 15:   /* Compute R3 and R5 */

 17:   PetscCall(VecScale(qp->R3, dtmp));
 18:   PetscCall(VecScale(qp->R5, dtmp));
 19:   qp->pinfeas = dtmp * qp->pinfeas;

 21:   PetscCall(VecCopy(qp->S, tao->gradient));
 22:   PetscCall(VecAXPY(tao->gradient, -1.0, qp->Z));

 24:   PetscCall(MatMult(tao->hessian, tao->solution, qp->RHS));
 25:   PetscCall(VecScale(qp->RHS, -1.0));
 26:   PetscCall(VecAXPY(qp->RHS, -1.0, qp->C));
 27:   PetscCall(VecAXPY(tao->gradient, -1.0, qp->RHS));

 29:   PetscCall(VecNorm(tao->gradient, NORM_1, &qp->dinfeas));
 30:   qp->rnorm = (qp->dinfeas + qp->pinfeas) / (qp->m + qp->n);
 31:   PetscFunctionReturn(PETSC_SUCCESS);
 32: }

 34: static PetscErrorCode QPIPSetInitialPoint(TAO_BQPIP *qp, Tao tao)
 35: {
 36:   PetscReal two = 2.0, p01 = 1;
 37:   PetscReal gap1, gap2, fff, mu;

 39:   PetscFunctionBegin;
 40:   /* Compute function, Gradient R=Hx+b, and Hessian */
 41:   PetscCall(MatMult(tao->hessian, tao->solution, tao->gradient));
 42:   PetscCall(VecCopy(qp->C, qp->Work));
 43:   PetscCall(VecAXPY(qp->Work, 0.5, tao->gradient));
 44:   PetscCall(VecAXPY(tao->gradient, 1.0, qp->C));
 45:   PetscCall(VecDot(tao->solution, qp->Work, &fff));
 46:   qp->pobj = fff + qp->d;

 48:   PetscCheck(!PetscIsInfOrNanReal(qp->pobj), PETSC_COMM_SELF, PETSC_ERR_USER, "User provided data contains Inf or NaN");

 50:   /* Initialize slack vectors */
 51:   /* T = XU - X; G = X - XL */
 52:   PetscCall(VecCopy(qp->XU, qp->T));
 53:   PetscCall(VecAXPY(qp->T, -1.0, tao->solution));
 54:   PetscCall(VecCopy(tao->solution, qp->G));
 55:   PetscCall(VecAXPY(qp->G, -1.0, qp->XL));

 57:   PetscCall(VecSet(qp->GZwork, p01));
 58:   PetscCall(VecSet(qp->TSwork, p01));

 60:   PetscCall(VecPointwiseMax(qp->G, qp->G, qp->GZwork));
 61:   PetscCall(VecPointwiseMax(qp->T, qp->T, qp->TSwork));

 63:   /* Initialize Dual Variable Vectors */
 64:   PetscCall(VecCopy(qp->G, qp->Z));
 65:   PetscCall(VecReciprocal(qp->Z));

 67:   PetscCall(VecCopy(qp->T, qp->S));
 68:   PetscCall(VecReciprocal(qp->S));

 70:   PetscCall(MatMult(tao->hessian, qp->Work, qp->RHS));
 71:   PetscCall(VecAbs(qp->RHS));
 72:   PetscCall(VecSet(qp->Work, p01));
 73:   PetscCall(VecPointwiseMax(qp->RHS, qp->RHS, qp->Work));

 75:   PetscCall(VecPointwiseDivide(qp->RHS, tao->gradient, qp->RHS));
 76:   PetscCall(VecNorm(qp->RHS, NORM_1, &gap1));
 77:   mu = PetscMin(10.0, (gap1 + 10.0) / qp->m);

 79:   PetscCall(VecScale(qp->S, mu));
 80:   PetscCall(VecScale(qp->Z, mu));

 82:   PetscCall(VecSet(qp->TSwork, p01));
 83:   PetscCall(VecSet(qp->GZwork, p01));
 84:   PetscCall(VecPointwiseMax(qp->S, qp->S, qp->TSwork));
 85:   PetscCall(VecPointwiseMax(qp->Z, qp->Z, qp->GZwork));

 87:   qp->mu      = 0;
 88:   qp->dinfeas = 1.0;
 89:   qp->pinfeas = 1.0;
 90:   while ((qp->dinfeas + qp->pinfeas) / (qp->m + qp->n) >= qp->mu) {
 91:     PetscCall(VecScale(qp->G, two));
 92:     PetscCall(VecScale(qp->Z, two));
 93:     PetscCall(VecScale(qp->S, two));
 94:     PetscCall(VecScale(qp->T, two));

 96:     PetscCall(QPIPComputeResidual(qp, tao));

 98:     PetscCall(VecCopy(tao->solution, qp->R3));
 99:     PetscCall(VecAXPY(qp->R3, -1.0, qp->G));
100:     PetscCall(VecAXPY(qp->R3, -1.0, qp->XL));

102:     PetscCall(VecCopy(tao->solution, qp->R5));
103:     PetscCall(VecAXPY(qp->R5, 1.0, qp->T));
104:     PetscCall(VecAXPY(qp->R5, -1.0, qp->XU));

106:     PetscCall(VecNorm(qp->R3, NORM_INFINITY, &gap1));
107:     PetscCall(VecNorm(qp->R5, NORM_INFINITY, &gap2));
108:     qp->pinfeas = PetscMax(gap1, gap2);

110:     /* Compute the duality gap */
111:     PetscCall(VecDot(qp->G, qp->Z, &gap1));
112:     PetscCall(VecDot(qp->T, qp->S, &gap2));

114:     qp->gap  = gap1 + gap2;
115:     qp->dobj = qp->pobj - qp->gap;
116:     if (qp->m > 0) {
117:       qp->mu = qp->gap / (qp->m);
118:     } else {
119:       qp->mu = 0.0;
120:     }
121:     qp->rgap = qp->gap / (PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
122:   }
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: static PetscErrorCode QPIPStepLength(TAO_BQPIP *qp)
127: {
128:   PetscReal tstep1, tstep2, tstep3, tstep4, tstep;

130:   PetscFunctionBegin;
131:   /* Compute stepsize to the boundary */
132:   PetscCall(VecStepMax(qp->G, qp->DG, &tstep1));
133:   PetscCall(VecStepMax(qp->T, qp->DT, &tstep2));
134:   PetscCall(VecStepMax(qp->S, qp->DS, &tstep3));
135:   PetscCall(VecStepMax(qp->Z, qp->DZ, &tstep4));

137:   tstep           = PetscMin(tstep1, tstep2);
138:   qp->psteplength = PetscMin(0.95 * tstep, 1.0);

140:   tstep           = PetscMin(tstep3, tstep4);
141:   qp->dsteplength = PetscMin(0.95 * tstep, 1.0);

143:   qp->psteplength = PetscMin(qp->psteplength, qp->dsteplength);
144:   qp->dsteplength = qp->psteplength;
145:   PetscFunctionReturn(PETSC_SUCCESS);
146: }

148: static PetscErrorCode QPIPComputeNormFromCentralPath(TAO_BQPIP *qp, PetscReal *norm)
149: {
150:   PetscReal gap[2], mu[2], nmu;

152:   PetscFunctionBegin;
153:   PetscCall(VecPointwiseMult(qp->GZwork, qp->G, qp->Z));
154:   PetscCall(VecPointwiseMult(qp->TSwork, qp->T, qp->S));
155:   PetscCall(VecNorm(qp->TSwork, NORM_1, &mu[0]));
156:   PetscCall(VecNorm(qp->GZwork, NORM_1, &mu[1]));

158:   nmu = -(mu[0] + mu[1]) / qp->m;

160:   PetscCall(VecShift(qp->GZwork, nmu));
161:   PetscCall(VecShift(qp->TSwork, nmu));

163:   PetscCall(VecNorm(qp->GZwork, NORM_2, &gap[0]));
164:   PetscCall(VecNorm(qp->TSwork, NORM_2, &gap[1]));
165:   gap[0] *= gap[0];
166:   gap[1] *= gap[1];

168:   qp->pathnorm = PetscSqrtScalar(gap[0] + gap[1]);
169:   *norm        = qp->pathnorm;
170:   PetscFunctionReturn(PETSC_SUCCESS);
171: }

173: static PetscErrorCode QPIPComputeStepDirection(TAO_BQPIP *qp, Tao tao)
174: {
175:   PetscFunctionBegin;
176:   /* Calculate DG */
177:   PetscCall(VecCopy(tao->stepdirection, qp->DG));
178:   PetscCall(VecAXPY(qp->DG, 1.0, qp->R3));

180:   /* Calculate DT */
181:   PetscCall(VecCopy(tao->stepdirection, qp->DT));
182:   PetscCall(VecScale(qp->DT, -1.0));
183:   PetscCall(VecAXPY(qp->DT, -1.0, qp->R5));

185:   /* Calculate DZ */
186:   PetscCall(VecAXPY(qp->DZ, -1.0, qp->Z));
187:   PetscCall(VecPointwiseDivide(qp->GZwork, qp->DG, qp->G));
188:   PetscCall(VecPointwiseMult(qp->GZwork, qp->GZwork, qp->Z));
189:   PetscCall(VecAXPY(qp->DZ, -1.0, qp->GZwork));

191:   /* Calculate DS */
192:   PetscCall(VecAXPY(qp->DS, -1.0, qp->S));
193:   PetscCall(VecPointwiseDivide(qp->TSwork, qp->DT, qp->T));
194:   PetscCall(VecPointwiseMult(qp->TSwork, qp->TSwork, qp->S));
195:   PetscCall(VecAXPY(qp->DS, -1.0, qp->TSwork));
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: static PetscErrorCode TaoSetup_BQPIP(Tao tao)
200: {
201:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;

203:   PetscFunctionBegin;
204:   /* Set pointers to Data */
205:   PetscCall(VecGetSize(tao->solution, &qp->n));

207:   /* Allocate some arrays */
208:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
209:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));

211:   PetscCall(VecDuplicate(tao->solution, &qp->Work));
212:   PetscCall(VecDuplicate(tao->solution, &qp->XU));
213:   PetscCall(VecDuplicate(tao->solution, &qp->XL));
214:   PetscCall(VecDuplicate(tao->solution, &qp->HDiag));
215:   PetscCall(VecDuplicate(tao->solution, &qp->DiagAxpy));
216:   PetscCall(VecDuplicate(tao->solution, &qp->RHS));
217:   PetscCall(VecDuplicate(tao->solution, &qp->RHS2));
218:   PetscCall(VecDuplicate(tao->solution, &qp->C));

220:   PetscCall(VecDuplicate(tao->solution, &qp->G));
221:   PetscCall(VecDuplicate(tao->solution, &qp->DG));
222:   PetscCall(VecDuplicate(tao->solution, &qp->S));
223:   PetscCall(VecDuplicate(tao->solution, &qp->Z));
224:   PetscCall(VecDuplicate(tao->solution, &qp->DZ));
225:   PetscCall(VecDuplicate(tao->solution, &qp->GZwork));
226:   PetscCall(VecDuplicate(tao->solution, &qp->R3));

228:   PetscCall(VecDuplicate(tao->solution, &qp->T));
229:   PetscCall(VecDuplicate(tao->solution, &qp->DT));
230:   PetscCall(VecDuplicate(tao->solution, &qp->DS));
231:   PetscCall(VecDuplicate(tao->solution, &qp->TSwork));
232:   PetscCall(VecDuplicate(tao->solution, &qp->R5));
233:   qp->m = 2 * qp->n;
234:   PetscFunctionReturn(PETSC_SUCCESS);
235: }

237: static PetscErrorCode TaoSolve_BQPIP(Tao tao)
238: {
239:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
240:   PetscInt   its;
241:   PetscReal  d1, d2, ksptol, sigmamu;
242:   PetscReal  gnorm, dstep, pstep, step = 0;
243:   PetscReal  gap[4];
244:   PetscBool  getdiagop;

246:   PetscFunctionBegin;
247:   qp->dobj        = 0.0;
248:   qp->pobj        = 1.0;
249:   qp->gap         = 10.0;
250:   qp->rgap        = 1.0;
251:   qp->mu          = 1.0;
252:   qp->dinfeas     = 1.0;
253:   qp->psteplength = 0.0;
254:   qp->dsteplength = 0.0;

256:   /* TODO
257:      - Remove fixed variables and treat them correctly
258:      - Use index sets for the infinite versus finite bounds
259:      - Update remaining code for fixed and free variables
260:      - Fix inexact solves for predictor and corrector
261:   */

263:   /* Tighten infinite bounds, things break when we don't do this
264:     -- see test_bqpip.c
265:   */
266:   PetscCall(TaoComputeVariableBounds(tao));
267:   PetscCall(VecSet(qp->XU, 1.0e20));
268:   PetscCall(VecSet(qp->XL, -1.0e20));
269:   if (tao->XL) PetscCall(VecPointwiseMax(qp->XL, qp->XL, tao->XL));
270:   if (tao->XU) PetscCall(VecPointwiseMin(qp->XU, qp->XU, tao->XU));
271:   PetscCall(VecMedian(qp->XL, tao->solution, qp->XU, tao->solution));

273:   /* Evaluate gradient and Hessian at zero to get the correct values
274:      without contaminating them with numerical artifacts.
275:   */
276:   PetscCall(VecSet(qp->Work, 0));
277:   PetscCall(TaoComputeObjectiveAndGradient(tao, qp->Work, &qp->d, qp->C));
278:   PetscCall(TaoComputeHessian(tao, qp->Work, tao->hessian, tao->hessian_pre));
279:   PetscCall(MatHasOperation(tao->hessian, MATOP_GET_DIAGONAL, &getdiagop));
280:   if (getdiagop) PetscCall(MatGetDiagonal(tao->hessian, qp->HDiag));

282:   /* Initialize starting point and residuals */
283:   PetscCall(QPIPSetInitialPoint(qp, tao));
284:   PetscCall(QPIPComputeResidual(qp, tao));

286:   /* Enter main loop */
287:   tao->reason = TAO_CONTINUE_ITERATING;
288:   while (1) {
289:     /* Check Stopping Condition      */
290:     gnorm = PetscSqrtScalar(qp->gap + qp->dinfeas);
291:     PetscCall(TaoLogConvergenceHistory(tao, qp->pobj, gnorm, qp->pinfeas, tao->ksp_its));
292:     PetscCall(TaoMonitor(tao, tao->niter, qp->pobj, gnorm, qp->pinfeas, step));
293:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
294:     if (tao->reason != TAO_CONTINUE_ITERATING) break;
295:     /* Call general purpose update function */
296:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
297:     tao->niter++;
298:     tao->ksp_its = 0;

300:     /*
301:        Dual Infeasibility Direction should already be in the right
302:        hand side from computing the residuals
303:     */

305:     PetscCall(QPIPComputeNormFromCentralPath(qp, &d1));

307:     PetscCall(VecSet(qp->DZ, 0.0));
308:     PetscCall(VecSet(qp->DS, 0.0));

310:     /*
311:        Compute the Primal Infeasiblitiy RHS and the
312:        Diagonal Matrix to be added to H and store in Work
313:     */
314:     PetscCall(VecPointwiseDivide(qp->DiagAxpy, qp->Z, qp->G));
315:     PetscCall(VecPointwiseMult(qp->GZwork, qp->DiagAxpy, qp->R3));
316:     PetscCall(VecAXPY(qp->RHS, -1.0, qp->GZwork));

318:     PetscCall(VecPointwiseDivide(qp->TSwork, qp->S, qp->T));
319:     PetscCall(VecAXPY(qp->DiagAxpy, 1.0, qp->TSwork));
320:     PetscCall(VecPointwiseMult(qp->TSwork, qp->TSwork, qp->R5));
321:     PetscCall(VecAXPY(qp->RHS, -1.0, qp->TSwork));

323:     /*  Determine the solving tolerance */
324:     ksptol = qp->mu / 10.0;
325:     ksptol = PetscMin(ksptol, 0.001);
326:     PetscCall(KSPSetTolerances(tao->ksp, ksptol, 1e-30, 1e30, PetscMax(10, qp->n)));

328:     /* Shift the diagonals of the Hessian matrix */
329:     PetscCall(MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES));
330:     if (!getdiagop) {
331:       PetscCall(VecCopy(qp->DiagAxpy, qp->HDiag));
332:       PetscCall(VecScale(qp->HDiag, -1.0));
333:     }
334:     PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
335:     PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));

337:     PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));
338:     PetscCall(KSPSolve(tao->ksp, qp->RHS, tao->stepdirection));
339:     PetscCall(KSPGetIterationNumber(tao->ksp, &its));
340:     tao->ksp_its += its;
341:     tao->ksp_tot_its += its;

343:     /* Restore the true diagonal of the Hessian matrix */
344:     if (getdiagop) {
345:       PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES));
346:     } else {
347:       PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES));
348:     }
349:     PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
350:     PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));
351:     PetscCall(QPIPComputeStepDirection(qp, tao));
352:     PetscCall(QPIPStepLength(qp));

354:     /* Calculate New Residual R1 in Work vector */
355:     PetscCall(MatMult(tao->hessian, tao->stepdirection, qp->RHS2));
356:     PetscCall(VecAXPY(qp->RHS2, 1.0, qp->DS));
357:     PetscCall(VecAXPY(qp->RHS2, -1.0, qp->DZ));
358:     PetscCall(VecAYPX(qp->RHS2, qp->dsteplength, tao->gradient));

360:     PetscCall(VecNorm(qp->RHS2, NORM_2, &qp->dinfeas));
361:     PetscCall(VecDot(qp->DZ, qp->DG, gap));
362:     PetscCall(VecDot(qp->DS, qp->DT, gap + 1));

364:     qp->rnorm = (qp->dinfeas + qp->psteplength * qp->pinfeas) / (qp->m + qp->n);
365:     pstep     = qp->psteplength;
366:     step      = PetscMin(qp->psteplength, qp->dsteplength);
367:     sigmamu   = (pstep * pstep * (gap[0] + gap[1]) + (1 - pstep) * qp->gap) / qp->m;

369:     if (qp->predcorr && step < 0.9) {
370:       if (sigmamu < qp->mu) {
371:         sigmamu = sigmamu / qp->mu;
372:         sigmamu = sigmamu * sigmamu * sigmamu;
373:       } else {
374:         sigmamu = 1.0;
375:       }
376:       sigmamu = sigmamu * qp->mu;

378:       /* Compute Corrector Step */
379:       PetscCall(VecPointwiseMult(qp->DZ, qp->DG, qp->DZ));
380:       PetscCall(VecScale(qp->DZ, -1.0));
381:       PetscCall(VecShift(qp->DZ, sigmamu));
382:       PetscCall(VecPointwiseDivide(qp->DZ, qp->DZ, qp->G));

384:       PetscCall(VecPointwiseMult(qp->DS, qp->DS, qp->DT));
385:       PetscCall(VecScale(qp->DS, -1.0));
386:       PetscCall(VecShift(qp->DS, sigmamu));
387:       PetscCall(VecPointwiseDivide(qp->DS, qp->DS, qp->T));

389:       PetscCall(VecCopy(qp->DZ, qp->RHS2));
390:       PetscCall(VecAXPY(qp->RHS2, -1.0, qp->DS));
391:       PetscCall(VecAXPY(qp->RHS2, 1.0, qp->RHS));

393:       /* Approximately solve the linear system */
394:       PetscCall(MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES));
395:       if (!getdiagop) {
396:         PetscCall(VecCopy(qp->DiagAxpy, qp->HDiag));
397:         PetscCall(VecScale(qp->HDiag, -1.0));
398:       }
399:       PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
400:       PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));

402:       /* Solve using the previous tolerances that were set */
403:       PetscCall(KSPSolve(tao->ksp, qp->RHS2, tao->stepdirection));
404:       PetscCall(KSPGetIterationNumber(tao->ksp, &its));
405:       tao->ksp_its += its;
406:       tao->ksp_tot_its += its;

408:       if (getdiagop) {
409:         PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES));
410:       } else {
411:         PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES));
412:       }
413:       PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
414:       PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));
415:       PetscCall(QPIPComputeStepDirection(qp, tao));
416:       PetscCall(QPIPStepLength(qp));
417:     } /* End Corrector step */

419:     /* Take the step */
420:     dstep = qp->dsteplength;

422:     PetscCall(VecAXPY(qp->Z, dstep, qp->DZ));
423:     PetscCall(VecAXPY(qp->S, dstep, qp->DS));
424:     PetscCall(VecAXPY(tao->solution, dstep, tao->stepdirection));
425:     PetscCall(VecAXPY(qp->G, dstep, qp->DG));
426:     PetscCall(VecAXPY(qp->T, dstep, qp->DT));

428:     /* Compute Residuals */
429:     PetscCall(QPIPComputeResidual(qp, tao));

431:     /* Evaluate quadratic function */
432:     PetscCall(MatMult(tao->hessian, tao->solution, qp->Work));

434:     PetscCall(VecDot(tao->solution, qp->Work, &d1));
435:     PetscCall(VecDot(tao->solution, qp->C, &d2));
436:     PetscCall(VecDot(qp->G, qp->Z, gap));
437:     PetscCall(VecDot(qp->T, qp->S, gap + 1));

439:     /* Compute the duality gap */
440:     qp->pobj = d1 / 2.0 + d2 + qp->d;
441:     qp->gap  = gap[0] + gap[1];
442:     qp->dobj = qp->pobj - qp->gap;
443:     if (qp->m > 0) qp->mu = qp->gap / (qp->m);
444:     qp->rgap = qp->gap / (PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
445:   } /* END MAIN LOOP  */
446:   PetscFunctionReturn(PETSC_SUCCESS);
447: }

449: static PetscErrorCode TaoView_BQPIP(Tao tao, PetscViewer viewer)
450: {
451:   PetscFunctionBegin;
452:   PetscFunctionReturn(PETSC_SUCCESS);
453: }

455: static PetscErrorCode TaoSetFromOptions_BQPIP(Tao tao, PetscOptionItems *PetscOptionsObject)
456: {
457:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;

459:   PetscFunctionBegin;
460:   PetscOptionsHeadBegin(PetscOptionsObject, "Interior point method for bound constrained quadratic optimization");
461:   PetscCall(PetscOptionsInt("-tao_bqpip_predcorr", "Use a predictor-corrector method", "", qp->predcorr, &qp->predcorr, NULL));
462:   PetscOptionsHeadEnd();
463:   PetscCall(KSPSetFromOptions(tao->ksp));
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }

467: static PetscErrorCode TaoDestroy_BQPIP(Tao tao)
468: {
469:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;

471:   PetscFunctionBegin;
472:   if (tao->setupcalled) {
473:     PetscCall(VecDestroy(&qp->G));
474:     PetscCall(VecDestroy(&qp->DG));
475:     PetscCall(VecDestroy(&qp->Z));
476:     PetscCall(VecDestroy(&qp->DZ));
477:     PetscCall(VecDestroy(&qp->GZwork));
478:     PetscCall(VecDestroy(&qp->R3));
479:     PetscCall(VecDestroy(&qp->S));
480:     PetscCall(VecDestroy(&qp->DS));
481:     PetscCall(VecDestroy(&qp->T));

483:     PetscCall(VecDestroy(&qp->DT));
484:     PetscCall(VecDestroy(&qp->TSwork));
485:     PetscCall(VecDestroy(&qp->R5));
486:     PetscCall(VecDestroy(&qp->HDiag));
487:     PetscCall(VecDestroy(&qp->Work));
488:     PetscCall(VecDestroy(&qp->XL));
489:     PetscCall(VecDestroy(&qp->XU));
490:     PetscCall(VecDestroy(&qp->DiagAxpy));
491:     PetscCall(VecDestroy(&qp->RHS));
492:     PetscCall(VecDestroy(&qp->RHS2));
493:     PetscCall(VecDestroy(&qp->C));
494:   }
495:   PetscCall(KSPDestroy(&tao->ksp));
496:   PetscCall(PetscFree(tao->data));
497:   PetscFunctionReturn(PETSC_SUCCESS);
498: }

500: static PetscErrorCode TaoComputeDual_BQPIP(Tao tao, Vec DXL, Vec DXU)
501: {
502:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;

504:   PetscFunctionBegin;
505:   PetscCall(VecCopy(qp->Z, DXL));
506:   PetscCall(VecCopy(qp->S, DXU));
507:   PetscCall(VecScale(DXU, -1.0));
508:   PetscFunctionReturn(PETSC_SUCCESS);
509: }

511: /*MC
512:  TAOBQPIP - interior-point method for quadratic programs with
513:     box constraints.

515:  Notes:
516:     This algorithms solves quadratic problems only, the Hessian will
517:         only be computed once.

519:  Options Database Keys:
520: . -tao_bqpip_predcorr - use a predictor/corrector method

522:   Level: beginner
523: M*/

525: PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao)
526: {
527:   TAO_BQPIP *qp;

529:   PetscFunctionBegin;
530:   PetscCall(PetscNew(&qp));

532:   tao->ops->setup          = TaoSetup_BQPIP;
533:   tao->ops->solve          = TaoSolve_BQPIP;
534:   tao->ops->view           = TaoView_BQPIP;
535:   tao->ops->setfromoptions = TaoSetFromOptions_BQPIP;
536:   tao->ops->destroy        = TaoDestroy_BQPIP;
537:   tao->ops->computedual    = TaoComputeDual_BQPIP;

539:   /* Override default settings (unless already changed) */
540:   if (!tao->max_it_changed) tao->max_it = 100;
541:   if (!tao->max_funcs_changed) tao->max_funcs = 500;
542: #if defined(PETSC_USE_REAL_SINGLE)
543:   if (!tao->catol_changed) tao->catol = 1e-6;
544: #else
545:   if (!tao->catol_changed) tao->catol = 1e-12;
546: #endif

548:   /* Initialize pointers and variables */
549:   qp->n = 0;
550:   qp->m = 0;

552:   qp->predcorr = 1;
553:   tao->data    = (void *)qp;

555:   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
556:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
557:   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
558:   PetscCall(KSPSetType(tao->ksp, KSPCG));
559:   PetscCall(KSPSetTolerances(tao->ksp, 1e-14, 1e-30, 1e30, PetscMax(10, qp->n)));
560:   PetscFunctionReturn(PETSC_SUCCESS);
561: }