xref: /petsc/src/tao/quadratic/impls/bqpip/bqpip.c (revision bd89dbf26d8a5efecb980364933175da61864cd7) 
1 /*
2     This file implements a Mehrotra predictor-corrector method for
3     bound-constrained quadratic programs.
4 
5  */
6 
7 #include <../src/tao/quadratic/impls/bqpip/bqpipimpl.h>
8 #include <petscksp.h>
9 
QPIPComputeResidual(TAO_BQPIP * qp,Tao tao)10 static PetscErrorCode QPIPComputeResidual(TAO_BQPIP *qp, Tao tao)
11 {
12   PetscReal dtmp = 1.0 - qp->psteplength;
13 
14   PetscFunctionBegin;
15   /* Compute R3 and R5 */
16 
17   PetscCall(VecScale(qp->R3, dtmp));
18   PetscCall(VecScale(qp->R5, dtmp));
19   qp->pinfeas = dtmp * qp->pinfeas;
20 
21   PetscCall(VecCopy(qp->S, tao->gradient));
22   PetscCall(VecAXPY(tao->gradient, -1.0, qp->Z));
23 
24   PetscCall(MatMult(tao->hessian, tao->solution, qp->RHS));
25   PetscCall(VecScale(qp->RHS, -1.0));
26   PetscCall(VecAXPY(qp->RHS, -1.0, qp->C));
27   PetscCall(VecAXPY(tao->gradient, -1.0, qp->RHS));
28 
29   PetscCall(VecNorm(tao->gradient, NORM_1, &qp->dinfeas));
30   qp->rnorm = (qp->dinfeas + qp->pinfeas) / (qp->m + qp->n);
31   PetscFunctionReturn(PETSC_SUCCESS);
32 }
33 
QPIPSetInitialPoint(TAO_BQPIP * qp,Tao tao)34 static PetscErrorCode QPIPSetInitialPoint(TAO_BQPIP *qp, Tao tao)
35 {
36   PetscReal two = 2.0, p01 = 1;
37   PetscReal gap1, gap2, fff, mu;
38 
39   PetscFunctionBegin;
40   /* Compute function, Gradient R=Hx+b, and Hessian */
41   PetscCall(MatMult(tao->hessian, tao->solution, tao->gradient));
42   PetscCall(VecCopy(qp->C, qp->Work));
43   PetscCall(VecAXPY(qp->Work, 0.5, tao->gradient));
44   PetscCall(VecAXPY(tao->gradient, 1.0, qp->C));
45   PetscCall(VecDot(tao->solution, qp->Work, &fff));
46   qp->pobj = fff + qp->d;
47 
48   PetscCheck(!PetscIsInfOrNanReal(qp->pobj), PETSC_COMM_SELF, PETSC_ERR_USER, "User provided data contains infinity or NaN");
49 
50   /* Initialize slack vectors */
51   /* T = XU - X; G = X - XL */
52   PetscCall(VecCopy(qp->XU, qp->T));
53   PetscCall(VecAXPY(qp->T, -1.0, tao->solution));
54   PetscCall(VecCopy(tao->solution, qp->G));
55   PetscCall(VecAXPY(qp->G, -1.0, qp->XL));
56 
57   PetscCall(VecSet(qp->GZwork, p01));
58   PetscCall(VecSet(qp->TSwork, p01));
59 
60   PetscCall(VecPointwiseMax(qp->G, qp->G, qp->GZwork));
61   PetscCall(VecPointwiseMax(qp->T, qp->T, qp->TSwork));
62 
63   /* Initialize Dual Variable Vectors */
64   PetscCall(VecCopy(qp->G, qp->Z));
65   PetscCall(VecReciprocal(qp->Z));
66 
67   PetscCall(VecCopy(qp->T, qp->S));
68   PetscCall(VecReciprocal(qp->S));
69 
70   PetscCall(MatMult(tao->hessian, qp->Work, qp->RHS));
71   PetscCall(VecAbs(qp->RHS));
72   PetscCall(VecSet(qp->Work, p01));
73   PetscCall(VecPointwiseMax(qp->RHS, qp->RHS, qp->Work));
74 
75   PetscCall(VecPointwiseDivide(qp->RHS, tao->gradient, qp->RHS));
76   PetscCall(VecNorm(qp->RHS, NORM_1, &gap1));
77   mu = PetscMin(10.0, (gap1 + 10.0) / qp->m);
78 
79   PetscCall(VecScale(qp->S, mu));
80   PetscCall(VecScale(qp->Z, mu));
81 
82   PetscCall(VecSet(qp->TSwork, p01));
83   PetscCall(VecSet(qp->GZwork, p01));
84   PetscCall(VecPointwiseMax(qp->S, qp->S, qp->TSwork));
85   PetscCall(VecPointwiseMax(qp->Z, qp->Z, qp->GZwork));
86 
87   qp->mu      = 0;
88   qp->dinfeas = 1.0;
89   qp->pinfeas = 1.0;
90   while ((qp->dinfeas + qp->pinfeas) / (qp->m + qp->n) >= qp->mu) {
91     PetscCall(VecScale(qp->G, two));
92     PetscCall(VecScale(qp->Z, two));
93     PetscCall(VecScale(qp->S, two));
94     PetscCall(VecScale(qp->T, two));
95 
96     PetscCall(QPIPComputeResidual(qp, tao));
97 
98     PetscCall(VecCopy(tao->solution, qp->R3));
99     PetscCall(VecAXPY(qp->R3, -1.0, qp->G));
100     PetscCall(VecAXPY(qp->R3, -1.0, qp->XL));
101 
102     PetscCall(VecCopy(tao->solution, qp->R5));
103     PetscCall(VecAXPY(qp->R5, 1.0, qp->T));
104     PetscCall(VecAXPY(qp->R5, -1.0, qp->XU));
105 
106     PetscCall(VecNorm(qp->R3, NORM_INFINITY, &gap1));
107     PetscCall(VecNorm(qp->R5, NORM_INFINITY, &gap2));
108     qp->pinfeas = PetscMax(gap1, gap2);
109 
110     /* Compute the duality gap */
111     PetscCall(VecDot(qp->G, qp->Z, &gap1));
112     PetscCall(VecDot(qp->T, qp->S, &gap2));
113 
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   PetscFunctionReturn(PETSC_SUCCESS);
124 }
125 
QPIPStepLength(TAO_BQPIP * qp)126 static PetscErrorCode QPIPStepLength(TAO_BQPIP *qp)
127 {
128   PetscReal tstep1, tstep2, tstep3, tstep4, tstep;
129 
130   PetscFunctionBegin;
131   /* Compute stepsize to the boundary */
132   PetscCall(VecStepMax(qp->G, qp->DG, &tstep1));
133   PetscCall(VecStepMax(qp->T, qp->DT, &tstep2));
134   PetscCall(VecStepMax(qp->S, qp->DS, &tstep3));
135   PetscCall(VecStepMax(qp->Z, qp->DZ, &tstep4));
136 
137   tstep           = PetscMin(tstep1, tstep2);
138   qp->psteplength = PetscMin(0.95 * tstep, 1.0);
139 
140   tstep           = PetscMin(tstep3, tstep4);
141   qp->dsteplength = PetscMin(0.95 * tstep, 1.0);
142 
143   qp->psteplength = PetscMin(qp->psteplength, qp->dsteplength);
144   qp->dsteplength = qp->psteplength;
145   PetscFunctionReturn(PETSC_SUCCESS);
146 }
147 
QPIPComputeNormFromCentralPath(TAO_BQPIP * qp,PetscReal * norm)148 static PetscErrorCode QPIPComputeNormFromCentralPath(TAO_BQPIP *qp, PetscReal *norm)
149 {
150   PetscReal gap[2], mu[2], nmu;
151 
152   PetscFunctionBegin;
153   PetscCall(VecPointwiseMult(qp->GZwork, qp->G, qp->Z));
154   PetscCall(VecPointwiseMult(qp->TSwork, qp->T, qp->S));
155   PetscCall(VecNorm(qp->TSwork, NORM_1, &mu[0]));
156   PetscCall(VecNorm(qp->GZwork, NORM_1, &mu[1]));
157 
158   nmu = -(mu[0] + mu[1]) / qp->m;
159 
160   PetscCall(VecShift(qp->GZwork, nmu));
161   PetscCall(VecShift(qp->TSwork, nmu));
162 
163   PetscCall(VecNorm(qp->GZwork, NORM_2, &gap[0]));
164   PetscCall(VecNorm(qp->TSwork, NORM_2, &gap[1]));
165   gap[0] *= gap[0];
166   gap[1] *= gap[1];
167 
168   qp->pathnorm = PetscSqrtScalar(gap[0] + gap[1]);
169   *norm        = qp->pathnorm;
170   PetscFunctionReturn(PETSC_SUCCESS);
171 }
172 
QPIPComputeStepDirection(TAO_BQPIP * qp,Tao tao)173 static PetscErrorCode QPIPComputeStepDirection(TAO_BQPIP *qp, Tao tao)
174 {
175   PetscFunctionBegin;
176   /* Calculate DG */
177   PetscCall(VecCopy(tao->stepdirection, qp->DG));
178   PetscCall(VecAXPY(qp->DG, 1.0, qp->R3));
179 
180   /* Calculate DT */
181   PetscCall(VecCopy(tao->stepdirection, qp->DT));
182   PetscCall(VecScale(qp->DT, -1.0));
183   PetscCall(VecAXPY(qp->DT, -1.0, qp->R5));
184 
185   /* Calculate DZ */
186   PetscCall(VecAXPY(qp->DZ, -1.0, qp->Z));
187   PetscCall(VecPointwiseDivide(qp->GZwork, qp->DG, qp->G));
188   PetscCall(VecPointwiseMult(qp->GZwork, qp->GZwork, qp->Z));
189   PetscCall(VecAXPY(qp->DZ, -1.0, qp->GZwork));
190 
191   /* Calculate DS */
192   PetscCall(VecAXPY(qp->DS, -1.0, qp->S));
193   PetscCall(VecPointwiseDivide(qp->TSwork, qp->DT, qp->T));
194   PetscCall(VecPointwiseMult(qp->TSwork, qp->TSwork, qp->S));
195   PetscCall(VecAXPY(qp->DS, -1.0, qp->TSwork));
196   PetscFunctionReturn(PETSC_SUCCESS);
197 }
198 
TaoSetup_BQPIP(Tao tao)199 static PetscErrorCode TaoSetup_BQPIP(Tao tao)
200 {
201   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
202 
203   PetscFunctionBegin;
204   /* Set pointers to Data */
205   PetscCall(VecGetSize(tao->solution, &qp->n));
206 
207   /* Allocate some arrays */
208   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
209   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
210 
211   PetscCall(VecDuplicate(tao->solution, &qp->Work));
212   PetscCall(VecDuplicate(tao->solution, &qp->XU));
213   PetscCall(VecDuplicate(tao->solution, &qp->XL));
214   PetscCall(VecDuplicate(tao->solution, &qp->HDiag));
215   PetscCall(VecDuplicate(tao->solution, &qp->DiagAxpy));
216   PetscCall(VecDuplicate(tao->solution, &qp->RHS));
217   PetscCall(VecDuplicate(tao->solution, &qp->RHS2));
218   PetscCall(VecDuplicate(tao->solution, &qp->C));
219 
220   PetscCall(VecDuplicate(tao->solution, &qp->G));
221   PetscCall(VecDuplicate(tao->solution, &qp->DG));
222   PetscCall(VecDuplicate(tao->solution, &qp->S));
223   PetscCall(VecDuplicate(tao->solution, &qp->Z));
224   PetscCall(VecDuplicate(tao->solution, &qp->DZ));
225   PetscCall(VecDuplicate(tao->solution, &qp->GZwork));
226   PetscCall(VecDuplicate(tao->solution, &qp->R3));
227 
228   PetscCall(VecDuplicate(tao->solution, &qp->T));
229   PetscCall(VecDuplicate(tao->solution, &qp->DT));
230   PetscCall(VecDuplicate(tao->solution, &qp->DS));
231   PetscCall(VecDuplicate(tao->solution, &qp->TSwork));
232   PetscCall(VecDuplicate(tao->solution, &qp->R5));
233   qp->m = 2 * qp->n;
234   PetscFunctionReturn(PETSC_SUCCESS);
235 }
236 
TaoSolve_BQPIP(Tao tao)237 static PetscErrorCode TaoSolve_BQPIP(Tao tao)
238 {
239   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
240   PetscInt   its;
241   PetscReal  d1, d2, ksptol, sigmamu;
242   PetscReal  gnorm, dstep, pstep, step = 0;
243   PetscReal  gap[4];
244   PetscBool  getdiagop;
245 
246   PetscFunctionBegin;
247   qp->dobj        = 0.0;
248   qp->pobj        = 1.0;
249   qp->gap         = 10.0;
250   qp->rgap        = 1.0;
251   qp->mu          = 1.0;
252   qp->dinfeas     = 1.0;
253   qp->psteplength = 0.0;
254   qp->dsteplength = 0.0;
255 
256   /* TODO
257      - Remove fixed variables and treat them correctly
258      - Use index sets for the infinite versus finite bounds
259      - Update remaining code for fixed and free variables
260      - Fix inexact solves for predictor and corrector
261   */
262 
263   /* Tighten infinite bounds, things break when we don't do this
264     -- see test_bqpip.c
265   */
266   PetscCall(TaoComputeVariableBounds(tao));
267   PetscCall(VecSet(qp->XU, 1.0e20));
268   PetscCall(VecSet(qp->XL, -1.0e20));
269   if (tao->XL) PetscCall(VecPointwiseMax(qp->XL, qp->XL, tao->XL));
270   if (tao->XU) PetscCall(VecPointwiseMin(qp->XU, qp->XU, tao->XU));
271   PetscCall(VecMedian(qp->XL, tao->solution, qp->XU, tao->solution));
272 
273   /* Evaluate gradient and Hessian at zero to get the correct values
274      without contaminating them with numerical artifacts.
275   */
276   PetscCall(VecSet(qp->Work, 0));
277   PetscCall(TaoComputeObjectiveAndGradient(tao, qp->Work, &qp->d, qp->C));
278   PetscCall(TaoComputeHessian(tao, qp->Work, tao->hessian, tao->hessian_pre));
279   PetscCall(MatHasOperation(tao->hessian, MATOP_GET_DIAGONAL, &getdiagop));
280   if (getdiagop) PetscCall(MatGetDiagonal(tao->hessian, qp->HDiag));
281 
282   /* Initialize starting point and residuals */
283   PetscCall(QPIPSetInitialPoint(qp, tao));
284   PetscCall(QPIPComputeResidual(qp, tao));
285 
286   /* Enter main loop */
287   tao->reason = TAO_CONTINUE_ITERATING;
288   while (1) {
289     /* Check Stopping Condition      */
290     gnorm = PetscSqrtScalar(qp->gap + qp->dinfeas);
291     PetscCall(TaoLogConvergenceHistory(tao, qp->pobj, gnorm, qp->pinfeas, tao->ksp_its));
292     PetscCall(TaoMonitor(tao, tao->niter, qp->pobj, gnorm, qp->pinfeas, step));
293     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
294     if (tao->reason != TAO_CONTINUE_ITERATING) break;
295     /* Call general purpose update function */
296     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
297     tao->niter++;
298     tao->ksp_its = 0;
299 
300     /*
301        Dual Infeasibility Direction should already be in the right
302        hand side from computing the residuals
303     */
304 
305     PetscCall(QPIPComputeNormFromCentralPath(qp, &d1));
306 
307     PetscCall(VecSet(qp->DZ, 0.0));
308     PetscCall(VecSet(qp->DS, 0.0));
309 
310     /*
311        Compute the Primal Infeasiblitiy RHS and the
312        Diagonal Matrix to be added to H and store in Work
313     */
314     PetscCall(VecPointwiseDivide(qp->DiagAxpy, qp->Z, qp->G));
315     PetscCall(VecPointwiseMult(qp->GZwork, qp->DiagAxpy, qp->R3));
316     PetscCall(VecAXPY(qp->RHS, -1.0, qp->GZwork));
317 
318     PetscCall(VecPointwiseDivide(qp->TSwork, qp->S, qp->T));
319     PetscCall(VecAXPY(qp->DiagAxpy, 1.0, qp->TSwork));
320     PetscCall(VecPointwiseMult(qp->TSwork, qp->TSwork, qp->R5));
321     PetscCall(VecAXPY(qp->RHS, -1.0, qp->TSwork));
322 
323     /*  Determine the solving tolerance */
324     ksptol = qp->mu / 10.0;
325     ksptol = PetscMin(ksptol, 0.001);
326     PetscCall(KSPSetTolerances(tao->ksp, ksptol, 1e-30, 1e30, PetscMax(10, qp->n)));
327 
328     /* Shift the diagonals of the Hessian matrix */
329     PetscCall(MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES));
330     if (!getdiagop) {
331       PetscCall(VecCopy(qp->DiagAxpy, qp->HDiag));
332       PetscCall(VecScale(qp->HDiag, -1.0));
333     }
334     PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
335     PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));
336 
337     PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));
338     PetscCall(KSPSolve(tao->ksp, qp->RHS, tao->stepdirection));
339     PetscCall(KSPGetIterationNumber(tao->ksp, &its));
340     tao->ksp_its += its;
341     tao->ksp_tot_its += its;
342 
343     /* Restore the true diagonal of the Hessian matrix */
344     if (getdiagop) {
345       PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES));
346     } else {
347       PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES));
348     }
349     PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
350     PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));
351     PetscCall(QPIPComputeStepDirection(qp, tao));
352     PetscCall(QPIPStepLength(qp));
353 
354     /* Calculate New Residual R1 in Work vector */
355     PetscCall(MatMult(tao->hessian, tao->stepdirection, qp->RHS2));
356     PetscCall(VecAXPY(qp->RHS2, 1.0, qp->DS));
357     PetscCall(VecAXPY(qp->RHS2, -1.0, qp->DZ));
358     PetscCall(VecAYPX(qp->RHS2, qp->dsteplength, tao->gradient));
359 
360     PetscCall(VecNorm(qp->RHS2, NORM_2, &qp->dinfeas));
361     PetscCall(VecDot(qp->DZ, qp->DG, gap));
362     PetscCall(VecDot(qp->DS, qp->DT, gap + 1));
363 
364     qp->rnorm = (qp->dinfeas + qp->psteplength * qp->pinfeas) / (qp->m + qp->n);
365     pstep     = qp->psteplength;
366     step      = PetscMin(qp->psteplength, qp->dsteplength);
367     sigmamu   = (pstep * pstep * (gap[0] + gap[1]) + (1 - pstep) * qp->gap) / qp->m;
368 
369     if (qp->predcorr && step < 0.9) {
370       if (sigmamu < qp->mu) {
371         sigmamu = sigmamu / qp->mu;
372         sigmamu = sigmamu * sigmamu * sigmamu;
373       } else {
374         sigmamu = 1.0;
375       }
376       sigmamu = sigmamu * qp->mu;
377 
378       /* Compute Corrector Step */
379       PetscCall(VecPointwiseMult(qp->DZ, qp->DG, qp->DZ));
380       PetscCall(VecScale(qp->DZ, -1.0));
381       PetscCall(VecShift(qp->DZ, sigmamu));
382       PetscCall(VecPointwiseDivide(qp->DZ, qp->DZ, qp->G));
383 
384       PetscCall(VecPointwiseMult(qp->DS, qp->DS, qp->DT));
385       PetscCall(VecScale(qp->DS, -1.0));
386       PetscCall(VecShift(qp->DS, sigmamu));
387       PetscCall(VecPointwiseDivide(qp->DS, qp->DS, qp->T));
388 
389       PetscCall(VecCopy(qp->DZ, qp->RHS2));
390       PetscCall(VecAXPY(qp->RHS2, -1.0, qp->DS));
391       PetscCall(VecAXPY(qp->RHS2, 1.0, qp->RHS));
392 
393       /* Approximately solve the linear system */
394       PetscCall(MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES));
395       if (!getdiagop) {
396         PetscCall(VecCopy(qp->DiagAxpy, qp->HDiag));
397         PetscCall(VecScale(qp->HDiag, -1.0));
398       }
399       PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
400       PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));
401 
402       /* Solve using the previous tolerances that were set */
403       PetscCall(KSPSolve(tao->ksp, qp->RHS2, tao->stepdirection));
404       PetscCall(KSPGetIterationNumber(tao->ksp, &its));
405       tao->ksp_its += its;
406       tao->ksp_tot_its += its;
407 
408       if (getdiagop) {
409         PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES));
410       } else {
411         PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES));
412       }
413       PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
414       PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));
415       PetscCall(QPIPComputeStepDirection(qp, tao));
416       PetscCall(QPIPStepLength(qp));
417     } /* End Corrector step */
418 
419     /* Take the step */
420     dstep = qp->dsteplength;
421 
422     PetscCall(VecAXPY(qp->Z, dstep, qp->DZ));
423     PetscCall(VecAXPY(qp->S, dstep, qp->DS));
424     PetscCall(VecAXPY(tao->solution, dstep, tao->stepdirection));
425     PetscCall(VecAXPY(qp->G, dstep, qp->DG));
426     PetscCall(VecAXPY(qp->T, dstep, qp->DT));
427 
428     /* Compute Residuals */
429     PetscCall(QPIPComputeResidual(qp, tao));
430 
431     /* Evaluate quadratic function */
432     PetscCall(MatMult(tao->hessian, tao->solution, qp->Work));
433 
434     PetscCall(VecDot(tao->solution, qp->Work, &d1));
435     PetscCall(VecDot(tao->solution, qp->C, &d2));
436     PetscCall(VecDot(qp->G, qp->Z, gap));
437     PetscCall(VecDot(qp->T, qp->S, gap + 1));
438 
439     /* Compute the duality gap */
440     qp->pobj = d1 / 2.0 + d2 + qp->d;
441     qp->gap  = gap[0] + gap[1];
442     qp->dobj = qp->pobj - qp->gap;
443     if (qp->m > 0) qp->mu = qp->gap / (qp->m);
444     qp->rgap = qp->gap / (PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
445   } /* END MAIN LOOP  */
446   PetscFunctionReturn(PETSC_SUCCESS);
447 }
448 
TaoView_BQPIP(Tao tao,PetscViewer viewer)449 static PetscErrorCode TaoView_BQPIP(Tao tao, PetscViewer viewer)
450 {
451   PetscFunctionBegin;
452   PetscFunctionReturn(PETSC_SUCCESS);
453 }
454 
TaoSetFromOptions_BQPIP(Tao tao,PetscOptionItems PetscOptionsObject)455 static PetscErrorCode TaoSetFromOptions_BQPIP(Tao tao, PetscOptionItems PetscOptionsObject)
456 {
457   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
458 
459   PetscFunctionBegin;
460   PetscOptionsHeadBegin(PetscOptionsObject, "Interior point method for bound constrained quadratic optimization");
461   PetscCall(PetscOptionsInt("-tao_bqpip_predcorr", "Use a predictor-corrector method", "", qp->predcorr, &qp->predcorr, NULL));
462   PetscOptionsHeadEnd();
463   PetscCall(KSPSetFromOptions(tao->ksp));
464   PetscFunctionReturn(PETSC_SUCCESS);
465 }
466 
TaoDestroy_BQPIP(Tao tao)467 static PetscErrorCode TaoDestroy_BQPIP(Tao tao)
468 {
469   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
470 
471   PetscFunctionBegin;
472   if (tao->setupcalled) {
473     PetscCall(VecDestroy(&qp->G));
474     PetscCall(VecDestroy(&qp->DG));
475     PetscCall(VecDestroy(&qp->Z));
476     PetscCall(VecDestroy(&qp->DZ));
477     PetscCall(VecDestroy(&qp->GZwork));
478     PetscCall(VecDestroy(&qp->R3));
479     PetscCall(VecDestroy(&qp->S));
480     PetscCall(VecDestroy(&qp->DS));
481     PetscCall(VecDestroy(&qp->T));
482 
483     PetscCall(VecDestroy(&qp->DT));
484     PetscCall(VecDestroy(&qp->TSwork));
485     PetscCall(VecDestroy(&qp->R5));
486     PetscCall(VecDestroy(&qp->HDiag));
487     PetscCall(VecDestroy(&qp->Work));
488     PetscCall(VecDestroy(&qp->XL));
489     PetscCall(VecDestroy(&qp->XU));
490     PetscCall(VecDestroy(&qp->DiagAxpy));
491     PetscCall(VecDestroy(&qp->RHS));
492     PetscCall(VecDestroy(&qp->RHS2));
493     PetscCall(VecDestroy(&qp->C));
494   }
495   PetscCall(KSPDestroy(&tao->ksp));
496   PetscCall(PetscFree(tao->data));
497   PetscFunctionReturn(PETSC_SUCCESS);
498 }
499 
TaoComputeDual_BQPIP(Tao tao,Vec DXL,Vec DXU)500 static PetscErrorCode TaoComputeDual_BQPIP(Tao tao, Vec DXL, Vec DXU)
501 {
502   TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
503 
504   PetscFunctionBegin;
505   PetscCall(VecCopy(qp->Z, DXL));
506   PetscCall(VecCopy(qp->S, DXU));
507   PetscCall(VecScale(DXU, -1.0));
508   PetscFunctionReturn(PETSC_SUCCESS);
509 }
510 
511 /*MC
512  TAOBQPIP - interior-point method for quadratic programs with
513     box constraints.
514 
515  Notes:
516     This algorithms solves quadratic problems only, the Hessian will
517         only be computed once.
518 
519  Options Database Keys:
520 . -tao_bqpip_predcorr - use a predictor/corrector method
521 
522   Level: beginner
523 M*/
524 
TaoCreate_BQPIP(Tao tao)525 PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao)
526 {
527   TAO_BQPIP *qp;
528 
529   PetscFunctionBegin;
530   PetscCall(PetscNew(&qp));
531 
532   tao->ops->setup          = TaoSetup_BQPIP;
533   tao->ops->solve          = TaoSolve_BQPIP;
534   tao->ops->view           = TaoView_BQPIP;
535   tao->ops->setfromoptions = TaoSetFromOptions_BQPIP;
536   tao->ops->destroy        = TaoDestroy_BQPIP;
537   tao->ops->computedual    = TaoComputeDual_BQPIP;
538 
539   /* Override default settings (unless already changed) */
540   PetscCall(TaoParametersInitialize(tao));
541   PetscObjectParameterSetDefault(tao, max_it, 100);
542   PetscObjectParameterSetDefault(tao, max_funcs, 500);
543   PetscObjectParameterSetDefault(tao, catol, PetscDefined(USE_REAL_SINGLE) ? 1e-6 : 1e-12);
544 
545   /* Initialize pointers and variables */
546   qp->n = 0;
547   qp->m = 0;
548 
549   qp->predcorr = 1;
550   tao->data    = (void *)qp;
551 
552   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
553   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
554   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
555   PetscCall(KSPSetType(tao->ksp, KSPCG));
556   PetscCall(KSPSetTolerances(tao->ksp, 1e-14, 1e-30, 1e30, PetscMax(10, qp->n)));
557   PetscFunctionReturn(PETSC_SUCCESS);
558 }
559