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:   /* Compute R3 and R5 */

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

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

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

 28:   VecNorm(tao->gradient, NORM_1, &qp->dinfeas);
 29:   qp->rnorm = (qp->dinfeas + qp->pinfeas) / (qp->m + qp->n);
 30:   return 0;
 31: }

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

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


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

 55:   VecSet(qp->GZwork, p01);
 56:   VecSet(qp->TSwork, p01);

 58:   VecPointwiseMax(qp->G, qp->G, qp->GZwork);
 59:   VecPointwiseMax(qp->T, qp->T, qp->TSwork);

 61:   /* Initialize Dual Variable Vectors */
 62:   VecCopy(qp->G, qp->Z);
 63:   VecReciprocal(qp->Z);

 65:   VecCopy(qp->T, qp->S);
 66:   VecReciprocal(qp->S);

 68:   MatMult(tao->hessian, qp->Work, qp->RHS);
 69:   VecAbs(qp->RHS);
 70:   VecSet(qp->Work, p01);
 71:   VecPointwiseMax(qp->RHS, qp->RHS, qp->Work);

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

 77:   VecScale(qp->S, mu);
 78:   VecScale(qp->Z, mu);

 80:   VecSet(qp->TSwork, p01);
 81:   VecSet(qp->GZwork, p01);
 82:   VecPointwiseMax(qp->S, qp->S, qp->TSwork);
 83:   VecPointwiseMax(qp->Z, qp->Z, qp->GZwork);

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

 94:     QPIPComputeResidual(qp, tao);

 96:     VecCopy(tao->solution, qp->R3);
 97:     VecAXPY(qp->R3, -1.0, qp->G);
 98:     VecAXPY(qp->R3, -1.0, qp->XL);

100:     VecCopy(tao->solution, qp->R5);
101:     VecAXPY(qp->R5, 1.0, qp->T);
102:     VecAXPY(qp->R5, -1.0, qp->XU);

104:     VecNorm(qp->R3, NORM_INFINITY, &gap1);
105:     VecNorm(qp->R5, NORM_INFINITY, &gap2);
106:     qp->pinfeas = PetscMax(gap1, gap2);

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

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

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

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

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

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

140:   qp->psteplength = PetscMin(qp->psteplength, qp->dsteplength);
141:   qp->dsteplength = qp->psteplength;
142:   return 0;
143: }

145: static PetscErrorCode QPIPComputeNormFromCentralPath(TAO_BQPIP *qp, PetscReal *norm)
146: {
147:   PetscReal gap[2], mu[2], nmu;

149:   VecPointwiseMult(qp->GZwork, qp->G, qp->Z);
150:   VecPointwiseMult(qp->TSwork, qp->T, qp->S);
151:   VecNorm(qp->TSwork, NORM_1, &mu[0]);
152:   VecNorm(qp->GZwork, NORM_1, &mu[1]);

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

156:   VecShift(qp->GZwork, nmu);
157:   VecShift(qp->TSwork, nmu);

159:   VecNorm(qp->GZwork, NORM_2, &gap[0]);
160:   VecNorm(qp->TSwork, NORM_2, &gap[1]);
161:   gap[0] *= gap[0];
162:   gap[1] *= gap[1];

164:   qp->pathnorm = PetscSqrtScalar(gap[0] + gap[1]);
165:   *norm        = qp->pathnorm;
166:   return 0;
167: }

169: static PetscErrorCode QPIPComputeStepDirection(TAO_BQPIP *qp, Tao tao)
170: {
171:   /* Calculate DG */
172:   VecCopy(tao->stepdirection, qp->DG);
173:   VecAXPY(qp->DG, 1.0, qp->R3);

175:   /* Calculate DT */
176:   VecCopy(tao->stepdirection, qp->DT);
177:   VecScale(qp->DT, -1.0);
178:   VecAXPY(qp->DT, -1.0, qp->R5);

180:   /* Calculate DZ */
181:   VecAXPY(qp->DZ, -1.0, qp->Z);
182:   VecPointwiseDivide(qp->GZwork, qp->DG, qp->G);
183:   VecPointwiseMult(qp->GZwork, qp->GZwork, qp->Z);
184:   VecAXPY(qp->DZ, -1.0, qp->GZwork);

186:   /* Calculate DS */
187:   VecAXPY(qp->DS, -1.0, qp->S);
188:   VecPointwiseDivide(qp->TSwork, qp->DT, qp->T);
189:   VecPointwiseMult(qp->TSwork, qp->TSwork, qp->S);
190:   VecAXPY(qp->DS, -1.0, qp->TSwork);
191:   return 0;
192: }

194: static PetscErrorCode TaoSetup_BQPIP(Tao tao)
195: {
196:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;

198:   /* Set pointers to Data */
199:   VecGetSize(tao->solution, &qp->n);

201:   /* Allocate some arrays */
202:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
203:   if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);

205:   VecDuplicate(tao->solution, &qp->Work);
206:   VecDuplicate(tao->solution, &qp->XU);
207:   VecDuplicate(tao->solution, &qp->XL);
208:   VecDuplicate(tao->solution, &qp->HDiag);
209:   VecDuplicate(tao->solution, &qp->DiagAxpy);
210:   VecDuplicate(tao->solution, &qp->RHS);
211:   VecDuplicate(tao->solution, &qp->RHS2);
212:   VecDuplicate(tao->solution, &qp->C);

214:   VecDuplicate(tao->solution, &qp->G);
215:   VecDuplicate(tao->solution, &qp->DG);
216:   VecDuplicate(tao->solution, &qp->S);
217:   VecDuplicate(tao->solution, &qp->Z);
218:   VecDuplicate(tao->solution, &qp->DZ);
219:   VecDuplicate(tao->solution, &qp->GZwork);
220:   VecDuplicate(tao->solution, &qp->R3);

222:   VecDuplicate(tao->solution, &qp->T);
223:   VecDuplicate(tao->solution, &qp->DT);
224:   VecDuplicate(tao->solution, &qp->DS);
225:   VecDuplicate(tao->solution, &qp->TSwork);
226:   VecDuplicate(tao->solution, &qp->R5);
227:   qp->m = 2 * qp->n;
228:   return 0;
229: }

231: static PetscErrorCode TaoSolve_BQPIP(Tao tao)
232: {
233:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
234:   PetscInt   its;
235:   PetscReal  d1, d2, ksptol, sigmamu;
236:   PetscReal  gnorm, dstep, pstep, step = 0;
237:   PetscReal  gap[4];
238:   PetscBool  getdiagop;

240:   qp->dobj        = 0.0;
241:   qp->pobj        = 1.0;
242:   qp->gap         = 10.0;
243:   qp->rgap        = 1.0;
244:   qp->mu          = 1.0;
245:   qp->dinfeas     = 1.0;
246:   qp->psteplength = 0.0;
247:   qp->dsteplength = 0.0;

249:   /* TODO
250:      - Remove fixed variables and treat them correctly
251:      - Use index sets for the infinite versus finite bounds
252:      - Update remaining code for fixed and free variables
253:      - Fix inexact solves for predictor and corrector
254:   */

256:   /* Tighten infinite bounds, things break when we don't do this
257:     -- see test_bqpip.c
258:   */
259:   TaoComputeVariableBounds(tao);
260:   VecSet(qp->XU, 1.0e20);
261:   VecSet(qp->XL, -1.0e20);
262:   if (tao->XL) VecPointwiseMax(qp->XL, qp->XL, tao->XL);
263:   if (tao->XU) VecPointwiseMin(qp->XU, qp->XU, tao->XU);
264:   VecMedian(qp->XL, tao->solution, qp->XU, tao->solution);

266:   /* Evaluate gradient and Hessian at zero to get the correct values
267:      without contaminating them with numerical artifacts.
268:   */
269:   VecSet(qp->Work, 0);
270:   TaoComputeObjectiveAndGradient(tao, qp->Work, &qp->d, qp->C);
271:   TaoComputeHessian(tao, qp->Work, tao->hessian, tao->hessian_pre);
272:   MatHasOperation(tao->hessian, MATOP_GET_DIAGONAL, &getdiagop);
273:   if (getdiagop) MatGetDiagonal(tao->hessian, qp->HDiag);

275:   /* Initialize starting point and residuals */
276:   QPIPSetInitialPoint(qp, tao);
277:   QPIPComputeResidual(qp, tao);

279:   /* Enter main loop */
280:   tao->reason = TAO_CONTINUE_ITERATING;
281:   while (1) {
282:     /* Check Stopping Condition      */
283:     gnorm = PetscSqrtScalar(qp->gap + qp->dinfeas);
284:     TaoLogConvergenceHistory(tao, qp->pobj, gnorm, qp->pinfeas, tao->ksp_its);
285:     TaoMonitor(tao, tao->niter, qp->pobj, gnorm, qp->pinfeas, step);
286:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
287:     if (tao->reason != TAO_CONTINUE_ITERATING) break;
288:     /* Call general purpose update function */
289:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
290:     tao->niter++;
291:     tao->ksp_its = 0;

293:     /*
294:        Dual Infeasibility Direction should already be in the right
295:        hand side from computing the residuals
296:     */

298:     QPIPComputeNormFromCentralPath(qp, &d1);

300:     VecSet(qp->DZ, 0.0);
301:     VecSet(qp->DS, 0.0);

303:     /*
304:        Compute the Primal Infeasiblitiy RHS and the
305:        Diagonal Matrix to be added to H and store in Work
306:     */
307:     VecPointwiseDivide(qp->DiagAxpy, qp->Z, qp->G);
308:     VecPointwiseMult(qp->GZwork, qp->DiagAxpy, qp->R3);
309:     VecAXPY(qp->RHS, -1.0, qp->GZwork);

311:     VecPointwiseDivide(qp->TSwork, qp->S, qp->T);
312:     VecAXPY(qp->DiagAxpy, 1.0, qp->TSwork);
313:     VecPointwiseMult(qp->TSwork, qp->TSwork, qp->R5);
314:     VecAXPY(qp->RHS, -1.0, qp->TSwork);

316:     /*  Determine the solving tolerance */
317:     ksptol = qp->mu / 10.0;
318:     ksptol = PetscMin(ksptol, 0.001);
319:     KSPSetTolerances(tao->ksp, ksptol, 1e-30, 1e30, PetscMax(10, qp->n));

321:     /* Shift the diagonals of the Hessian matrix */
322:     MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES);
323:     if (!getdiagop) {
324:       VecCopy(qp->DiagAxpy, qp->HDiag);
325:       VecScale(qp->HDiag, -1.0);
326:     }
327:     MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY);
328:     MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY);

330:     KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);
331:     KSPSolve(tao->ksp, qp->RHS, tao->stepdirection);
332:     KSPGetIterationNumber(tao->ksp, &its);
333:     tao->ksp_its += its;
334:     tao->ksp_tot_its += its;

336:     /* Restore the true diagonal of the Hessian matrix */
337:     if (getdiagop) {
338:       MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES);
339:     } else {
340:       MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES);
341:     }
342:     MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY);
343:     MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY);
344:     QPIPComputeStepDirection(qp, tao);
345:     QPIPStepLength(qp);

347:     /* Calculate New Residual R1 in Work vector */
348:     MatMult(tao->hessian, tao->stepdirection, qp->RHS2);
349:     VecAXPY(qp->RHS2, 1.0, qp->DS);
350:     VecAXPY(qp->RHS2, -1.0, qp->DZ);
351:     VecAYPX(qp->RHS2, qp->dsteplength, tao->gradient);

353:     VecNorm(qp->RHS2, NORM_2, &qp->dinfeas);
354:     VecDot(qp->DZ, qp->DG, gap);
355:     VecDot(qp->DS, qp->DT, gap + 1);

357:     qp->rnorm = (qp->dinfeas + qp->psteplength * qp->pinfeas) / (qp->m + qp->n);
358:     pstep     = qp->psteplength;
359:     step      = PetscMin(qp->psteplength, qp->dsteplength);
360:     sigmamu   = (pstep * pstep * (gap[0] + gap[1]) + (1 - pstep) * qp->gap) / qp->m;

362:     if (qp->predcorr && step < 0.9) {
363:       if (sigmamu < qp->mu) {
364:         sigmamu = sigmamu / qp->mu;
365:         sigmamu = sigmamu * sigmamu * sigmamu;
366:       } else {
367:         sigmamu = 1.0;
368:       }
369:       sigmamu = sigmamu * qp->mu;

371:       /* Compute Corrector Step */
372:       VecPointwiseMult(qp->DZ, qp->DG, qp->DZ);
373:       VecScale(qp->DZ, -1.0);
374:       VecShift(qp->DZ, sigmamu);
375:       VecPointwiseDivide(qp->DZ, qp->DZ, qp->G);

377:       VecPointwiseMult(qp->DS, qp->DS, qp->DT);
378:       VecScale(qp->DS, -1.0);
379:       VecShift(qp->DS, sigmamu);
380:       VecPointwiseDivide(qp->DS, qp->DS, qp->T);

382:       VecCopy(qp->DZ, qp->RHS2);
383:       VecAXPY(qp->RHS2, -1.0, qp->DS);
384:       VecAXPY(qp->RHS2, 1.0, qp->RHS);

386:       /* Approximately solve the linear system */
387:       MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES);
388:       if (!getdiagop) {
389:         VecCopy(qp->DiagAxpy, qp->HDiag);
390:         VecScale(qp->HDiag, -1.0);
391:       }
392:       MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY);
393:       MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY);

395:       /* Solve using the previous tolerances that were set */
396:       KSPSolve(tao->ksp, qp->RHS2, tao->stepdirection);
397:       KSPGetIterationNumber(tao->ksp, &its);
398:       tao->ksp_its += its;
399:       tao->ksp_tot_its += its;

401:       if (getdiagop) {
402:         MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES);
403:       } else {
404:         MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES);
405:       }
406:       MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY);
407:       MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY);
408:       QPIPComputeStepDirection(qp, tao);
409:       QPIPStepLength(qp);
410:     } /* End Corrector step */

412:     /* Take the step */
413:     dstep = qp->dsteplength;

415:     VecAXPY(qp->Z, dstep, qp->DZ);
416:     VecAXPY(qp->S, dstep, qp->DS);
417:     VecAXPY(tao->solution, dstep, tao->stepdirection);
418:     VecAXPY(qp->G, dstep, qp->DG);
419:     VecAXPY(qp->T, dstep, qp->DT);

421:     /* Compute Residuals */
422:     QPIPComputeResidual(qp, tao);

424:     /* Evaluate quadratic function */
425:     MatMult(tao->hessian, tao->solution, qp->Work);

427:     VecDot(tao->solution, qp->Work, &d1);
428:     VecDot(tao->solution, qp->C, &d2);
429:     VecDot(qp->G, qp->Z, gap);
430:     VecDot(qp->T, qp->S, gap + 1);

432:     /* Compute the duality gap */
433:     qp->pobj = d1 / 2.0 + d2 + qp->d;
434:     qp->gap  = gap[0] + gap[1];
435:     qp->dobj = qp->pobj - qp->gap;
436:     if (qp->m > 0) qp->mu = qp->gap / (qp->m);
437:     qp->rgap = qp->gap / (PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
438:   } /* END MAIN LOOP  */
439:   return 0;
440: }

442: static PetscErrorCode TaoView_BQPIP(Tao tao, PetscViewer viewer)
443: {
444:   return 0;
445: }

447: static PetscErrorCode TaoSetFromOptions_BQPIP(Tao tao, PetscOptionItems *PetscOptionsObject)
448: {
449:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;

451:   PetscOptionsHeadBegin(PetscOptionsObject, "Interior point method for bound constrained quadratic optimization");
452:   PetscOptionsInt("-tao_bqpip_predcorr", "Use a predictor-corrector method", "", qp->predcorr, &qp->predcorr, NULL);
453:   PetscOptionsHeadEnd();
454:   KSPSetFromOptions(tao->ksp);
455:   return 0;
456: }

458: static PetscErrorCode TaoDestroy_BQPIP(Tao tao)
459: {
460:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;

462:   if (tao->setupcalled) {
463:     VecDestroy(&qp->G);
464:     VecDestroy(&qp->DG);
465:     VecDestroy(&qp->Z);
466:     VecDestroy(&qp->DZ);
467:     VecDestroy(&qp->GZwork);
468:     VecDestroy(&qp->R3);
469:     VecDestroy(&qp->S);
470:     VecDestroy(&qp->DS);
471:     VecDestroy(&qp->T);

473:     VecDestroy(&qp->DT);
474:     VecDestroy(&qp->TSwork);
475:     VecDestroy(&qp->R5);
476:     VecDestroy(&qp->HDiag);
477:     VecDestroy(&qp->Work);
478:     VecDestroy(&qp->XL);
479:     VecDestroy(&qp->XU);
480:     VecDestroy(&qp->DiagAxpy);
481:     VecDestroy(&qp->RHS);
482:     VecDestroy(&qp->RHS2);
483:     VecDestroy(&qp->C);
484:   }
485:   KSPDestroy(&tao->ksp);
486:   PetscFree(tao->data);
487:   return 0;
488: }

490: static PetscErrorCode TaoComputeDual_BQPIP(Tao tao, Vec DXL, Vec DXU)
491: {
492:   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;

494:   VecCopy(qp->Z, DXL);
495:   VecCopy(qp->S, DXU);
496:   VecScale(DXU, -1.0);
497:   return 0;
498: }

500: /*MC
501:  TAOBQPIP - interior-point method for quadratic programs with
502:     box constraints.

504:  Notes:
505:     This algorithms solves quadratic problems only, the Hessian will
506:         only be computed once.

508:  Options Database Keys:
509: . -tao_bqpip_predcorr - use a predictor/corrector method

511:   Level: beginner
512: M*/

514: PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao)
515: {
516:   TAO_BQPIP *qp;

518:   PetscNew(&qp);

520:   tao->ops->setup          = TaoSetup_BQPIP;
521:   tao->ops->solve          = TaoSolve_BQPIP;
522:   tao->ops->view           = TaoView_BQPIP;
523:   tao->ops->setfromoptions = TaoSetFromOptions_BQPIP;
524:   tao->ops->destroy        = TaoDestroy_BQPIP;
525:   tao->ops->computedual    = TaoComputeDual_BQPIP;

527:   /* Override default settings (unless already changed) */
528:   if (!tao->max_it_changed) tao->max_it = 100;
529:   if (!tao->max_funcs_changed) tao->max_funcs = 500;
530: #if defined(PETSC_USE_REAL_SINGLE)
531:   if (!tao->catol_changed) tao->catol = 1e-6;
532: #else
533:   if (!tao->catol_changed) tao->catol = 1e-12;
534: #endif

536:   /* Initialize pointers and variables */
537:   qp->n = 0;
538:   qp->m = 0;

540:   qp->predcorr = 1;
541:   tao->data    = (void *)qp;

543:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
544:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
545:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
546:   KSPSetType(tao->ksp, KSPCG);
547:   KSPSetTolerances(tao->ksp, 1e-14, 1e-30, 1e30, PetscMax(10, qp->n));
548:   return 0;
549: }