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:   PetscCall(TaoParametersInitialize(tao));
    541:   PetscObjectParameterSetDefault(tao, max_it, 100);
    542:   PetscObjectParameterSetDefault(tao, max_funcs, 500);
    543:   PetscObjectParameterSetDefault(tao, catol, PetscDefined(USE_REAL_SINGLE) ? 1e-6 : 1e-12);

    545:   /* Initialize pointers and variables */
    546:   qp->n = 0;
    547:   qp->m = 0;

    549:   qp->predcorr = 1;
    550:   tao->data    = (void *)qp;

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