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;qp->dinfeas=1.0;qp->pinfeas=1.0;
 86:   while ((qp->dinfeas+qp->pinfeas)/(qp->m+qp->n) >= qp->mu) {
 87:     VecScale(qp->G,two);
 88:     VecScale(qp->Z,two);
 89:     VecScale(qp->S,two);
 90:     VecScale(qp->T,two);

 92:     QPIPComputeResidual(qp,tao);

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

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

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

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

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

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

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

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

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

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

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

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

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

154:   VecShift(qp->GZwork,nmu);
155:   VecShift(qp->TSwork,nmu);

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

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

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

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

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

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

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

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

199:   /* Allocate some arrays */
200:   if (!tao->gradient) {
201:     VecDuplicate(tao->solution,&tao->gradient);
202:   }
203:   if (!tao->stepdirection) {
204:     VecDuplicate(tao->solution,&tao->stepdirection);
205:   }
206:   if (!tao->XL) {
207:     VecDuplicate(tao->solution,&tao->XL);
208:     VecSet(tao->XL,PETSC_NINFINITY);
209:   }
210:   if (!tao->XU) {
211:     VecDuplicate(tao->solution,&tao->XU);
212:     VecSet(tao->XU,PETSC_INFINITY);
213:   }

215:   VecDuplicate(tao->solution,&qp->Work);
216:   VecDuplicate(tao->solution,&qp->XU);
217:   VecDuplicate(tao->solution,&qp->XL);
218:   VecDuplicate(tao->solution,&qp->HDiag);
219:   VecDuplicate(tao->solution,&qp->DiagAxpy);
220:   VecDuplicate(tao->solution,&qp->RHS);
221:   VecDuplicate(tao->solution,&qp->RHS2);
222:   VecDuplicate(tao->solution,&qp->C);

224:   VecDuplicate(tao->solution,&qp->G);
225:   VecDuplicate(tao->solution,&qp->DG);
226:   VecDuplicate(tao->solution,&qp->S);
227:   VecDuplicate(tao->solution,&qp->Z);
228:   VecDuplicate(tao->solution,&qp->DZ);
229:   VecDuplicate(tao->solution,&qp->GZwork);
230:   VecDuplicate(tao->solution,&qp->R3);

232:   VecDuplicate(tao->solution,&qp->T);
233:   VecDuplicate(tao->solution,&qp->DT);
234:   VecDuplicate(tao->solution,&qp->DS);
235:   VecDuplicate(tao->solution,&qp->TSwork);
236:   VecDuplicate(tao->solution,&qp->R5);
237:   qp->m=2*qp->n;
238:   return 0;
239: }

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

250:   qp->dobj        = 0.0;
251:   qp->pobj        = 1.0;
252:   qp->gap         = 10.0;
253:   qp->rgap        = 1.0;
254:   qp->mu          = 1.0;
255:   qp->dinfeas     = 1.0;
256:   qp->psteplength = 0.0;
257:   qp->dsteplength = 0.0;

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

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

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

287:   /* Initialize starting point and residuals */
288:   QPIPSetInitialPoint(qp,tao);
289:   QPIPComputeResidual(qp,tao);

291:   /* Enter main loop */
292:   tao->reason = TAO_CONTINUE_ITERATING;
293:   while (1) {

295:     /* Check Stopping Condition      */
296:     gnorm = PetscSqrtScalar(qp->gap + qp->dinfeas);
297:     TaoLogConvergenceHistory(tao,qp->pobj,gnorm,qp->pinfeas,tao->ksp_its);
298:     TaoMonitor(tao,tao->niter,qp->pobj,gnorm,qp->pinfeas,step);
299:     (*tao->ops->convergencetest)(tao,tao->cnvP);
300:     if (tao->reason != TAO_CONTINUE_ITERATING) break;
301:     /* Call general purpose update function */
302:     if (tao->ops->update) {
303:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
304:     }
305:     tao->niter++;
306:     tao->ksp_its = 0;

308:     /*
309:        Dual Infeasibility Direction should already be in the right
310:        hand side from computing the residuals
311:     */

313:     QPIPComputeNormFromCentralPath(qp,&d1);

315:     VecSet(qp->DZ,0.0);
316:     VecSet(qp->DS,0.0);

318:     /*
319:        Compute the Primal Infeasiblitiy RHS and the
320:        Diagonal Matrix to be added to H and store in Work
321:     */
322:     VecPointwiseDivide(qp->DiagAxpy,qp->Z,qp->G);
323:     VecPointwiseMult(qp->GZwork,qp->DiagAxpy,qp->R3);
324:     VecAXPY(qp->RHS,-1.0,qp->GZwork);

326:     VecPointwiseDivide(qp->TSwork,qp->S,qp->T);
327:     VecAXPY(qp->DiagAxpy,1.0,qp->TSwork);
328:     VecPointwiseMult(qp->TSwork,qp->TSwork,qp->R5);
329:     VecAXPY(qp->RHS,-1.0,qp->TSwork);

331:     /*  Determine the solving tolerance */
332:     ksptol = qp->mu/10.0;
333:     ksptol = PetscMin(ksptol,0.001);
334:     KSPSetTolerances(tao->ksp,ksptol,1e-30,1e30,PetscMax(10,qp->n));

336:     /* Shift the diagonals of the Hessian matrix */
337:     MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES);
338:     if (!getdiagop) {
339:       VecCopy(qp->DiagAxpy,qp->HDiag);
340:       VecScale(qp->HDiag,-1.0);
341:     }
342:     MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
343:     MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);

345:     KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre);
346:     KSPSolve(tao->ksp,qp->RHS,tao->stepdirection);
347:     KSPGetIterationNumber(tao->ksp,&its);
348:     tao->ksp_its += its;
349:     tao->ksp_tot_its += its;

351:     /* Restore the true diagonal of the Hessian matrix */
352:     if (getdiagop) {
353:       MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES);
354:     } else {
355:       MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES);
356:     }
357:     MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
358:     MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
359:     QPIPComputeStepDirection(qp,tao);
360:     QPIPStepLength(qp);

362:     /* Calculate New Residual R1 in Work vector */
363:     MatMult(tao->hessian,tao->stepdirection,qp->RHS2);
364:     VecAXPY(qp->RHS2,1.0,qp->DS);
365:     VecAXPY(qp->RHS2,-1.0,qp->DZ);
366:     VecAYPX(qp->RHS2,qp->dsteplength,tao->gradient);

368:     VecNorm(qp->RHS2,NORM_2,&qp->dinfeas);
369:     VecDot(qp->DZ,qp->DG,gap);
370:     VecDot(qp->DS,qp->DT,gap+1);

372:     qp->rnorm = (qp->dinfeas+qp->psteplength*qp->pinfeas)/(qp->m+qp->n);
373:     pstep     = qp->psteplength;
374:     step      = PetscMin(qp->psteplength,qp->dsteplength);
375:     sigmamu   = (pstep*pstep*(gap[0]+gap[1]) + (1 - pstep)*qp->gap)/qp->m;

377:     if (qp->predcorr && step < 0.9) {
378:       if (sigmamu < qp->mu) {
379:         sigmamu = sigmamu/qp->mu;
380:         sigmamu = sigmamu*sigmamu*sigmamu;
381:       } else {
382:         sigmamu = 1.0;
383:       }
384:       sigmamu = sigmamu*qp->mu;

386:       /* Compute Corrector Step */
387:       VecPointwiseMult(qp->DZ,qp->DG,qp->DZ);
388:       VecScale(qp->DZ,-1.0);
389:       VecShift(qp->DZ,sigmamu);
390:       VecPointwiseDivide(qp->DZ,qp->DZ,qp->G);

392:       VecPointwiseMult(qp->DS,qp->DS,qp->DT);
393:       VecScale(qp->DS,-1.0);
394:       VecShift(qp->DS,sigmamu);
395:       VecPointwiseDivide(qp->DS,qp->DS,qp->T);

397:       VecCopy(qp->DZ,qp->RHS2);
398:       VecAXPY(qp->RHS2,-1.0,qp->DS);
399:       VecAXPY(qp->RHS2,1.0,qp->RHS);

401:       /* Approximately solve the linear system */
402:       MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES);
403:       if (!getdiagop) {
404:         VecCopy(qp->DiagAxpy,qp->HDiag);
405:         VecScale(qp->HDiag,-1.0);
406:       }
407:       MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
408:       MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);

410:       /* Solve using the previous tolerances that were set */
411:       KSPSolve(tao->ksp,qp->RHS2,tao->stepdirection);
412:       KSPGetIterationNumber(tao->ksp,&its);
413:       tao->ksp_its += its;
414:       tao->ksp_tot_its += its;

416:       if (getdiagop) {
417:         MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES);
418:       } else {
419:         MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES);
420:       }
421:       MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
422:       MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
423:       QPIPComputeStepDirection(qp,tao);
424:       QPIPStepLength(qp);
425:     }  /* End Corrector step */

427:     /* Take the step */
428:     dstep = qp->dsteplength;

430:     VecAXPY(qp->Z,dstep,qp->DZ);
431:     VecAXPY(qp->S,dstep,qp->DS);
432:     VecAXPY(tao->solution,dstep,tao->stepdirection);
433:     VecAXPY(qp->G,dstep,qp->DG);
434:     VecAXPY(qp->T,dstep,qp->DT);

436:     /* Compute Residuals */
437:     QPIPComputeResidual(qp,tao);

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

442:     VecDot(tao->solution,qp->Work,&d1);
443:     VecDot(tao->solution,qp->C,&d2);
444:     VecDot(qp->G,qp->Z,gap);
445:     VecDot(qp->T,qp->S,gap+1);

447:     /* Compute the duality gap */
448:     qp->pobj = d1/2.0 + d2+qp->d;
449:     qp->gap  = gap[0]+gap[1];
450:     qp->dobj = qp->pobj - qp->gap;
451:     if (qp->m > 0) {
452:       qp->mu = qp->gap/(qp->m);
453:     }
454:     qp->rgap = qp->gap/(PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
455:   }  /* END MAIN LOOP  */
456:   return 0;
457: }

459: static PetscErrorCode TaoView_BQPIP(Tao tao,PetscViewer viewer)
460: {
461:   return 0;
462: }

464: static PetscErrorCode TaoSetFromOptions_BQPIP(PetscOptionItems *PetscOptionsObject,Tao tao)
465: {
466:   TAO_BQPIP      *qp = (TAO_BQPIP*)tao->data;

468:   PetscOptionsHead(PetscOptionsObject,"Interior point method for bound constrained quadratic optimization");
469:   PetscOptionsInt("-tao_bqpip_predcorr","Use a predictor-corrector method","",qp->predcorr,&qp->predcorr,NULL);
470:   PetscOptionsTail();
471:   KSPSetFromOptions(tao->ksp);
472:   return 0;
473: }

475: static PetscErrorCode TaoDestroy_BQPIP(Tao tao)
476: {
477:   TAO_BQPIP      *qp = (TAO_BQPIP*)tao->data;

479:   if (tao->setupcalled) {
480:     VecDestroy(&qp->G);
481:     VecDestroy(&qp->DG);
482:     VecDestroy(&qp->Z);
483:     VecDestroy(&qp->DZ);
484:     VecDestroy(&qp->GZwork);
485:     VecDestroy(&qp->R3);
486:     VecDestroy(&qp->S);
487:     VecDestroy(&qp->DS);
488:     VecDestroy(&qp->T);

490:     VecDestroy(&qp->DT);
491:     VecDestroy(&qp->TSwork);
492:     VecDestroy(&qp->R5);
493:     VecDestroy(&qp->HDiag);
494:     VecDestroy(&qp->Work);
495:     VecDestroy(&qp->XL);
496:     VecDestroy(&qp->XU);
497:     VecDestroy(&qp->DiagAxpy);
498:     VecDestroy(&qp->RHS);
499:     VecDestroy(&qp->RHS2);
500:     VecDestroy(&qp->C);
501:   }
502:   PetscFree(tao->data);
503:   return 0;
504: }

506: static PetscErrorCode TaoComputeDual_BQPIP(Tao tao,Vec DXL,Vec DXU)
507: {
508:   TAO_BQPIP       *qp = (TAO_BQPIP*)tao->data;

510:   VecCopy(qp->Z,DXL);
511:   VecCopy(qp->S,DXU);
512:   VecScale(DXU,-1.0);
513:   return 0;
514: }

516: /*MC
517:  TAOBQPIP - interior-point method for quadratic programs with
518:     box constraints.

520:  Notes:
521:     This algorithms solves quadratic problems only, the Hessian will
522:         only be computed once.

524:  Options Database Keys:
525: . -tao_bqpip_predcorr - use a predictor/corrector method

527:   Level: beginner
528: M*/

530: PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao)
531: {
532:   TAO_BQPIP      *qp;

534:   PetscNewLog(tao,&qp);

536:   tao->ops->setup = TaoSetup_BQPIP;
537:   tao->ops->solve = TaoSolve_BQPIP;
538:   tao->ops->view = TaoView_BQPIP;
539:   tao->ops->setfromoptions = TaoSetFromOptions_BQPIP;
540:   tao->ops->destroy = TaoDestroy_BQPIP;
541:   tao->ops->computedual = TaoComputeDual_BQPIP;

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

552:   /* Initialize pointers and variables */
553:   qp->n = 0;
554:   qp->m = 0;

556:   qp->predcorr = 1;
557:   tao->data    = (void*)qp;

559:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
560:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
561:   KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
562:   KSPSetType(tao->ksp,KSPCG);
563:   KSPSetTolerances(tao->ksp,1e-14,1e-30,1e30,PetscMax(10,qp->n));
564:   return 0;
565: }