Actual source code: bqpip.c

  1: #include "bqpip.h"
  2: #include "petscksp.h"

  6: static PetscErrorCode TaoSetUp_BQPIP(TaoSolver tao)
  7: {
  8:   TAO_BQPIP *qp =(TAO_BQPIP*)tao->data;

 12:   /* Set pointers to Data */
 13:   VecGetSize(tao->solution,&qp->n); 
 14:   KSPSetTolerances(tao->ksp, PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, qp->n); 

 16:   /* Allocate some arrays */
 17:   if (!tao->gradient) {
 18:       VecDuplicate(tao->solution, &tao->gradient);
 19:       
 20:   }
 21:   if (!tao->stepdirection) {
 22:       VecDuplicate(tao->solution, &tao->stepdirection); 
 23:       
 24:   }
 25:   if (!tao->XL) {
 26:       VecDuplicate(tao->solution, &tao->XL); 
 27:       VecSet(tao->XL, -1.0e-20); 
 28:   }
 29:   if (!tao->XU) {
 30:       VecDuplicate(tao->solution, &tao->XU); 
 31:       VecSet(tao->XU, 1.0e20); 
 32:   }


 35:   VecDuplicate(tao->solution, &qp->Work); 
 36:   VecDuplicate(tao->solution, &qp->XU); 
 37:   VecDuplicate(tao->solution, &qp->XL); 
 38:   VecDuplicate(tao->solution, &qp->HDiag); 
 39:   VecDuplicate(tao->solution, &qp->DiagAxpy); 
 40:   VecDuplicate(tao->solution, &qp->RHS); 
 41:   VecDuplicate(tao->solution, &qp->RHS2); 
 42:   VecDuplicate(tao->solution, &qp->C0); 

 44:   VecDuplicate(tao->solution, &qp->G); 
 45:   VecDuplicate(tao->solution, &qp->DG); 
 46:   VecDuplicate(tao->solution, &qp->S); 
 47:   VecDuplicate(tao->solution, &qp->Z); 
 48:   VecDuplicate(tao->solution, &qp->DZ); 
 49:   VecDuplicate(tao->solution, &qp->GZwork); 
 50:   VecDuplicate(tao->solution, &qp->R3); 

 52:   VecDuplicate(tao->solution, &qp->T); 
 53:   VecDuplicate(tao->solution, &qp->DT); 
 54:   VecDuplicate(tao->solution, &qp->DS); 
 55:   VecDuplicate(tao->solution, &qp->TSwork); 
 56:   VecDuplicate(tao->solution, &qp->R5); 


 59:   qp->m=2*qp->n;

 61:   return(0);
 62: }

 66: static PetscErrorCode  QPIPSetInitialPoint(TAO_BQPIP *qp, TaoSolver tao)
 67: {
 68:   PetscErrorCode       ierr;
 69:   PetscReal    two=2.0,p01=1;
 70:   PetscReal    gap1,gap2,fff,mu;

 73:   /* Compute function, Gradient R=Hx+b, and Hessian */
 74:   TaoComputeVariableBounds(tao); 
 75:   VecMedian(qp->XL, tao->solution, qp->XU, tao->solution); 
 76:   MatMult(tao->hessian, tao->solution, tao->gradient); 
 77:   VecCopy(qp->C0, qp->Work); 
 78:   VecAXPY(qp->Work, 0.5, tao->gradient); 
 79:   VecAXPY(tao->gradient, 1.0, qp->C0); 
 80:   VecDot(tao->solution, qp->Work, &fff); 
 81:   qp->pobj = fff + qp->c;

 83:   /* Initialize Primal Vectors */
 84:   /* T = XU - X; G = X - XL */
 85:   VecCopy(qp->XU, qp->T); 
 86:   VecAXPY(qp->T, -1.0, tao->solution); 
 87:   VecCopy(tao->solution, qp->G); 
 88:   VecAXPY(qp->G, -1.0, qp->XL); 

 90:   VecSet(qp->GZwork, p01); 
 91:   VecSet(qp->TSwork, p01); 

 93:   VecPointwiseMax(qp->G, qp->G, qp->GZwork); 
 94:   VecPointwiseMax(qp->T, qp->T, qp->TSwork); 

 96:   /* Initialize Dual Variable Vectors */
 97:   VecCopy(qp->G, qp->Z); 
 98:   VecReciprocal(qp->Z); 

100:   VecCopy(qp->T, qp->S); 
101:   VecReciprocal(qp->S); 

103:   MatMult(tao->hessian, qp->Work, qp->RHS); 
104:   VecAbs(qp->RHS); 
105:   VecSet(qp->Work, p01); 
106:   VecPointwiseMax(qp->RHS, qp->RHS, qp->Work); 

108:   VecPointwiseDivide(qp->RHS, tao->gradient, qp->RHS); 
109:   VecNorm(qp->RHS, NORM_1, &gap1); 
110:   mu = PetscMin(10.0,(gap1+10.0)/qp->m);

112:   VecScale(qp->S, mu); 
113:   VecScale(qp->Z, mu); 

115:   VecSet(qp->TSwork, p01); 
116:   VecSet(qp->GZwork, p01); 
117:   VecPointwiseMax(qp->S, qp->S, qp->TSwork); 
118:   VecPointwiseMax(qp->Z, qp->Z, qp->GZwork); 

120:   qp->mu=0;qp->dinfeas=1.0;qp->pinfeas=1.0;
121:   while ( (qp->dinfeas+qp->pinfeas)/(qp->m+qp->n) >= qp->mu ){

123:     VecScale(qp->G, two); 
124:     VecScale(qp->Z, two); 
125:     VecScale(qp->S, two); 
126:     VecScale(qp->T, two); 

128:     QPIPComputeResidual(qp,tao); 

130:     VecCopy(tao->solution, qp->R3); 
131:     VecAXPY(qp->R3, -1.0, qp->G); 
132:     VecAXPY(qp->R3, -1.0, qp->XL); 

134:     VecCopy(tao->solution, qp->R5); 
135:     VecAXPY(qp->R5, 1.0, qp->T); 
136:     VecAXPY(qp->R5, -1.0, qp->XU); 

138:     VecNorm(qp->R3, NORM_INFINITY, &gap1); 
139:     VecNorm(qp->R5, NORM_INFINITY, &gap2); 
140:     qp->pinfeas=PetscMax(gap1,gap2);
141:     
142:     /* Compute the duality gap */
143:     VecDot(qp->G, qp->Z, &gap1); 
144:     VecDot(qp->T, qp->S, &gap2); 
145:     
146:     qp->gap = (gap1+gap2);
147:     qp->dobj = qp->pobj - qp->gap;
148:     if (qp->m>0) qp->mu=qp->gap/(qp->m); else qp->mu=0.0;
149:     qp->rgap=qp->gap/( PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0 );
150:   }

152:   return(0);
153: }


158: static PetscErrorCode TaoDestroy_BQPIP(TaoSolver tao)
159: {
160:   TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;

163:   /* Free allocated memory in BQPIP structure */
165:   if (tao->setupcalled) {
166:     VecDestroy(&qp->G); 
167:     VecDestroy(&qp->DG); 
168:     VecDestroy(&qp->Z); 
169:     VecDestroy(&qp->DZ); 
170:     VecDestroy(&qp->GZwork); 
171:     VecDestroy(&qp->R3); 
172:     VecDestroy(&qp->S); 
173:     VecDestroy(&qp->DS); 
174:     VecDestroy(&qp->T); 

176:     VecDestroy(&qp->DT); 
177:     VecDestroy(&qp->TSwork); 
178:     VecDestroy(&qp->R5); 
179:     VecDestroy(&qp->HDiag); 
180:     VecDestroy(&qp->Work); 
181:     VecDestroy(&qp->XL); 
182:     VecDestroy(&qp->XU); 
183:     VecDestroy(&qp->DiagAxpy); 
184:     VecDestroy(&qp->RHS); 
185:     VecDestroy(&qp->RHS2); 
186:     VecDestroy(&qp->C0); 
187:   }
188:   PetscFree(tao->data); 
189:   tao->data = PETSC_NULL;

191:   return(0);
192: }

196: static PetscErrorCode TaoSolve_BQPIP(TaoSolver tao)
197: {
198:   TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
200:   PetscInt    iter=0,its;
201:   PetscReal    d1,d2,ksptol,sigma;
202:   PetscReal    sigmamu;
203:   PetscReal    dstep,pstep,step=0;
204:   PetscReal    gap[4];
205:   TaoSolverTerminationReason reason;
206:   MatStructure matflag;
207:   
209:   qp->dobj           = 0.0;
210:   qp->pobj           = 1.0;
211:   qp->gap            = 10.0;
212:   qp->rgap           = 1.0;
213:   qp->mu             = 1.0;
214:   qp->sigma          = 1.0;
215:   qp->dinfeas        = 1.0;
216:   qp->psteplength    = 0.0;
217:   qp->dsteplength    = 0.0;


220:   /* Tighten infinite bounds, things break when we don't do this
221:     -- see test_bqpip.c 
222:   */
223:   VecSet(qp->XU,1.0e20); 
224:   VecSet(qp->XL,-1.0e20); 
225:   VecPointwiseMax(qp->XL,qp->XL,tao->XL); 
226:   VecPointwiseMin(qp->XU,qp->XU,tao->XU); 

228:   TaoComputeObjectiveAndGradient(tao,tao->solution,&qp->c,qp->C0); 
229:   TaoComputeHessian(tao,tao->solution,&tao->hessian,&tao->hessian_pre,&matflag); 
230:   MatMult(tao->hessian, tao->solution, qp->Work); 
231:   VecDot(tao->solution, qp->Work, &d1); 
232:   VecAXPY(qp->C0, -1.0, qp->Work); 
233:   VecDot(qp->C0, tao->solution, &d2); 
234:   qp->c -= (d1/2.0+d2);
235:   MatGetDiagonal(tao->hessian, qp->HDiag); 

237:   QPIPSetInitialPoint(qp,tao); 
238:   QPIPComputeResidual(qp,tao); 
239:   
240:   /* Enter main loop */
241:   while (1){

243:     /* Check Stopping Condition      */
244:     TaoMonitor(tao,iter++,qp->pobj,PetscSqrtScalar(qp->gap + qp->dinfeas),
245:                             qp->pinfeas, step, &reason); 
246:     if (reason != TAO_CONTINUE_ITERATING) break;

248:     /* 
249:        Dual Infeasibility Direction should already be in the right
250:        hand side from computing the residuals 
251:     */

253:     QPIPComputeNormFromCentralPath(qp,&d1); 

255:     if (iter > 0 && (qp->rnorm>5*qp->mu || d1*d1>qp->m*qp->mu*qp->mu) ) {
256:       sigma=1.0;sigmamu=qp->mu;
257:       sigma=0.0;sigmamu=0;
258:     } else {
259:       sigma=0.0;sigmamu=0;
260:     }
261:     VecSet(qp->DZ, sigmamu); 
262:     VecSet(qp->DS, sigmamu); 

264:     if (sigmamu !=0){
265:       VecPointwiseDivide(qp->DZ, qp->DZ, qp->G); 
266:       VecPointwiseDivide(qp->DS, qp->DS, qp->T); 
267:       VecCopy(qp->DZ,qp->RHS2); 
268:       VecAXPY(qp->RHS2, 1.0, qp->DS); 
269:     } else {
270:       VecZeroEntries(qp->RHS2); 
271:     }


274:     /* 
275:        Compute the Primal Infeasiblitiy RHS and the 
276:        Diagonal Matrix to be added to H and store in Work 
277:     */
278:     VecPointwiseDivide(qp->DiagAxpy, qp->Z, qp->G); 
279:     VecPointwiseMult(qp->GZwork, qp->DiagAxpy, qp->R3); 
280:     VecAXPY(qp->RHS, -1.0, qp->GZwork); 

282:     VecPointwiseDivide(qp->TSwork, qp->S, qp->T); 
283:     VecAXPY(qp->DiagAxpy, 1.0, qp->TSwork); 
284:     VecPointwiseMult(qp->TSwork, qp->TSwork, qp->R5); 
285:     VecAXPY(qp->RHS, -1.0, qp->TSwork); 
286:     VecAXPY(qp->RHS2, 1.0, qp->RHS); 

288:     /*  Determine the solving tolerance */
289:     ksptol = qp->mu/10.0;
290:     ksptol = PetscMin(ksptol,0.001);

292:     MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES); 
293:     
294:     KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre, matflag); 
295:     KSPSolve(tao->ksp, qp->RHS, tao->stepdirection); 
296:     KSPGetIterationNumber(tao->ksp,&its); 
297:     tao->ksp_its+=its;

299:     VecScale(qp->DiagAxpy, -1.0); 
300:     MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES); 
301:     VecScale(qp->DiagAxpy, -1.0); 
302:     QPComputeStepDirection(qp,tao); 
303:     QPStepLength(qp);  

305:     /* Calculate New Residual R1 in Work vector */
306:     MatMult(tao->hessian, tao->stepdirection, qp->RHS2); 
307:     VecAXPY(qp->RHS2, 1.0, qp->DS); 
308:     VecAXPY(qp->RHS2, -1.0, qp->DZ); 
309:     VecAYPX(qp->RHS2, qp->dsteplength, tao->gradient); 

311:     VecNorm(qp->RHS2, NORM_2, &qp->dinfeas); 
312:     VecDot(qp->DZ, qp->DG, gap); 
313:     VecDot(qp->DS, qp->DT, gap+1); 

315:     qp->rnorm=(qp->dinfeas+qp->psteplength*qp->pinfeas)/(qp->m+qp->n);
316:     pstep = qp->psteplength; dstep = qp->dsteplength;    
317:     step = PetscMin(qp->psteplength,qp->dsteplength);
318:     sigmamu= ( pstep*pstep*(gap[0]+gap[1]) +
319:                (1 - pstep + pstep*sigma)*qp->gap  )/qp->m;

321:     if (qp->predcorr && step < 0.9){
322:       if (sigmamu < qp->mu){ 
323:         sigmamu=sigmamu/qp->mu;
324:         sigmamu=sigmamu*sigmamu*sigmamu;
325:       } else {sigmamu = 1.0;}
326:       sigmamu = sigmamu*qp->mu;
327:       
328:       /* Compute Corrector Step */
329:       VecPointwiseMult(qp->DZ, qp->DG, qp->DZ); 
330:       VecScale(qp->DZ, -1.0); 
331:       VecShift(qp->DZ, sigmamu); 
332:       VecPointwiseDivide(qp->DZ, qp->DZ, qp->G); 
333:       
334:       VecPointwiseMult(qp->DS, qp->DS, qp->DT); 
335:       VecScale(qp->DS, -1.0); 
336:       VecShift(qp->DS, sigmamu); 
337:       VecPointwiseDivide(qp->DS, qp->DS, qp->T); 

339:       VecCopy(qp->DZ, qp->RHS2); 
340:       VecAXPY(qp->RHS2, -1.0, qp->DS); 
341:       VecAXPY(qp->RHS2, 1.0, qp->RHS); 

343:       /* Approximately solve the linear system */
344:       MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES); 
345:       KSPSolve(tao->ksp, qp->RHS2, tao->stepdirection); 
346:       KSPGetIterationNumber(tao->ksp,&its); 
347:       tao->ksp_its+=its;
348:       
349:       MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES); 
350:       QPComputeStepDirection(qp,tao); 
351:       QPStepLength(qp); 

353:     }  /* End Corrector step */


356:     /* Take the step */
357:     pstep = qp->psteplength; dstep = qp->dsteplength;

359:     VecAXPY(qp->Z, dstep, qp->DZ); 
360:     VecAXPY(qp->S, dstep, qp->DS); 
361:     VecAXPY(tao->solution, dstep, tao->stepdirection); 
362:     VecAXPY(qp->G, dstep, qp->DG); 
363:     VecAXPY(qp->T, dstep, qp->DT); 

365:     /* Compute Residuals */
366:     QPIPComputeResidual(qp,tao); 

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

371:     VecDot(tao->solution, qp->Work, &d1); 
372:     VecDot(tao->solution, qp->C0, &d2); 
373:     VecDot(qp->G, qp->Z, gap); 
374:     VecDot(qp->T, qp->S, gap+1); 

376:     qp->pobj=d1/2.0 + d2+qp->c;
377:     /* Compute the duality gap */
378:     qp->gap = (gap[0]+gap[1]);
379:     qp->dobj = qp->pobj - qp->gap;
380:     if (qp->m>0) qp->mu=qp->gap/(qp->m);
381:     qp->rgap=qp->gap/( PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0 );
382:   }  /* END MAIN LOOP  */

384:   return(0);
385: }

389: static PetscErrorCode QPComputeStepDirection(TAO_BQPIP *qp, TaoSolver tao)
390: {


395:   /* Calculate DG */
396:   VecCopy(tao->stepdirection, qp->DG); 
397:   VecAXPY(qp->DG, 1.0, qp->R3); 

399:   /* Calculate DT */
400:   VecCopy(tao->stepdirection, qp->DT); 
401:   VecScale(qp->DT, -1.0); 
402:   VecAXPY(qp->DT, -1.0, qp->R5); 


405:   /* Calculate DZ */
406:   VecAXPY(qp->DZ, -1.0, qp->Z); 
407:   VecPointwiseDivide(qp->GZwork, qp->DG, qp->G); 
408:   VecPointwiseMult(qp->GZwork, qp->GZwork, qp->Z); 
409:   VecAXPY(qp->DZ, -1.0, qp->GZwork); 

411:   /* Calculate DS */
412:   VecAXPY(qp->DS, -1.0, qp->S); 
413:   VecPointwiseDivide(qp->TSwork, qp->DT, qp->T); 
414:   VecPointwiseMult(qp->TSwork, qp->TSwork, qp->S); 
415:   VecAXPY(qp->DS, -1.0, qp->TSwork); 

417:   return(0);
418: }

422: static PetscErrorCode QPIPComputeResidual(TAO_BQPIP *qp, TaoSolver tao)
423: {
425:   PetscReal dtmp = 1.0 - qp->psteplength;


429:   /* Compute R3 and R5 */

431:   VecScale(qp->R3, dtmp); 
432:   VecScale(qp->R5, dtmp); 
433:   qp->pinfeas=dtmp*qp->pinfeas;


436:   VecCopy(qp->S, tao->gradient); 
437:   VecAXPY(tao->gradient, -1.0, qp->Z); 
438:   
439:   MatMult(tao->hessian, tao->solution, qp->RHS); 
440:   VecScale(qp->RHS, -1.0); 
441:   VecAXPY(qp->RHS, -1.0, qp->C0); 
442:   VecAXPY(tao->gradient, -1.0, qp->RHS); 

444:   VecNorm(tao->gradient, NORM_1, &qp->dinfeas); 
445:   qp->rnorm=(qp->dinfeas+qp->pinfeas)/(qp->m+qp->n);
446:   
447:   return(0);
448: }

452: static PetscErrorCode QPStepLength(TAO_BQPIP *qp)
453: {
454:   PetscReal tstep1,tstep2,tstep3,tstep4,tstep;

458:   /* Compute stepsize to the boundary */
459:   VecStepMax(qp->G, qp->DG, &tstep1); 
460:   VecStepMax(qp->T, qp->DT, &tstep2); 
461:   VecStepMax(qp->S, qp->DS, &tstep3); 
462:   VecStepMax(qp->Z, qp->DZ, &tstep4); 

464:   
465:   tstep = PetscMin(tstep1,tstep2);
466:   qp->psteplength = PetscMin(0.95*tstep,1.0);

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

471:   qp->psteplength = PetscMin(qp->psteplength,qp->dsteplength);
472:   qp->dsteplength = qp->psteplength;

474:   return(0);
475: }


480: PetscErrorCode TaoComputeDual_BQPIP(TaoSolver tao,Vec DXL, Vec DXU)
481: {
482:   TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
483:   PetscErrorCode       ierr;

486:   VecCopy(qp->Z, DXL); 
487:   VecCopy(qp->S, DXU); 
488:   VecScale(DXU, -1.0); 
489:   return(0);
490: }


495: PetscErrorCode QPIPComputeNormFromCentralPath(TAO_BQPIP *qp, PetscReal *norm)
496: {
497:   PetscErrorCode       ierr;
498:   PetscReal    gap[2],mu[2], nmu;
499:   
501:   VecPointwiseMult(qp->GZwork, qp->G, qp->Z); 
502:   VecPointwiseMult(qp->TSwork, qp->T, qp->S); 
503:   VecNorm(qp->TSwork, NORM_1, &mu[0]); 
504:   VecNorm(qp->GZwork, NORM_1, &mu[1]); 

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

508:   VecShift(qp->GZwork,nmu); 
509:   VecShift(qp->TSwork,nmu); 

511:   VecNorm(qp->GZwork,NORM_2,&gap[0]); 
512:   VecNorm(qp->TSwork,NORM_2,&gap[1]); 
513:   gap[0]*=gap[0];
514:   gap[1]*=gap[1];


517:   qp->pathnorm=PetscSqrtScalar( (gap[0]+gap[1]) );
518:   *norm=qp->pathnorm;

520:   return(0);
521: }


526: static PetscErrorCode TaoSetFromOptions_BQPIP(TaoSolver tao) 
527: {
528:   TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
529:   PetscErrorCode       ierr;

532:   PetscOptionsHead("Interior point method for bound constrained quadratic optimization");
533:   PetscOptionsInt("-predcorr","Use a predictor-corrector method","",qp->predcorr,&qp->predcorr,0);
534:   
535:   PetscOptionsTail();
536:   KSPSetFromOptions(tao->ksp); 
537:   return(0);
538: }


543: static PetscErrorCode TaoView_BQPIP(TaoSolver tao, PetscViewer viewer)
544: {
546:   return(0);
547: }


550: /* --------------------------------------------------------- */
554: PetscErrorCode TaoCreate_BQPIP(TaoSolver tao)
555: {
556:   TAO_BQPIP *qp;
557:   PetscErrorCode       ierr;

560:   PetscNewLog(tao, TAO_BQPIP, &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:   tao->max_it=100;
569:   tao->max_funcs = 500;
570:   tao->fatol=1e-12;
571:   tao->frtol=1e-12;
572:   tao->gatol=1e-12;
573:   tao->grtol=1e-12;
574:   tao->catol=1e-12;


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

582:   qp->predcorr       = 1;


585:   tao->data = (void*)qp;

587:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp); 
588:   KSPSetType(tao->ksp, KSPCG); 
589:   KSPSetTolerances(tao->ksp, 1e-14, 1e-30, 1e30, qp->n); 


592:   return(0);
593: }