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