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