Actual source code: bqpip.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  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: {
 13:   PetscReal      dtmp = 1.0 - qp->psteplength;

 16:   /* Compute R3 and R5 */

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

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

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

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

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

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

 50:   if (PetscIsInfOrNanReal(qp->pobj)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER, "User provided data contains Inf or NaN");

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

 59:   VecSet(qp->GZwork,p01);
 60:   VecSet(qp->TSwork,p01);

 62:   VecPointwiseMax(qp->G,qp->G,qp->GZwork);
 63:   VecPointwiseMax(qp->T,qp->T,qp->TSwork);

 65:   /* Initialize Dual Variable Vectors */
 66:   VecCopy(qp->G,qp->Z);
 67:   VecReciprocal(qp->Z);

 69:   VecCopy(qp->T,qp->S);
 70:   VecReciprocal(qp->S);

 72:   MatMult(tao->hessian,qp->Work,qp->RHS);
 73:   VecAbs(qp->RHS);
 74:   VecSet(qp->Work,p01);
 75:   VecPointwiseMax(qp->RHS,qp->RHS,qp->Work);

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

 81:   VecScale(qp->S,mu);
 82:   VecScale(qp->Z,mu);

 84:   VecSet(qp->TSwork,p01);
 85:   VecSet(qp->GZwork,p01);
 86:   VecPointwiseMax(qp->S,qp->S,qp->TSwork);
 87:   VecPointwiseMax(qp->Z,qp->Z,qp->GZwork);

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

 96:     QPIPComputeResidual(qp,tao);

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

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

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

110:     /* Compute the duality gap */
111:     VecDot(qp->G,qp->Z,&gap1);
112:     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:   return(0);
124: }

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

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

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

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

144:   qp->psteplength = PetscMin(qp->psteplength,qp->dsteplength);
145:   qp->dsteplength = qp->psteplength;
146:   return(0);
147: }

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

155:   VecPointwiseMult(qp->GZwork,qp->G,qp->Z);
156:   VecPointwiseMult(qp->TSwork,qp->T,qp->S);
157:   VecNorm(qp->TSwork,NORM_1,&mu[0]);
158:   VecNorm(qp->GZwork,NORM_1,&mu[1]);

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

162:   VecShift(qp->GZwork,nmu);
163:   VecShift(qp->TSwork,nmu);

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

170:   qp->pathnorm=PetscSqrtScalar(gap[0]+gap[1]);
171:   *norm=qp->pathnorm;
172:   return(0);
173: }

175: static PetscErrorCode QPIPComputeStepDirection(TAO_BQPIP *qp,Tao tao)
176: {

180:   /* Calculate DG */
181:   VecCopy(tao->stepdirection,qp->DG);
182:   VecAXPY(qp->DG,1.0,qp->R3);

184:   /* Calculate DT */
185:   VecCopy(tao->stepdirection,qp->DT);
186:   VecScale(qp->DT,-1.0);
187:   VecAXPY(qp->DT,-1.0,qp->R5);

189:   /* Calculate DZ */
190:   VecAXPY(qp->DZ,-1.0,qp->Z);
191:   VecPointwiseDivide(qp->GZwork,qp->DG,qp->G);
192:   VecPointwiseMult(qp->GZwork,qp->GZwork,qp->Z);
193:   VecAXPY(qp->DZ,-1.0,qp->GZwork);

195:   /* Calculate DS */
196:   VecAXPY(qp->DS,-1.0,qp->S);
197:   VecPointwiseDivide(qp->TSwork,qp->DT,qp->T);
198:   VecPointwiseMult(qp->TSwork,qp->TSwork,qp->S);
199:   VecAXPY(qp->DS,-1.0,qp->TSwork);
200:   return(0);
201: }

203: static PetscErrorCode TaoSetup_BQPIP(Tao tao)
204: {
205:   TAO_BQPIP      *qp =(TAO_BQPIP*)tao->data;

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

212:   /* Allocate some arrays */
213:   if (!tao->gradient) {
214:     VecDuplicate(tao->solution,&tao->gradient);
215:   }
216:   if (!tao->stepdirection) {
217:     VecDuplicate(tao->solution,&tao->stepdirection);
218:   }
219:   if (!tao->XL) {
220:     VecDuplicate(tao->solution,&tao->XL);
221:     VecSet(tao->XL,PETSC_NINFINITY);
222:   }
223:   if (!tao->XU) {
224:     VecDuplicate(tao->solution,&tao->XU);
225:     VecSet(tao->XU,PETSC_INFINITY);
226:   }

228:   VecDuplicate(tao->solution,&qp->Work);
229:   VecDuplicate(tao->solution,&qp->XU);
230:   VecDuplicate(tao->solution,&qp->XL);
231:   VecDuplicate(tao->solution,&qp->HDiag);
232:   VecDuplicate(tao->solution,&qp->DiagAxpy);
233:   VecDuplicate(tao->solution,&qp->RHS);
234:   VecDuplicate(tao->solution,&qp->RHS2);
235:   VecDuplicate(tao->solution,&qp->C);

237:   VecDuplicate(tao->solution,&qp->G);
238:   VecDuplicate(tao->solution,&qp->DG);
239:   VecDuplicate(tao->solution,&qp->S);
240:   VecDuplicate(tao->solution,&qp->Z);
241:   VecDuplicate(tao->solution,&qp->DZ);
242:   VecDuplicate(tao->solution,&qp->GZwork);
243:   VecDuplicate(tao->solution,&qp->R3);

245:   VecDuplicate(tao->solution,&qp->T);
246:   VecDuplicate(tao->solution,&qp->DT);
247:   VecDuplicate(tao->solution,&qp->DS);
248:   VecDuplicate(tao->solution,&qp->TSwork);
249:   VecDuplicate(tao->solution,&qp->R5);
250:   qp->m=2*qp->n;
251:   return(0);
252: }

254: static PetscErrorCode TaoSolve_BQPIP(Tao tao)
255: {
256:   TAO_BQPIP          *qp = (TAO_BQPIP*)tao->data;
257:   PetscErrorCode     ierr;
258:   PetscInt           its;
259:   PetscReal          d1,d2,ksptol,sigmamu;
260:   PetscReal          gnorm,dstep,pstep,step=0;
261:   PetscReal          gap[4];
262:   PetscBool          getdiagop;

265:   qp->dobj        = 0.0;
266:   qp->pobj        = 1.0;
267:   qp->gap         = 10.0;
268:   qp->rgap        = 1.0;
269:   qp->mu          = 1.0;
270:   qp->dinfeas     = 1.0;
271:   qp->psteplength = 0.0;
272:   qp->dsteplength = 0.0;

274:   /* TODO
275:      - Remove fixed variables and treat them correctly
276:      - Use index sets for the infinite versus finite bounds
277:      - Update remaining code for fixed and free variables
278:      - Fix inexact solves for predictor and corrector
279:   */

281:   /* Tighten infinite bounds, things break when we don't do this
282:     -- see test_bqpip.c
283:   */
284:   TaoComputeVariableBounds(tao);
285:   VecSet(qp->XU,1.0e20);
286:   VecSet(qp->XL,-1.0e20);
287:   VecPointwiseMax(qp->XL,qp->XL,tao->XL);
288:   VecPointwiseMin(qp->XU,qp->XU,tao->XU);
289:   VecMedian(qp->XL,tao->solution,qp->XU,tao->solution);

291:   /* Evaluate gradient and Hessian at zero to get the correct values
292:      without contaminating them with numerical artifacts.
293:   */
294:   VecSet(qp->Work,0);
295:   TaoComputeObjectiveAndGradient(tao,qp->Work,&qp->d,qp->C);
296:   TaoComputeHessian(tao,qp->Work,tao->hessian,tao->hessian_pre);
297:   MatHasOperation(tao->hessian,MATOP_GET_DIAGONAL,&getdiagop);
298:   if (getdiagop) {
299:     MatGetDiagonal(tao->hessian,qp->HDiag);
300:   }

302:   /* Initialize starting point and residuals */
303:   QPIPSetInitialPoint(qp,tao);
304:   QPIPComputeResidual(qp,tao);

306:   /* Enter main loop */
307:   tao->reason = TAO_CONTINUE_ITERATING;
308:   while (1) {

310:     /* Check Stopping Condition      */
311:     gnorm = PetscSqrtScalar(qp->gap + qp->dinfeas);
312:     TaoLogConvergenceHistory(tao,qp->pobj,gnorm,qp->pinfeas,tao->ksp_its);
313:     TaoMonitor(tao,tao->niter,qp->pobj,gnorm,qp->pinfeas,step);
314:     (*tao->ops->convergencetest)(tao,tao->cnvP);
315:     if (tao->reason != TAO_CONTINUE_ITERATING) break;
316:     /* Call general purpose update function */
317:     if (tao->ops->update) {
318:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
319:     }
320:     tao->niter++;
321:     tao->ksp_its = 0;

323:     /*
324:        Dual Infeasibility Direction should already be in the right
325:        hand side from computing the residuals
326:     */

328:     QPIPComputeNormFromCentralPath(qp,&d1);

330:     VecSet(qp->DZ,0.0);
331:     VecSet(qp->DS,0.0);

333:     /*
334:        Compute the Primal Infeasiblitiy RHS and the
335:        Diagonal Matrix to be added to H and store in Work
336:     */
337:     VecPointwiseDivide(qp->DiagAxpy,qp->Z,qp->G);
338:     VecPointwiseMult(qp->GZwork,qp->DiagAxpy,qp->R3);
339:     VecAXPY(qp->RHS,-1.0,qp->GZwork);

341:     VecPointwiseDivide(qp->TSwork,qp->S,qp->T);
342:     VecAXPY(qp->DiagAxpy,1.0,qp->TSwork);
343:     VecPointwiseMult(qp->TSwork,qp->TSwork,qp->R5);
344:     VecAXPY(qp->RHS,-1.0,qp->TSwork);

346:     /*  Determine the solving tolerance */
347:     ksptol = qp->mu/10.0;
348:     ksptol = PetscMin(ksptol,0.001);
349:     KSPSetTolerances(tao->ksp,ksptol,1e-30,1e30,PetscMax(10,qp->n));

351:     /* Shift the diagonals of the Hessian matrix */
352:     MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES);
353:     if (!getdiagop) {
354:       VecCopy(qp->DiagAxpy,qp->HDiag);
355:       VecScale(qp->HDiag,-1.0);
356:     }
357:     MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
358:     MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);

360:     KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre);
361:     KSPSolve(tao->ksp,qp->RHS,tao->stepdirection);
362:     KSPGetIterationNumber(tao->ksp,&its);
363:     tao->ksp_its += its;
364:     tao->ksp_tot_its += its;

366:     /* Restore the true diagonal of the Hessian matrix */
367:     if (getdiagop) {
368:       MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES);
369:     } else {
370:       MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES);
371:     }
372:     MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
373:     MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
374:     QPIPComputeStepDirection(qp,tao);
375:     QPIPStepLength(qp);

377:     /* Calculate New Residual R1 in Work vector */
378:     MatMult(tao->hessian,tao->stepdirection,qp->RHS2);
379:     VecAXPY(qp->RHS2,1.0,qp->DS);
380:     VecAXPY(qp->RHS2,-1.0,qp->DZ);
381:     VecAYPX(qp->RHS2,qp->dsteplength,tao->gradient);

383:     VecNorm(qp->RHS2,NORM_2,&qp->dinfeas);
384:     VecDot(qp->DZ,qp->DG,gap);
385:     VecDot(qp->DS,qp->DT,gap+1);

387:     qp->rnorm = (qp->dinfeas+qp->psteplength*qp->pinfeas)/(qp->m+qp->n);
388:     pstep     = qp->psteplength;
389:     step      = PetscMin(qp->psteplength,qp->dsteplength);
390:     sigmamu   = (pstep*pstep*(gap[0]+gap[1]) + (1 - pstep)*qp->gap)/qp->m;

392:     if (qp->predcorr && step < 0.9) {
393:       if (sigmamu < qp->mu) {
394:         sigmamu = sigmamu/qp->mu;
395:         sigmamu = sigmamu*sigmamu*sigmamu;
396:       } else {
397:         sigmamu = 1.0;
398:       }
399:       sigmamu = sigmamu*qp->mu;

401:       /* Compute Corrector Step */
402:       VecPointwiseMult(qp->DZ,qp->DG,qp->DZ);
403:       VecScale(qp->DZ,-1.0);
404:       VecShift(qp->DZ,sigmamu);
405:       VecPointwiseDivide(qp->DZ,qp->DZ,qp->G);

407:       VecPointwiseMult(qp->DS,qp->DS,qp->DT);
408:       VecScale(qp->DS,-1.0);
409:       VecShift(qp->DS,sigmamu);
410:       VecPointwiseDivide(qp->DS,qp->DS,qp->T);

412:       VecCopy(qp->DZ,qp->RHS2);
413:       VecAXPY(qp->RHS2,-1.0,qp->DS);
414:       VecAXPY(qp->RHS2,1.0,qp->RHS);

416:       /* Approximately solve the linear system */
417:       MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES);
418:       if (!getdiagop) {
419:         VecCopy(qp->DiagAxpy,qp->HDiag);
420:         VecScale(qp->HDiag,-1.0);
421:       }
422:       MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
423:       MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);

425:       /* Solve using the previous tolerances that were set */
426:       KSPSolve(tao->ksp,qp->RHS2,tao->stepdirection);
427:       KSPGetIterationNumber(tao->ksp,&its);
428:       tao->ksp_its += its;
429:       tao->ksp_tot_its += its;

431:       if (getdiagop) {
432:         MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES);
433:       } else {
434:         MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES);
435:       }
436:       MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
437:       MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
438:       QPIPComputeStepDirection(qp,tao);
439:       QPIPStepLength(qp);
440:     }  /* End Corrector step */


443:     /* Take the step */
444:     dstep = qp->dsteplength;

446:     VecAXPY(qp->Z,dstep,qp->DZ);
447:     VecAXPY(qp->S,dstep,qp->DS);
448:     VecAXPY(tao->solution,dstep,tao->stepdirection);
449:     VecAXPY(qp->G,dstep,qp->DG);
450:     VecAXPY(qp->T,dstep,qp->DT);

452:     /* Compute Residuals */
453:     QPIPComputeResidual(qp,tao);

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

458:     VecDot(tao->solution,qp->Work,&d1);
459:     VecDot(tao->solution,qp->C,&d2);
460:     VecDot(qp->G,qp->Z,gap);
461:     VecDot(qp->T,qp->S,gap+1);

463:     /* Compute the duality gap */
464:     qp->pobj = d1/2.0 + d2+qp->d;
465:     qp->gap  = gap[0]+gap[1];
466:     qp->dobj = qp->pobj - qp->gap;
467:     if (qp->m > 0) {
468:       qp->mu = qp->gap/(qp->m);
469:     }
470:     qp->rgap = qp->gap/(PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
471:   }  /* END MAIN LOOP  */
472:   return(0);
473: }

475: static PetscErrorCode TaoView_BQPIP(Tao tao,PetscViewer viewer)
476: {
478:   return(0);
479: }

481: static PetscErrorCode TaoSetFromOptions_BQPIP(PetscOptionItems *PetscOptionsObject,Tao tao)
482: {
483:   TAO_BQPIP      *qp = (TAO_BQPIP*)tao->data;

487:   PetscOptionsHead(PetscOptionsObject,"Interior point method for bound constrained quadratic optimization");
488:   PetscOptionsInt("-tao_bqpip_predcorr","Use a predictor-corrector method","",qp->predcorr,&qp->predcorr,NULL);
489:   PetscOptionsTail();
490:   KSPSetFromOptions(tao->ksp);
491:   return(0);
492: }

494: static PetscErrorCode TaoDestroy_BQPIP(Tao tao)
495: {
496:   TAO_BQPIP      *qp = (TAO_BQPIP*)tao->data;

500:   if (tao->setupcalled) {
501:     VecDestroy(&qp->G);
502:     VecDestroy(&qp->DG);
503:     VecDestroy(&qp->Z);
504:     VecDestroy(&qp->DZ);
505:     VecDestroy(&qp->GZwork);
506:     VecDestroy(&qp->R3);
507:     VecDestroy(&qp->S);
508:     VecDestroy(&qp->DS);
509:     VecDestroy(&qp->T);

511:     VecDestroy(&qp->DT);
512:     VecDestroy(&qp->TSwork);
513:     VecDestroy(&qp->R5);
514:     VecDestroy(&qp->HDiag);
515:     VecDestroy(&qp->Work);
516:     VecDestroy(&qp->XL);
517:     VecDestroy(&qp->XU);
518:     VecDestroy(&qp->DiagAxpy);
519:     VecDestroy(&qp->RHS);
520:     VecDestroy(&qp->RHS2);
521:     VecDestroy(&qp->C);
522:   }
523:   PetscFree(tao->data);
524:   return(0);
525: }

527: static PetscErrorCode TaoComputeDual_BQPIP(Tao tao,Vec DXL,Vec DXU)
528: {
529:   TAO_BQPIP       *qp = (TAO_BQPIP*)tao->data;
530:   PetscErrorCode  ierr;

533:   VecCopy(qp->Z,DXL);
534:   VecCopy(qp->S,DXU);
535:   VecScale(DXU,-1.0);
536:   return(0);
537: }

539: /*MC
540:  TAOBQPIP - interior-point method for quadratic programs with
541:     box constraints.

543:  Notes:
544:     This algorithms solves quadratic problems only, the Hessian will
545:         only be computed once.

547:  Options Database Keys:
548: . -tao_bqpip_predcorr - use a predictor/corrector method

550:   Level: beginner
551: M*/

553: PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao)
554: {
555:   TAO_BQPIP      *qp;

559:   PetscNewLog(tao,&qp);

561:   tao->ops->setup = TaoSetup_BQPIP;
562:   tao->ops->solve = TaoSolve_BQPIP;
563:   tao->ops->view = TaoView_BQPIP;
564:   tao->ops->setfromoptions = TaoSetFromOptions_BQPIP;
565:   tao->ops->destroy = TaoDestroy_BQPIP;
566:   tao->ops->computedual = TaoComputeDual_BQPIP;

568:   /* Override default settings (unless already changed) */
569:   if (!tao->max_it_changed) tao->max_it=100;
570:   if (!tao->max_funcs_changed) tao->max_funcs = 500;
571: #if defined(PETSC_USE_REAL_SINGLE)
572:   if (!tao->catol_changed) tao->catol=1e-6;
573: #else
574:   if (!tao->catol_changed) tao->catol=1e-12;
575: #endif

577:   /* Initialize pointers and variables */
578:   qp->n = 0;
579:   qp->m = 0;

581:   qp->predcorr = 1;
582:   tao->data    = (void*)qp;

584:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
585:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
586:   KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
587:   KSPSetType(tao->ksp,KSPCG);
588:   KSPSetTolerances(tao->ksp,1e-14,1e-30,1e30,PetscMax(10,qp->n));
589:   return(0);
590: }