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: {
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: }