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: {
12: PetscReal dtmp = 1.0 - qp->psteplength;
14: PetscFunctionBegin;
15: /* Compute R3 and R5 */
17: PetscCall(VecScale(qp->R3, dtmp));
18: PetscCall(VecScale(qp->R5, dtmp));
19: qp->pinfeas = dtmp * qp->pinfeas;
21: PetscCall(VecCopy(qp->S, tao->gradient));
22: PetscCall(VecAXPY(tao->gradient, -1.0, qp->Z));
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));
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: }
34: static PetscErrorCode QPIPSetInitialPoint(TAO_BQPIP *qp, Tao tao)
35: {
36: PetscReal two = 2.0, p01 = 1;
37: PetscReal gap1, gap2, fff, mu;
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;
48: PetscCheck(!PetscIsInfOrNanReal(qp->pobj), PETSC_COMM_SELF, PETSC_ERR_USER, "User provided data contains infinity or NaN");
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));
57: PetscCall(VecSet(qp->GZwork, p01));
58: PetscCall(VecSet(qp->TSwork, p01));
60: PetscCall(VecPointwiseMax(qp->G, qp->G, qp->GZwork));
61: PetscCall(VecPointwiseMax(qp->T, qp->T, qp->TSwork));
63: /* Initialize Dual Variable Vectors */
64: PetscCall(VecCopy(qp->G, qp->Z));
65: PetscCall(VecReciprocal(qp->Z));
67: PetscCall(VecCopy(qp->T, qp->S));
68: PetscCall(VecReciprocal(qp->S));
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));
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);
79: PetscCall(VecScale(qp->S, mu));
80: PetscCall(VecScale(qp->Z, mu));
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));
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));
96: PetscCall(QPIPComputeResidual(qp, tao));
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));
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));
106: PetscCall(VecNorm(qp->R3, NORM_INFINITY, &gap1));
107: PetscCall(VecNorm(qp->R5, NORM_INFINITY, &gap2));
108: qp->pinfeas = PetscMax(gap1, gap2);
110: /* Compute the duality gap */
111: PetscCall(VecDot(qp->G, qp->Z, &gap1));
112: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode QPIPStepLength(TAO_BQPIP *qp)
127: {
128: PetscReal tstep1, tstep2, tstep3, tstep4, tstep;
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));
137: tstep = PetscMin(tstep1, tstep2);
138: qp->psteplength = PetscMin(0.95 * tstep, 1.0);
140: tstep = PetscMin(tstep3, tstep4);
141: qp->dsteplength = PetscMin(0.95 * tstep, 1.0);
143: qp->psteplength = PetscMin(qp->psteplength, qp->dsteplength);
144: qp->dsteplength = qp->psteplength;
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: static PetscErrorCode QPIPComputeNormFromCentralPath(TAO_BQPIP *qp, PetscReal *norm)
149: {
150: PetscReal gap[2], mu[2], nmu;
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]));
158: nmu = -(mu[0] + mu[1]) / qp->m;
160: PetscCall(VecShift(qp->GZwork, nmu));
161: PetscCall(VecShift(qp->TSwork, nmu));
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];
168: qp->pathnorm = PetscSqrtScalar(gap[0] + gap[1]);
169: *norm = qp->pathnorm;
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
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));
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));
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));
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: }
199: static PetscErrorCode TaoSetup_BQPIP(Tao tao)
200: {
201: TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
203: PetscFunctionBegin;
204: /* Set pointers to Data */
205: PetscCall(VecGetSize(tao->solution, &qp->n));
207: /* Allocate some arrays */
208: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
209: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
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));
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));
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: }
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;
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;
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: */
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));
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));
282: /* Initialize starting point and residuals */
283: PetscCall(QPIPSetInitialPoint(qp, tao));
284: PetscCall(QPIPComputeResidual(qp, tao));
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;
300: /*
301: Dual Infeasibility Direction should already be in the right
302: hand side from computing the residuals
303: */
305: PetscCall(QPIPComputeNormFromCentralPath(qp, &d1));
307: PetscCall(VecSet(qp->DZ, 0.0));
308: PetscCall(VecSet(qp->DS, 0.0));
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));
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));
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)));
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));
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;
343: /* Restore the true diagonal of the Hessian matrix */
344: if (getdiagop) PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES));
345: else PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES));
346: PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
347: PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));
348: PetscCall(QPIPComputeStepDirection(qp, tao));
349: PetscCall(QPIPStepLength(qp));
351: /* Calculate New Residual R1 in Work vector */
352: PetscCall(MatMult(tao->hessian, tao->stepdirection, qp->RHS2));
353: PetscCall(VecAXPY(qp->RHS2, 1.0, qp->DS));
354: PetscCall(VecAXPY(qp->RHS2, -1.0, qp->DZ));
355: PetscCall(VecAYPX(qp->RHS2, qp->dsteplength, tao->gradient));
357: PetscCall(VecNorm(qp->RHS2, NORM_2, &qp->dinfeas));
358: PetscCall(VecDot(qp->DZ, qp->DG, gap));
359: PetscCall(VecDot(qp->DS, qp->DT, gap + 1));
361: qp->rnorm = (qp->dinfeas + qp->psteplength * qp->pinfeas) / (qp->m + qp->n);
362: pstep = qp->psteplength;
363: step = PetscMin(qp->psteplength, qp->dsteplength);
364: sigmamu = (pstep * pstep * (gap[0] + gap[1]) + (1 - pstep) * qp->gap) / qp->m;
366: if (qp->predcorr && step < 0.9) {
367: if (sigmamu < qp->mu) {
368: sigmamu = sigmamu / qp->mu;
369: sigmamu = sigmamu * sigmamu * sigmamu;
370: } else {
371: sigmamu = 1.0;
372: }
373: sigmamu = sigmamu * qp->mu;
375: /* Compute Corrector Step */
376: PetscCall(VecPointwiseMult(qp->DZ, qp->DG, qp->DZ));
377: PetscCall(VecScale(qp->DZ, -1.0));
378: PetscCall(VecShift(qp->DZ, sigmamu));
379: PetscCall(VecPointwiseDivide(qp->DZ, qp->DZ, qp->G));
381: PetscCall(VecPointwiseMult(qp->DS, qp->DS, qp->DT));
382: PetscCall(VecScale(qp->DS, -1.0));
383: PetscCall(VecShift(qp->DS, sigmamu));
384: PetscCall(VecPointwiseDivide(qp->DS, qp->DS, qp->T));
386: PetscCall(VecCopy(qp->DZ, qp->RHS2));
387: PetscCall(VecAXPY(qp->RHS2, -1.0, qp->DS));
388: PetscCall(VecAXPY(qp->RHS2, 1.0, qp->RHS));
390: /* Approximately solve the linear system */
391: PetscCall(MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES));
392: if (!getdiagop) {
393: PetscCall(VecCopy(qp->DiagAxpy, qp->HDiag));
394: PetscCall(VecScale(qp->HDiag, -1.0));
395: }
396: PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
397: PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));
399: /* Solve using the previous tolerances that were set */
400: PetscCall(KSPSolve(tao->ksp, qp->RHS2, tao->stepdirection));
401: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
402: tao->ksp_its += its;
403: tao->ksp_tot_its += its;
405: if (getdiagop) PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES));
406: else PetscCall(MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES));
407: PetscCall(MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY));
408: PetscCall(MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY));
409: PetscCall(QPIPComputeStepDirection(qp, tao));
410: PetscCall(QPIPStepLength(qp));
411: } /* End Corrector step */
413: /* Take the step */
414: dstep = qp->dsteplength;
416: PetscCall(VecAXPY(qp->Z, dstep, qp->DZ));
417: PetscCall(VecAXPY(qp->S, dstep, qp->DS));
418: PetscCall(VecAXPY(tao->solution, dstep, tao->stepdirection));
419: PetscCall(VecAXPY(qp->G, dstep, qp->DG));
420: PetscCall(VecAXPY(qp->T, dstep, qp->DT));
422: /* Compute Residuals */
423: PetscCall(QPIPComputeResidual(qp, tao));
425: /* Evaluate quadratic function */
426: PetscCall(MatMult(tao->hessian, tao->solution, qp->Work));
428: PetscCall(VecDot(tao->solution, qp->Work, &d1));
429: PetscCall(VecDot(tao->solution, qp->C, &d2));
430: PetscCall(VecDot(qp->G, qp->Z, gap));
431: PetscCall(VecDot(qp->T, qp->S, gap + 1));
433: /* Compute the duality gap */
434: qp->pobj = d1 / 2.0 + d2 + qp->d;
435: qp->gap = gap[0] + gap[1];
436: qp->dobj = qp->pobj - qp->gap;
437: if (qp->m > 0) qp->mu = qp->gap / (qp->m);
438: qp->rgap = qp->gap / (PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
439: } /* END MAIN LOOP */
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: static PetscErrorCode TaoView_BQPIP(Tao tao, PetscViewer viewer)
444: {
445: PetscFunctionBegin;
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: static PetscErrorCode TaoSetFromOptions_BQPIP(Tao tao, PetscOptionItems PetscOptionsObject)
450: {
451: TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
453: PetscFunctionBegin;
454: PetscOptionsHeadBegin(PetscOptionsObject, "Interior point method for bound constrained quadratic optimization");
455: PetscCall(PetscOptionsInt("-tao_bqpip_predcorr", "Use a predictor-corrector method", "", qp->predcorr, &qp->predcorr, NULL));
456: PetscOptionsHeadEnd();
457: PetscCall(KSPSetFromOptions(tao->ksp));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: static PetscErrorCode TaoDestroy_BQPIP(Tao tao)
462: {
463: TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
465: PetscFunctionBegin;
466: if (tao->setupcalled) {
467: PetscCall(VecDestroy(&qp->G));
468: PetscCall(VecDestroy(&qp->DG));
469: PetscCall(VecDestroy(&qp->Z));
470: PetscCall(VecDestroy(&qp->DZ));
471: PetscCall(VecDestroy(&qp->GZwork));
472: PetscCall(VecDestroy(&qp->R3));
473: PetscCall(VecDestroy(&qp->S));
474: PetscCall(VecDestroy(&qp->DS));
475: PetscCall(VecDestroy(&qp->T));
477: PetscCall(VecDestroy(&qp->DT));
478: PetscCall(VecDestroy(&qp->TSwork));
479: PetscCall(VecDestroy(&qp->R5));
480: PetscCall(VecDestroy(&qp->HDiag));
481: PetscCall(VecDestroy(&qp->Work));
482: PetscCall(VecDestroy(&qp->XL));
483: PetscCall(VecDestroy(&qp->XU));
484: PetscCall(VecDestroy(&qp->DiagAxpy));
485: PetscCall(VecDestroy(&qp->RHS));
486: PetscCall(VecDestroy(&qp->RHS2));
487: PetscCall(VecDestroy(&qp->C));
488: }
489: PetscCall(KSPDestroy(&tao->ksp));
490: PetscCall(PetscFree(tao->data));
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: static PetscErrorCode TaoComputeDual_BQPIP(Tao tao, Vec DXL, Vec DXU)
495: {
496: TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
498: PetscFunctionBegin;
499: PetscCall(VecCopy(qp->Z, DXL));
500: PetscCall(VecCopy(qp->S, DXU));
501: PetscCall(VecScale(DXU, -1.0));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: /*MC
506: TAOBQPIP - interior-point method for quadratic programs with
507: box constraints.
509: Options Database Key:
510: . -tao_bqpip_predcorr - use a predictor/corrector method
512: Level: beginner
514: Note:
515: This algorithm solves quadratic problems only, the Hessian will only be computed once.
517: .seealso: `Tao`, `TaoType`, `TAOGPCG`, `TAOTRON`
518: M*/
520: PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao)
521: {
522: TAO_BQPIP *qp;
524: PetscFunctionBegin;
525: PetscCall(PetscNew(&qp));
527: tao->ops->setup = TaoSetup_BQPIP;
528: tao->ops->solve = TaoSolve_BQPIP;
529: tao->ops->view = TaoView_BQPIP;
530: tao->ops->setfromoptions = TaoSetFromOptions_BQPIP;
531: tao->ops->destroy = TaoDestroy_BQPIP;
532: tao->ops->computedual = TaoComputeDual_BQPIP;
533: tao->uses_gradient = PETSC_TRUE;
534: tao->uses_hessian_matrices = PETSC_TRUE;
536: /* Override default settings (unless already changed) */
537: PetscCall(TaoParametersInitialize(tao));
538: PetscObjectParameterSetDefault(tao, max_it, 100);
539: PetscObjectParameterSetDefault(tao, max_funcs, 500);
540: PetscObjectParameterSetDefault(tao, catol, PetscDefined(USE_REAL_SINGLE) ? 1e-6 : 1e-12);
542: /* Initialize pointers and variables */
543: qp->n = 0;
544: qp->m = 0;
546: qp->predcorr = 1;
547: tao->data = (void *)qp;
549: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
550: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
551: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
552: PetscCall(KSPSetType(tao->ksp, KSPCG));
553: PetscCall(KSPSetTolerances(tao->ksp, 1e-14, 1e-30, 1e30, PetscMax(10, qp->n)));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }