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