Actual source code: bqpip.c
petsc-3.6.4 2016-04-12
1: #include <../src/tao/bound/impls/bqpip/bqpip.h>
2: #include <petscksp.h>
6: static PetscErrorCode TaoSetUp_BQPIP(Tao tao)
7: {
8: TAO_BQPIP *qp =(TAO_BQPIP*)tao->data;
12: /* Set pointers to Data */
13: VecGetSize(tao->solution,&qp->n);
15: /* Allocate some arrays */
16: if (!tao->gradient) {
17: VecDuplicate(tao->solution, &tao->gradient);
18: }
19: if (!tao->stepdirection) {
20: VecDuplicate(tao->solution, &tao->stepdirection);
21: }
22: if (!tao->XL) {
23: VecDuplicate(tao->solution, &tao->XL);
24: VecSet(tao->XL, -1.0e-20);
25: }
26: if (!tao->XU) {
27: VecDuplicate(tao->solution, &tao->XU);
28: VecSet(tao->XU, 1.0e20);
29: }
31: VecDuplicate(tao->solution, &qp->Work);
32: VecDuplicate(tao->solution, &qp->XU);
33: VecDuplicate(tao->solution, &qp->XL);
34: VecDuplicate(tao->solution, &qp->HDiag);
35: VecDuplicate(tao->solution, &qp->DiagAxpy);
36: VecDuplicate(tao->solution, &qp->RHS);
37: VecDuplicate(tao->solution, &qp->RHS2);
38: VecDuplicate(tao->solution, &qp->C0);
40: VecDuplicate(tao->solution, &qp->G);
41: VecDuplicate(tao->solution, &qp->DG);
42: VecDuplicate(tao->solution, &qp->S);
43: VecDuplicate(tao->solution, &qp->Z);
44: VecDuplicate(tao->solution, &qp->DZ);
45: VecDuplicate(tao->solution, &qp->GZwork);
46: VecDuplicate(tao->solution, &qp->R3);
48: VecDuplicate(tao->solution, &qp->T);
49: VecDuplicate(tao->solution, &qp->DT);
50: VecDuplicate(tao->solution, &qp->DS);
51: VecDuplicate(tao->solution, &qp->TSwork);
52: VecDuplicate(tao->solution, &qp->R5);
53: qp->m=2*qp->n;
54: return(0);
55: }
59: static PetscErrorCode QPIPSetInitialPoint(TAO_BQPIP *qp, Tao tao)
60: {
62: PetscReal two=2.0,p01=1;
63: PetscReal gap1,gap2,fff,mu;
66: /* Compute function, Gradient R=Hx+b, and Hessian */
67: TaoComputeVariableBounds(tao);
68: VecMedian(qp->XL, tao->solution, qp->XU, tao->solution);
69: MatMult(tao->hessian, tao->solution, tao->gradient);
70: VecCopy(qp->C0, qp->Work);
71: VecAXPY(qp->Work, 0.5, tao->gradient);
72: VecAXPY(tao->gradient, 1.0, qp->C0);
73: VecDot(tao->solution, qp->Work, &fff);
74: qp->pobj = fff + qp->c;
76: /* Initialize Primal Vectors */
77: /* T = XU - X; G = X - XL */
78: VecCopy(qp->XU, qp->T);
79: VecAXPY(qp->T, -1.0, tao->solution);
80: VecCopy(tao->solution, qp->G);
81: VecAXPY(qp->G, -1.0, qp->XL);
83: VecSet(qp->GZwork, p01);
84: VecSet(qp->TSwork, p01);
86: VecPointwiseMax(qp->G, qp->G, qp->GZwork);
87: VecPointwiseMax(qp->T, qp->T, qp->TSwork);
89: /* Initialize Dual Variable Vectors */
90: VecCopy(qp->G, qp->Z);
91: VecReciprocal(qp->Z);
93: VecCopy(qp->T, qp->S);
94: VecReciprocal(qp->S);
96: MatMult(tao->hessian, qp->Work, qp->RHS);
97: VecAbs(qp->RHS);
98: VecSet(qp->Work, p01);
99: VecPointwiseMax(qp->RHS, qp->RHS, qp->Work);
101: VecPointwiseDivide(qp->RHS, tao->gradient, qp->RHS);
102: VecNorm(qp->RHS, NORM_1, &gap1);
103: mu = PetscMin(10.0,(gap1+10.0)/qp->m);
105: VecScale(qp->S, mu);
106: VecScale(qp->Z, mu);
108: VecSet(qp->TSwork, p01);
109: VecSet(qp->GZwork, p01);
110: VecPointwiseMax(qp->S, qp->S, qp->TSwork);
111: VecPointwiseMax(qp->Z, qp->Z, qp->GZwork);
113: qp->mu=0;qp->dinfeas=1.0;qp->pinfeas=1.0;
114: while ( (qp->dinfeas+qp->pinfeas)/(qp->m+qp->n) >= qp->mu ){
116: VecScale(qp->G, two);
117: VecScale(qp->Z, two);
118: VecScale(qp->S, two);
119: VecScale(qp->T, two);
121: QPIPComputeResidual(qp,tao);
123: VecCopy(tao->solution, qp->R3);
124: VecAXPY(qp->R3, -1.0, qp->G);
125: VecAXPY(qp->R3, -1.0, qp->XL);
127: VecCopy(tao->solution, qp->R5);
128: VecAXPY(qp->R5, 1.0, qp->T);
129: VecAXPY(qp->R5, -1.0, qp->XU);
131: VecNorm(qp->R3, NORM_INFINITY, &gap1);
132: VecNorm(qp->R5, NORM_INFINITY, &gap2);
133: qp->pinfeas=PetscMax(gap1,gap2);
135: /* Compute the duality gap */
136: VecDot(qp->G, qp->Z, &gap1);
137: VecDot(qp->T, qp->S, &gap2);
139: qp->gap = (gap1+gap2);
140: qp->dobj = qp->pobj - qp->gap;
141: if (qp->m>0) qp->mu=qp->gap/(qp->m); else qp->mu=0.0;
142: qp->rgap=qp->gap/( PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0 );
143: }
144: return(0);
145: }
149: static PetscErrorCode TaoDestroy_BQPIP(Tao tao)
150: {
151: TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
155: if (tao->setupcalled) {
156: VecDestroy(&qp->G);
157: VecDestroy(&qp->DG);
158: VecDestroy(&qp->Z);
159: VecDestroy(&qp->DZ);
160: VecDestroy(&qp->GZwork);
161: VecDestroy(&qp->R3);
162: VecDestroy(&qp->S);
163: VecDestroy(&qp->DS);
164: VecDestroy(&qp->T);
166: VecDestroy(&qp->DT);
167: VecDestroy(&qp->TSwork);
168: VecDestroy(&qp->R5);
169: VecDestroy(&qp->HDiag);
170: VecDestroy(&qp->Work);
171: VecDestroy(&qp->XL);
172: VecDestroy(&qp->XU);
173: VecDestroy(&qp->DiagAxpy);
174: VecDestroy(&qp->RHS);
175: VecDestroy(&qp->RHS2);
176: VecDestroy(&qp->C0);
177: }
178: PetscFree(tao->data);
179: return(0);
180: }
184: static PetscErrorCode TaoSolve_BQPIP(Tao tao)
185: {
186: TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
187: PetscErrorCode ierr;
188: PetscInt its;
189: PetscReal d1,d2,ksptol,sigma;
190: PetscReal sigmamu;
191: PetscReal dstep,pstep,step=0;
192: PetscReal gap[4];
193: TaoConvergedReason reason;
196: qp->dobj = 0.0;
197: qp->pobj = 1.0;
198: qp->gap = 10.0;
199: qp->rgap = 1.0;
200: qp->mu = 1.0;
201: qp->sigma = 1.0;
202: qp->dinfeas = 1.0;
203: qp->psteplength = 0.0;
204: qp->dsteplength = 0.0;
206: /* Tighten infinite bounds, things break when we don't do this
207: -- see test_bqpip.c
208: */
209: VecSet(qp->XU,1.0e20);
210: VecSet(qp->XL,-1.0e20);
211: VecPointwiseMax(qp->XL,qp->XL,tao->XL);
212: VecPointwiseMin(qp->XU,qp->XU,tao->XU);
214: TaoComputeObjectiveAndGradient(tao,tao->solution,&qp->c,qp->C0);
215: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
216: MatMult(tao->hessian, tao->solution, qp->Work);
217: VecDot(tao->solution, qp->Work, &d1);
218: VecAXPY(qp->C0, -1.0, qp->Work);
219: VecDot(qp->C0, tao->solution, &d2);
220: qp->c -= (d1/2.0+d2);
221: MatGetDiagonal(tao->hessian, qp->HDiag);
223: QPIPSetInitialPoint(qp,tao);
224: QPIPComputeResidual(qp,tao);
226: /* Enter main loop */
227: while (PETSC_TRUE){
229: /* Check Stopping Condition */
230: TaoMonitor(tao,tao->niter,qp->pobj,PetscSqrtScalar(qp->gap + qp->dinfeas),qp->pinfeas, step, &reason);
231: if (reason != TAO_CONTINUE_ITERATING) break;
232: tao->niter++;
233: tao->ksp_its=0;
235: /*
236: Dual Infeasibility Direction should already be in the right
237: hand side from computing the residuals
238: */
240: QPIPComputeNormFromCentralPath(qp,&d1);
242: if (tao->niter > 0 && (qp->rnorm>5*qp->mu || d1*d1>qp->m*qp->mu*qp->mu) ) {
243: sigma=1.0;sigmamu=qp->mu;
244: sigma=0.0;sigmamu=0;
245: } else {
246: sigma=0.0;sigmamu=0;
247: }
248: VecSet(qp->DZ, sigmamu);
249: VecSet(qp->DS, sigmamu);
251: if (sigmamu !=0){
252: VecPointwiseDivide(qp->DZ, qp->DZ, qp->G);
253: VecPointwiseDivide(qp->DS, qp->DS, qp->T);
254: VecCopy(qp->DZ,qp->RHS2);
255: VecAXPY(qp->RHS2, 1.0, qp->DS);
256: } else {
257: VecZeroEntries(qp->RHS2);
258: }
261: /*
262: Compute the Primal Infeasiblitiy RHS and the
263: Diagonal Matrix to be added to H and store in Work
264: */
265: VecPointwiseDivide(qp->DiagAxpy, qp->Z, qp->G);
266: VecPointwiseMult(qp->GZwork, qp->DiagAxpy, qp->R3);
267: VecAXPY(qp->RHS, -1.0, qp->GZwork);
269: VecPointwiseDivide(qp->TSwork, qp->S, qp->T);
270: VecAXPY(qp->DiagAxpy, 1.0, qp->TSwork);
271: VecPointwiseMult(qp->TSwork, qp->TSwork, qp->R5);
272: VecAXPY(qp->RHS, -1.0, qp->TSwork);
273: VecAXPY(qp->RHS2, 1.0, qp->RHS);
275: /* Determine the solving tolerance */
276: ksptol = qp->mu/10.0;
277: ksptol = PetscMin(ksptol,0.001);
279: MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES);
280: MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
281: MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
283: KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);
284: KSPSolve(tao->ksp, qp->RHS, tao->stepdirection);
285: KSPGetIterationNumber(tao->ksp,&its);
286: tao->ksp_its+=its;
287: tao->ksp_tot_its+=its;
289: VecScale(qp->DiagAxpy, -1.0);
290: MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES);
291: MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
292: MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
293: VecScale(qp->DiagAxpy, -1.0);
294: QPComputeStepDirection(qp,tao);
295: QPStepLength(qp);
297: /* Calculate New Residual R1 in Work vector */
298: MatMult(tao->hessian, tao->stepdirection, qp->RHS2);
299: VecAXPY(qp->RHS2, 1.0, qp->DS);
300: VecAXPY(qp->RHS2, -1.0, qp->DZ);
301: VecAYPX(qp->RHS2, qp->dsteplength, tao->gradient);
303: VecNorm(qp->RHS2, NORM_2, &qp->dinfeas);
304: VecDot(qp->DZ, qp->DG, gap);
305: VecDot(qp->DS, qp->DT, gap+1);
307: qp->rnorm=(qp->dinfeas+qp->psteplength*qp->pinfeas)/(qp->m+qp->n);
308: pstep = qp->psteplength; dstep = qp->dsteplength;
309: step = PetscMin(qp->psteplength,qp->dsteplength);
310: sigmamu= ( pstep*pstep*(gap[0]+gap[1]) +
311: (1 - pstep + pstep*sigma)*qp->gap )/qp->m;
313: if (qp->predcorr && step < 0.9){
314: if (sigmamu < qp->mu){
315: sigmamu=sigmamu/qp->mu;
316: sigmamu=sigmamu*sigmamu*sigmamu;
317: } else {sigmamu = 1.0;}
318: sigmamu = sigmamu*qp->mu;
320: /* Compute Corrector Step */
321: VecPointwiseMult(qp->DZ, qp->DG, qp->DZ);
322: VecScale(qp->DZ, -1.0);
323: VecShift(qp->DZ, sigmamu);
324: VecPointwiseDivide(qp->DZ, qp->DZ, qp->G);
326: VecPointwiseMult(qp->DS, qp->DS, qp->DT);
327: VecScale(qp->DS, -1.0);
328: VecShift(qp->DS, sigmamu);
329: VecPointwiseDivide(qp->DS, qp->DS, qp->T);
331: VecCopy(qp->DZ, qp->RHS2);
332: VecAXPY(qp->RHS2, -1.0, qp->DS);
333: VecAXPY(qp->RHS2, 1.0, qp->RHS);
335: /* Approximately solve the linear system */
336: MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES);
337: MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
338: MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
339: KSPSolve(tao->ksp, qp->RHS2, tao->stepdirection);
340: KSPGetIterationNumber(tao->ksp,&its);
341: tao->ksp_its+=its;
342: tao->ksp_tot_its+=its;
344: MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES);
345: MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
346: MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
347: QPComputeStepDirection(qp,tao);
348: QPStepLength(qp);
350: } /* End Corrector step */
353: /* Take the step */
354: pstep = qp->psteplength; dstep = qp->dsteplength;
356: VecAXPY(qp->Z, dstep, qp->DZ);
357: VecAXPY(qp->S, dstep, qp->DS);
358: VecAXPY(tao->solution, dstep, tao->stepdirection);
359: VecAXPY(qp->G, dstep, qp->DG);
360: VecAXPY(qp->T, dstep, qp->DT);
362: /* Compute Residuals */
363: QPIPComputeResidual(qp,tao);
365: /* Evaluate quadratic function */
366: MatMult(tao->hessian, tao->solution, qp->Work);
368: VecDot(tao->solution, qp->Work, &d1);
369: VecDot(tao->solution, qp->C0, &d2);
370: VecDot(qp->G, qp->Z, gap);
371: VecDot(qp->T, qp->S, gap+1);
373: qp->pobj=d1/2.0 + d2+qp->c;
374: /* Compute the duality gap */
375: qp->gap = (gap[0]+gap[1]);
376: qp->dobj = qp->pobj - qp->gap;
377: if (qp->m>0) qp->mu=qp->gap/(qp->m);
378: qp->rgap=qp->gap/( PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0 );
379: } /* END MAIN LOOP */
381: return(0);
382: }
386: static PetscErrorCode QPComputeStepDirection(TAO_BQPIP *qp, Tao tao)
387: {
392: /* Calculate DG */
393: VecCopy(tao->stepdirection, qp->DG);
394: VecAXPY(qp->DG, 1.0, qp->R3);
396: /* Calculate DT */
397: VecCopy(tao->stepdirection, qp->DT);
398: VecScale(qp->DT, -1.0);
399: VecAXPY(qp->DT, -1.0, qp->R5);
402: /* Calculate DZ */
403: VecAXPY(qp->DZ, -1.0, qp->Z);
404: VecPointwiseDivide(qp->GZwork, qp->DG, qp->G);
405: VecPointwiseMult(qp->GZwork, qp->GZwork, qp->Z);
406: VecAXPY(qp->DZ, -1.0, qp->GZwork);
408: /* Calculate DS */
409: VecAXPY(qp->DS, -1.0, qp->S);
410: VecPointwiseDivide(qp->TSwork, qp->DT, qp->T);
411: VecPointwiseMult(qp->TSwork, qp->TSwork, qp->S);
412: VecAXPY(qp->DS, -1.0, qp->TSwork);
414: return(0);
415: }
419: static PetscErrorCode QPIPComputeResidual(TAO_BQPIP *qp, Tao tao)
420: {
422: PetscReal dtmp = 1.0 - qp->psteplength;
426: /* Compute R3 and R5 */
428: VecScale(qp->R3, dtmp);
429: VecScale(qp->R5, dtmp);
430: qp->pinfeas=dtmp*qp->pinfeas;
433: VecCopy(qp->S, tao->gradient);
434: VecAXPY(tao->gradient, -1.0, qp->Z);
436: MatMult(tao->hessian, tao->solution, qp->RHS);
437: VecScale(qp->RHS, -1.0);
438: VecAXPY(qp->RHS, -1.0, qp->C0);
439: VecAXPY(tao->gradient, -1.0, qp->RHS);
441: VecNorm(tao->gradient, NORM_1, &qp->dinfeas);
442: qp->rnorm=(qp->dinfeas+qp->pinfeas)/(qp->m+qp->n);
444: return(0);
445: }
449: static PetscErrorCode QPStepLength(TAO_BQPIP *qp)
450: {
451: PetscReal tstep1,tstep2,tstep3,tstep4,tstep;
455: /* Compute stepsize to the boundary */
456: VecStepMax(qp->G, qp->DG, &tstep1);
457: VecStepMax(qp->T, qp->DT, &tstep2);
458: VecStepMax(qp->S, qp->DS, &tstep3);
459: VecStepMax(qp->Z, qp->DZ, &tstep4);
462: tstep = PetscMin(tstep1,tstep2);
463: qp->psteplength = PetscMin(0.95*tstep,1.0);
465: tstep = PetscMin(tstep3,tstep4);
466: qp->dsteplength = PetscMin(0.95*tstep,1.0);
468: qp->psteplength = PetscMin(qp->psteplength,qp->dsteplength);
469: qp->dsteplength = qp->psteplength;
471: return(0);
472: }
477: PetscErrorCode TaoComputeDual_BQPIP(Tao tao,Vec DXL, Vec DXU)
478: {
479: TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
480: PetscErrorCode ierr;
483: VecCopy(qp->Z, DXL);
484: VecCopy(qp->S, DXU);
485: VecScale(DXU, -1.0);
486: return(0);
487: }
492: PetscErrorCode QPIPComputeNormFromCentralPath(TAO_BQPIP *qp, PetscReal *norm)
493: {
494: PetscErrorCode ierr;
495: PetscReal gap[2],mu[2], nmu;
498: VecPointwiseMult(qp->GZwork, qp->G, qp->Z);
499: VecPointwiseMult(qp->TSwork, qp->T, qp->S);
500: VecNorm(qp->TSwork, NORM_1, &mu[0]);
501: VecNorm(qp->GZwork, NORM_1, &mu[1]);
503: nmu=-(mu[0]+mu[1])/qp->m;
505: VecShift(qp->GZwork,nmu);
506: VecShift(qp->TSwork,nmu);
508: VecNorm(qp->GZwork,NORM_2,&gap[0]);
509: VecNorm(qp->TSwork,NORM_2,&gap[1]);
510: gap[0]*=gap[0];
511: gap[1]*=gap[1];
514: qp->pathnorm=PetscSqrtScalar( (gap[0]+gap[1]) );
515: *norm=qp->pathnorm;
517: return(0);
518: }
522: static PetscErrorCode TaoSetFromOptions_BQPIP(PetscOptions *PetscOptionsObject,Tao tao)
523: {
524: TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
528: PetscOptionsHead(PetscOptionsObject,"Interior point method for bound constrained quadratic optimization");
529: PetscOptionsInt("-tao_bqpip_predcorr","Use a predictor-corrector method","",qp->predcorr,&qp->predcorr,0);
530: PetscOptionsTail();
531: KSPSetFromOptions(tao->ksp);
532: return(0);
533: }
537: static PetscErrorCode TaoView_BQPIP(Tao tao, PetscViewer viewer)
538: {
540: return(0);
541: }
543: /* --------------------------------------------------------- */
544: /*MC
545: TAOBQPIP - bounded quadratic interior point algorithm for quadratic
546: optimization with box constraints.
548: Notes: This algorithms solves quadratic problems only, the linear Hessian will
549: only be computed once.
551: Options Database Keys:
552: . -tao_bqpip_predcorr - use a predictor/corrector method
554: Level: beginner
555: M*/
559: PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao)
560: {
561: TAO_BQPIP *qp;
565: PetscNewLog(tao,&qp);
566: tao->ops->setup = TaoSetUp_BQPIP;
567: tao->ops->solve = TaoSolve_BQPIP;
568: tao->ops->view = TaoView_BQPIP;
569: tao->ops->setfromoptions = TaoSetFromOptions_BQPIP;
570: tao->ops->destroy = TaoDestroy_BQPIP;
571: tao->ops->computedual = TaoComputeDual_BQPIP;
573: /* Override default settings (unless already changed) */
574: if (!tao->max_it_changed) tao->max_it=100;
575: if (!tao->max_funcs_changed) tao->max_funcs = 500;
576: #if defined(PETSC_USE_REAL_SINGLE)
577: if (!tao->fatol_changed) tao->fatol=1e-6;
578: if (!tao->frtol_changed) tao->frtol=1e-6;
579: if (!tao->catol_changed) tao->catol=1e-6;
580: #else
581: if (!tao->fatol_changed) tao->fatol=1e-12;
582: if (!tao->frtol_changed) tao->frtol=1e-12;
583: if (!tao->catol_changed) tao->catol=1e-12;
584: #endif
586: /* Initialize pointers and variables */
587: qp->n = 0;
588: qp->m = 0;
589: qp->ksp_tol = 0.1;
591: qp->predcorr = 1;
592: tao->data = (void*)qp;
594: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
595: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
596: KSPSetType(tao->ksp, KSPCG);
597: KSPSetTolerances(tao->ksp, 1e-14, 1e-30, 1e30, PetscMax(10,qp->n));
598: return(0);
599: }