Actual source code: bqpip.c

  1: /*$Id$*/

  3: #include "bqpip.h"
  4: #include "taolinearsolver.h"
  5: static int TaoSetDown_BQPIP(TAO_SOLVER, void*);

  9: static int TaoSetUp_BQPIP(TAO_SOLVER tao, void *solver)
 10: {
 11:   TAO_BQPIP *qp =(TAO_BQPIP*)solver;
 12:   TaoInt       n, info;
 13:   TaoLinearSolver *ksp;
 14:   TaoFunctionBegin;

 16:   /* Set pointers to Data */
 17:   info = TaoGetHessian(tao,&qp->H);
 18:   info = TaoGetSolution(tao,&qp->XY);CHKERRQ(info);
 19:   info = qp->XY->GetDimension(&qp->n); CHKERRQ(info);

 21:   /* Allocate some arrays */
 22:   info = qp->XY->Clone(&qp->Work); CHKERRQ(info);
 23:   info = qp->XY->Clone(&qp->DXY); CHKERRQ(info);
 24:   info = qp->XY->Clone(&qp->HDiag); CHKERRQ(info);
 25:   info = qp->XY->Clone(&qp->DiagAxpy); CHKERRQ(info);
 26:   info = qp->XY->Clone(&qp->RHS); CHKERRQ(info);
 27:   info = qp->XY->Clone(&qp->RHS2); CHKERRQ(info);
 28:   info = qp->XY->Clone(&qp->C0); CHKERRQ(info);
 29:   info = qp->XY->Clone(&qp->R12); CHKERRQ(info);

 31:   info = qp->XY->Clone(&qp->XL); CHKERRQ(info);
 32:   info = qp->XL->Clone(&qp->G); CHKERRQ(info);
 33:   info = qp->XL->Clone(&qp->DG); CHKERRQ(info);
 34:   info = qp->XL->Clone(&qp->Z); CHKERRQ(info);
 35:   info = qp->XL->Clone(&qp->DZ); CHKERRQ(info);
 36:   info = qp->XL->Clone(&qp->GZwork); CHKERRQ(info);
 37:   info = qp->XL->Clone(&qp->R3); CHKERRQ(info);

 39:   info = qp->XY->Clone(&qp->XU); CHKERRQ(info);
 40:   info = qp->XU->Clone(&qp->T); CHKERRQ(info);
 41:   info = qp->XU->Clone(&qp->DT); CHKERRQ(info);
 42:   info = qp->XU->Clone(&qp->S); CHKERRQ(info);
 43:   info = qp->XU->Clone(&qp->DS); CHKERRQ(info);
 44:   info = qp->XU->Clone(&qp->TSwork); CHKERRQ(info);
 45:   info = qp->XU->Clone(&qp->R5); CHKERRQ(info);

 47:   info = TaoSetLagrangianGradientVector(tao,qp->R12);CHKERRQ(info);
 48:   info = TaoSetVariableBounds(tao,qp->XL,qp->XU);CHKERRQ(info);
 49:   info = TaoSetStepDirectionVector(tao,qp->DXY);CHKERRQ(info);

 51:   /* Register the Events */
 52:   info = qp->XY->GetDimension(&qp->n); CHKERRQ(info);

 54:   qp->m=0;
 55:   info=qp->G->GetDimension(&n); CHKERRQ(info); qp->m+=n;
 56:   info=qp->T->GetDimension(&n); CHKERRQ(info); qp->m+=n;

 58:   info = TaoCreateLinearSolver(tao,qp->H,300,&ksp); CHKERRQ(info);
 59:   info = ksp->SetTolerances(1e-14,1e-30,1e30,qp->n); CHKERRQ(info);

 61:   TaoFunctionReturn(0);
 62: }

 66: static int  QPIPSetInitialPoint(TAO_SOLVER tao, TAO_BQPIP *qp)
 67: {
 68:   int       info;
 69:   double    two=2.0,p01=1;
 70:   double    gap1,gap2,fff,mu;

 72:   TaoFunctionBegin;
 73:   /* Compute function, Gradient R=Hx+b, and Hessian */
 74:   info = qp->XY->Median(qp->XL,qp->XY,qp->XU); CHKERRQ(info);
 75:   info = qp->H->Multiply(qp->XY,qp->R12); CHKERRQ(info);

 77:   info = qp->Work->Waxpby(0.5,qp->R12,1.0,qp->C0);CHKERRQ(info);
 78:   info = qp->R12->Axpy(1.0,qp->C0);CHKERRQ(info);
 79:   info = qp->Work->Dot(qp->XY,&fff);CHKERRQ(info);
 80:   qp->pobj = fff + qp->c;

 82:   /* Initialize Primal Vectors */

 84:   info = qp->T->Waxpby(1.0,qp->XU,-1.0,qp->XY);CHKERRQ(info);
 85:   info = qp->G->Waxpby(1.0,qp->XY,-1.0,qp->XL);CHKERRQ(info);

 87:   info = qp->GZwork->SetToConstant(p01);CHKERRQ(info);
 88:   info = qp->TSwork->SetToConstant(p01);CHKERRQ(info);

 90:   info = qp->G->PointwiseMaximum(qp->G,qp->GZwork); CHKERRQ(info);
 91:   info = qp->T->PointwiseMaximum(qp->T,qp->TSwork); CHKERRQ(info);

 93:   /* Initialize Dual Variable Vectors */

 95:   info = qp->Z->CopyFrom(qp->G); CHKERRQ(info);
 96:   info = qp->Z->Reciprocal(); CHKERRQ(info);

 98:   info = qp->S->CopyFrom(qp->T); CHKERRQ(info);
 99:   info = qp->S->Reciprocal(); CHKERRQ(info);

101:   info = qp->H->Multiply(qp->Work,qp->RHS); CHKERRQ(info);
102:   info = qp->RHS->AbsoluteValue(); CHKERRQ(info);
103:   info = qp->Work->SetToConstant(p01);CHKERRQ(info);
104:   info = qp->RHS->PointwiseMaximum(qp->RHS,qp->Work); CHKERRQ(info);

106:   info = qp->RHS->PointwiseDivide(qp->R12,qp->RHS); CHKERRQ(info);
107:   info = qp->RHS->Norm1(&gap1); CHKERRQ(info);
108:   mu = TaoMin(10.0,(gap1+10.0)/qp->m);

110:   info = qp->S->Scale(mu); CHKERRQ(info);
111:   info = qp->Z->Scale(mu); CHKERRQ(info);

113:   info = qp->TSwork->SetToConstant(p01); CHKERRQ(info);
114:   info = qp->GZwork->SetToConstant(p01); CHKERRQ(info);
115:   info = qp->S->PointwiseMaximum(qp->S,qp->TSwork); CHKERRQ(info);
116:   info = qp->Z->PointwiseMaximum(qp->Z,qp->GZwork); CHKERRQ(info);

118:   qp->mu=0;qp->dinfeas=1.0;qp->pinfeas=1.0;
119:   while ( (qp->dinfeas+qp->pinfeas)/(qp->m+qp->n) >= qp->mu ){

121:     info=qp->G->Scale(two); CHKERRQ(info);
122:     info=qp->Z->Scale(two); CHKERRQ(info);
123:     info=qp->S->Scale(two); CHKERRQ(info);
124:     info=qp->T->Scale(two); CHKERRQ(info);

126:     info = QPIPComputeResidual(qp); CHKERRQ(info);
127:     
128:     info=qp->R3->Waxpby(1.0,qp->XY,-1.0,qp->G);CHKERRQ(info);
129:     info=qp->R3->Axpy(-1.0,qp->XL);CHKERRQ(info);

131:     info=qp->R5->Waxpby(1.0,qp->XY,1.0,qp->T);CHKERRQ(info);
132:     info=qp->R5->Axpy(-1.0,qp->XU);CHKERRQ(info);
133:     
134:     info=qp->R3->NormInfinity(&gap1);CHKERRQ(info);
135:     info=qp->R5->NormInfinity(&gap2);CHKERRQ(info);
136:     qp->pinfeas=TaoMax(gap1,gap2);
137:     
138:     /* Compute the duality gap */
139:     info=qp->G->Dot(qp->Z,&gap1);CHKERRQ(info);
140:     info=qp->T->Dot(qp->S,&gap2);CHKERRQ(info);
141:     
142:     qp->gap = (gap1+gap2);
143:     qp->dobj = qp->pobj - qp->gap;
144:     if (qp->m>0) qp->mu=qp->gap/(qp->m); else qp->mu=0.0;
145:     qp->rgap=qp->gap/( TaoAbsScalar(qp->dobj) + TaoAbsScalar(qp->pobj) + 1.0 );
146:   }

148:   TaoFunctionReturn(0);
149: }


154: static int TaoSetDown_BQPIP(TAO_SOLVER tao, void*solver)
155: {
156:   TAO_BQPIP *qp = (TAO_BQPIP*)solver;
157:   int info;

159:   /* Free allocated memory in GPCG structure */
160:   TaoFunctionBegin;

162:   info=TaoVecDestroy(qp->G);CHKERRQ(info);
163:   info=TaoVecDestroy(qp->DG);CHKERRQ(info);
164:   info=TaoVecDestroy(qp->Z);CHKERRQ(info);
165:   info=TaoVecDestroy(qp->DZ);CHKERRQ(info);
166:   info=TaoVecDestroy(qp->GZwork);CHKERRQ(info);
167:   info=TaoVecDestroy(qp->R3);CHKERRQ(info);
168:   info=TaoVecDestroy(qp->XL);CHKERRQ(info);
169:   
170:   info=TaoVecDestroy(qp->S);CHKERRQ(info);
171:   info=TaoVecDestroy(qp->DS);CHKERRQ(info);
172:   info=TaoVecDestroy(qp->T);CHKERRQ(info);
173:   info=TaoVecDestroy(qp->DT);CHKERRQ(info);
174:   info=TaoVecDestroy(qp->TSwork);CHKERRQ(info);
175:   info=TaoVecDestroy(qp->R5);CHKERRQ(info);
176:   info=TaoVecDestroy(qp->XU);CHKERRQ(info);
177:   
178:   info=TaoVecDestroy(qp->HDiag);CHKERRQ(info);
179:   info=TaoVecDestroy(qp->Work);CHKERRQ(info);
180:   info=TaoVecDestroy(qp->DiagAxpy);CHKERRQ(info);
181:   info=TaoVecDestroy(qp->RHS);CHKERRQ(info);
182:   info=TaoVecDestroy(qp->RHS2);CHKERRQ(info);
183:   info=TaoVecDestroy(qp->DXY);CHKERRQ(info);
184:   info=TaoVecDestroy(qp->C0);CHKERRQ(info);
185:   info=TaoVecDestroy(qp->R12);CHKERRQ(info);
186:   
187:   info = TaoDestroyLinearSolver(tao);CHKERRQ(info);

189:   TaoFunctionReturn(0);
190: }

194: static int TaoSolve_BQPIP(TAO_SOLVER tao, void *solver)
195: {
196:   TAO_BQPIP *qp = (TAO_BQPIP*)solver;
197:   int       info;
198:   TaoInt    iter=0;
199:   TaoInt    n;
200:   double    d1,d2,ksptol,sigma;
201:   double    sigmamu;
202:   double    dstep,pstep,step=0;
203:   double    gap[4];
204:   TaoTruth  kspsuccess;
205:   TaoTerminateReason reason;
206:   
207:   TaoFunctionBegin;

209:   info = TaoGetSolution(tao,&qp->XY);CHKERRQ(info);
210:   info = TaoEvaluateVariableBounds(tao,qp->XL,qp->XU); CHKERRQ(info);
211:   info = qp->XY->GetDimension(&n);CHKERRQ(info);

213:   info = TaoComputeFunctionGradient(tao,qp->XY,&qp->c,qp->C0); CHKERRQ(info);
214:   info = TaoComputeHessian(tao,qp->XY,qp->H); CHKERRQ(info);
215:   info = qp->H->Multiply(qp->XY,qp->Work); CHKERRQ(info);
216:   info = qp->XY->Dot(qp->Work,&d1); CHKERRQ(info);
217:   info = qp->C0->Axpy(-1.0,qp->Work); CHKERRQ(info);
218:   info = qp->XY->Dot(qp->C0,&d2); CHKERRQ(info);
219:   qp->c -= (d1/2.0+d2);

221:   info = qp->H->GetDiagonal(qp->HDiag); CHKERRQ(info);

223:   info = QPIPSetInitialPoint(tao,qp); CHKERRQ(info);
224:   info = QPIPComputeResidual(qp); CHKERRQ(info);
225:   
226:   /* Enter main loop */
227:   while (1){

229:     /* Check Stopping Condition      */
230:     info=TaoMonitor(tao,iter++,qp->pobj,sqrt(qp->gap+qp->dinfeas),qp->pinfeas,
231:                     step,&reason); CHKERRQ(info);
232:     if (reason != TAO_CONTINUE_ITERATING) break;

234:     /* 
235:        Dual Infeasibility Direction should already be in the right
236:        hand side from computing the residuals 
237:     */

239:     info = TAOComputeNormFromCentralPath_BQPIP(tao,&d1); CHKERRQ(info);

241:     if (iter > 0 && (qp->rnorm>5*qp->mu || d1*d1>qp->m*qp->mu*qp->mu) ) {
242:       sigma=1.0;sigmamu=qp->mu;
243:       sigma=0.0;sigmamu=0;
244:     } else {
245:       sigma=0.0;sigmamu=0;
246:     }
247:     info = qp->DZ->SetToConstant(sigmamu); CHKERRQ(info);
248:     info = qp->DS->SetToConstant(sigmamu); CHKERRQ(info);

250:     if (sigmamu !=0){
251:       info = qp->DZ->PointwiseDivide(qp->DZ,qp->G); CHKERRQ(info);
252:       info = qp->DS->PointwiseDivide(qp->DS,qp->T); CHKERRQ(info);
253:       info = qp->RHS2->Waxpby(1.0,qp->DZ,1.0,qp->DS); CHKERRQ(info);
254:     } else {
255:       info = qp->RHS2->SetToZero(); CHKERRQ(info);
256:     }


259:     /* 
260:        Compute the Primal Infeasiblitiy RHS and the 
261:        Diagonal Matrix to be added to H and store in Work 
262:     */
263:     info = qp->DiagAxpy->PointwiseDivide(qp->Z,qp->G); CHKERRQ(info);
264:     info = qp->GZwork->PointwiseMultiply(qp->DiagAxpy,qp->R3); CHKERRQ(info);
265:     info = qp->RHS->Axpy(-1.0,qp->GZwork); CHKERRQ(info);

267:     info = qp->TSwork->PointwiseDivide(qp->S,qp->T); CHKERRQ(info);
268:     info = qp->DiagAxpy->Axpy(1.0,qp->TSwork); CHKERRQ(info);
269:     info = qp->TSwork->PointwiseMultiply(qp->TSwork,qp->R5); CHKERRQ(info);
270:     info = qp->RHS->Axpy(-1.0,qp->TSwork); CHKERRQ(info);

272:     info = qp->RHS2->Axpy(1.0,qp->RHS); CHKERRQ(info);

274:     /*  Determine the solving tolerance */
275:     ksptol = qp->mu/10.0;
276:     ksptol = TaoMin(ksptol,0.001);

278:     info = qp->H->AddDiagonal(qp->DiagAxpy); CHKERRQ(info);

280:     info = TaoPreLinearSolve(tao,qp->H);CHKERRQ(info);
281:     info = TaoLinearSolve(tao,qp->H,qp->RHS,qp->DXY,&kspsuccess);CHKERRQ(info);
282:     

284:     info = qp->DiagAxpy->Negate(); CHKERRQ(info);
285:     info = qp->H->AddDiagonal(qp->DiagAxpy); CHKERRQ(info);
286:     info = qp->DiagAxpy->Negate(); CHKERRQ(info);
287:     info = QPComputeStepDirection(qp); CHKERRQ(info);
288:     info = QPStepLength(qp);  CHKERRQ(info);

290:     /* Calculate New Residual R1 in Work vector */
291:     info = qp->H->Multiply(qp->DXY,qp->RHS2); CHKERRQ(info);
292:     info = qp->RHS2->Axpy(1.0,qp->DS); CHKERRQ(info);
293:     info = qp->RHS2->Axpy(-1.0,qp->DZ); CHKERRQ(info);
294:     info = qp->RHS2->Aypx(qp->dsteplength,qp->R12); CHKERRQ(info);

296:     qp->RHS2->Norm2(&qp->dinfeas); CHKERRQ(info);
297:     qp->DG->Dot(qp->DZ, gap); CHKERRQ(info);
298:     qp->DT->Dot(qp->DS, gap+1); CHKERRQ(info);
299:  
300:     qp->rnorm=(qp->dinfeas+qp->psteplength*qp->pinfeas)/(qp->m+qp->n);
301:     pstep = qp->psteplength; dstep = qp->dsteplength;    
302:     step = TaoMin(qp->psteplength,qp->dsteplength);
303:     sigmamu= ( pstep*pstep*(gap[0]+gap[1]) +
304:                (1 - pstep + pstep*sigma)*qp->gap  )/qp->m;

306:     if (qp->predcorr && step < 0.9){
307:       if (sigmamu < qp->mu){ 
308:         sigmamu=sigmamu/qp->mu;
309:         sigmamu=sigmamu*sigmamu*sigmamu;
310:       } else {sigmamu = 1.0;}
311:       sigmamu = sigmamu*qp->mu;
312:       
313:       /* Compute Corrector Step */
314:       info = qp->DZ->PointwiseMultiply(qp->DG,qp->DZ); CHKERRQ(info);
315:       info = qp->DZ->Negate(); CHKERRQ(info);
316:       info = qp->DZ->AddConstant(sigmamu); CHKERRQ(info);
317:       info = qp->DZ->PointwiseDivide(qp->DZ,qp->G); CHKERRQ(info);

319:       info = qp->DS->PointwiseMultiply(qp->DS,qp->DT); CHKERRQ(info);
320:       info = qp->DS->Negate(); CHKERRQ(info);
321:       info = qp->DS->AddConstant(sigmamu); CHKERRQ(info);
322:       info = qp->DS->PointwiseDivide(qp->DS,qp->T); CHKERRQ(info);
323:       
324:       info = qp->RHS2->Waxpby(1.0,qp->DZ,-1.0,qp->DS); CHKERRQ(info);
325:       info = qp->RHS2->Axpy(1.0,qp->RHS); CHKERRQ(info);

327:       /* Approximately solve the linear system */
328:       
329:       info = qp->H->AddDiagonal(qp->DiagAxpy); CHKERRQ(info);
330:        
331:       info = TaoLinearSolve(tao,qp->H,qp->RHS2,qp->DXY,&kspsuccess);CHKERRQ(info);
332:       
333:       info = qp->H->SetDiagonal(qp->HDiag); CHKERRQ(info);
334:       info = QPComputeStepDirection(qp); CHKERRQ(info);
335:       info = QPStepLength(qp); CHKERRQ(info);

337:     }  /* End Corrector step */


340:     /* Take the step */
341:     pstep = qp->psteplength; dstep = qp->dsteplength;

343:     info = qp->Z->Axpy(dstep,qp->DZ); CHKERRQ(info);
344:     info = qp->S->Axpy(dstep,qp->DS); CHKERRQ(info);
345:     info = qp->XY->Axpy(dstep,qp->DXY); CHKERRQ(info);
346:     info = qp->G->Axpy(pstep,qp->DG); CHKERRQ(info);
347:     info = qp->T->Axpy(pstep,qp->DT); CHKERRQ(info);
348:     
349:     /* Compute Residuals */
350:     info = QPIPComputeResidual(qp); CHKERRQ(info);

352:     /* Evaluate quadratic function */
353:     info = qp->H->Multiply(qp->XY,qp->Work); CHKERRQ(info);

355:     info = qp->XY->Dot(qp->Work,&d1); CHKERRQ(info);
356:     info = qp->XY->Dot(qp->C0,&d2); CHKERRQ(info);
357:     info = qp->G->Dot(qp->Z,gap); CHKERRQ(info);
358:     info = qp->T->Dot(qp->S,gap+1); CHKERRQ(info);

360:     qp->pobj=d1/2.0 + d2+qp->c;
361:     /* Compute the duality gap */
362:     qp->gap = (gap[0]+gap[1]);
363:     qp->dobj = qp->pobj - qp->gap;
364:     if (qp->m>0) qp->mu=qp->gap/(qp->m);
365:     qp->rgap=qp->gap/( fabs(qp->dobj) + fabs(qp->pobj) + 1.0 );
366:   }  /* END MAIN LOOP  */

368:   TaoFunctionReturn(0);
369: }

373: static int QPComputeStepDirection(TAO_BQPIP *qp)
374: {
375:   int info;

377:   TaoFunctionBegin;

379:   /* Calculate DG */
380:   info = qp->DG->Waxpby(1.0,qp->DXY,1.0,qp->R3);CHKERRQ(info);

382:   /* Calculate DT */
383:   info = qp->DT->Waxpby(-1.0,qp->DXY,-1.0,qp->R5);CHKERRQ(info);

385:   /* Calculate DZ */
386:   
387:   info = qp->DZ->Axpy(-1.0,qp->Z);CHKERRQ(info);
388:   info = qp->GZwork->PointwiseDivide(qp->DG,qp->G); CHKERRQ(info);
389:   info = qp->GZwork->PointwiseMultiply(qp->GZwork,qp->Z); CHKERRQ(info);
390:   info = qp->DZ->Axpy(-1.0,qp->GZwork);CHKERRQ(info);

392:   /* Calculate DS */

394:   info = qp->DS->Axpy(-1.0,qp->S);CHKERRQ(info);
395:   info = qp->TSwork->PointwiseDivide(qp->DT,qp->T); CHKERRQ(info);
396:   info = qp->TSwork->PointwiseMultiply(qp->TSwork,qp->S); CHKERRQ(info);
397:   info = qp->DS->Axpy(-1.0,qp->TSwork);CHKERRQ(info);

399:   TaoFunctionReturn(0);
400: }

404: static int QPIPComputeResidual(TAO_BQPIP *qp)
405: {
406:   int info;
407:   double gap1,gap2,dtmp = 1.0 - qp->psteplength;

409:   TaoFunctionBegin;

411:   /* Compute R3 and R5 */

413:   if (1==1){

415:     info = qp->R3->Scale(dtmp);
416:     info = qp->R5->Scale(dtmp);
417:     qp->pinfeas=dtmp*qp->pinfeas;

419:   } else {

421:     info = qp->R3->Waxpby(1.0,qp->XY,-1.0,qp->XL);CHKERRQ(info);
422:     info = qp->R3->Axpy(-1.0,qp->G);CHKERRQ(info);

424:     info = qp->R5->Waxpby(1.0,qp->XY,-1.0,qp->XU);CHKERRQ(info);
425:     info = qp->R5->Axpy(1.0,qp->T);CHKERRQ(info);

427:     info = qp->R3->NormInfinity(&gap1);CHKERRQ(info);
428:     info = qp->R5->NormInfinity(&gap2);CHKERRQ(info);

430:     qp->pinfeas=TaoMax(gap1,gap2);

432:   }

434:   qp->R12->Waxpby(1.0,qp->S,-1.0,qp->Z);CHKERRQ(info);

436:   info = qp->H->Multiply(qp->XY,qp->RHS); CHKERRQ(info);
437:   info = qp->RHS->Negate();CHKERRQ(info);
438:   info = qp->RHS->Axpy(-1.0,qp->C0);CHKERRQ(info);
439:   info = qp->R12->Axpy(-1.0,qp->RHS);CHKERRQ(info);

441:   info = qp->R12->Norm1(&qp->dinfeas); CHKERRQ(info);
442:   qp->rnorm=(qp->dinfeas+qp->pinfeas)/(qp->m+qp->n);
443:   
444:   TaoFunctionReturn(0);
445: }

449: static int QPStepLength(TAO_BQPIP *qp)
450: {
451:   double tstep1,tstep2,tstep3,tstep4,tstep;
452:   int info;

454:   TaoFunctionBegin;
455:   /* Compute stepsize to the boundary */

457:   info = qp->G->StepMax(qp->DG,&tstep1); CHKERRQ(info);
458:   info = qp->T->StepMax(qp->DT,&tstep2); CHKERRQ(info);
459:   info = qp->S->StepMax(qp->DS,&tstep3); CHKERRQ(info);
460:   info = qp->Z->StepMax(qp->DZ,&tstep4); CHKERRQ(info);

462:   tstep = TaoMin(tstep1,tstep2);
463:   qp->psteplength = TaoMin(0.95*tstep,1.0);

465:   tstep = TaoMin(tstep3,tstep4);
466:   qp->dsteplength = TaoMin(0.95*tstep,1.0);

468:   qp->psteplength = TaoMin(qp->psteplength,qp->dsteplength);
469:   qp->dsteplength = qp->psteplength;

471:   TaoFunctionReturn(0);
472: }


477: int TaoGetDualVariables_BQPIP(TAO_SOLVER tao,TaoVec* DXL, TaoVec* DXU, void *solver)
478: {
479:   TAO_BQPIP *qp = (TAO_BQPIP*)solver;
480:   int       info;

482:   TaoFunctionBegin;
483:   info = DXL->CopyFrom(qp->Z); CHKERRQ(info);
484:   info = DXU->CopyFrom(qp->S); CHKERRQ(info);
485:   info = DXU->Negate(); CHKERRQ(info);
486:   TaoFunctionReturn(0);
487: }


492: int TAOComputeNormFromCentralPath_BQPIP(TAO_SOLVER tao, double *norm)
493: {
494:   TAO_BQPIP *qp;
495:   int       info;
496:   double    gap[2],mu[2], nmu;
497:   
498:   TaoFunctionBegin;
499:   info = TaoGetSolverContext(tao,"tao_bqpip",(void**)&qp); CHKERRQ(info);
500:   info = qp->GZwork->PointwiseMultiply(qp->G,qp->Z); CHKERRQ(info);
501:   info = qp->TSwork->PointwiseMultiply(qp->T,qp->S); CHKERRQ(info);
502:   info = qp->TSwork->Norm1(&mu[0]); CHKERRQ(info);
503:   info = qp->GZwork->Norm1(&mu[1]); CHKERRQ(info);

505:   nmu=-(mu[0]+mu[1])/qp->m;

507:   qp->GZwork->AddConstant(nmu); CHKERRQ(info);
508:   qp->TSwork->AddConstant(nmu); CHKERRQ(info);

510:   qp->GZwork->Norm2squared(&gap[0]); CHKERRQ(info);
511:   qp->TSwork->Norm2squared(&gap[1]); CHKERRQ(info);

513:   qp->pathnorm=sqrt( (gap[0]+gap[1]) );
514:   *norm=qp->pathnorm;

516:   TaoFunctionReturn(0);
517: }


522: static int TaoSetOptions_BQPIP(TAO_SOLVER tao, void *solver){

524:   TAO_BQPIP *qp = (TAO_BQPIP*)solver;
525:   int       info;
526:   TaoTruth  flg;

528:   TaoFunctionBegin;
529:   info = TaoOptionsHead("Interior point method for bound constrained quadratic optimization");CHKERRQ(info);
530:   info = TaoOptionInt("-predcorr","Use a predictor-corrector method","",qp->predcorr,&qp->predcorr,&flg);
531:   CHKERRQ(info);
532:   info = TaoOptionsTail();CHKERRQ(info);

534:   /*   info = TaoSetLinearSolverOptions(tao);CHKERRQ(info); */

536:   TaoFunctionReturn(0);
537: }


542: static int TaoView_BQPIP(TAO_SOLVER tao,void *solver){
543:   TaoFunctionBegin;
544:   TaoFunctionReturn(0);
545: }


548: /* --------------------------------------------------------- */
552: int TaoCreate_BQPIP(TAO_SOLVER tao)
553: {
554:   TAO_BQPIP *qp;
555:   int       info;

557:   TaoFunctionBegin;

559:   info = TaoNew(TAO_BQPIP,&qp); CHKERRQ(info);
560:   info = PetscLogObjectMemory(tao,sizeof(TAO_BQPIP)); CHKERRQ(info);

562:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_BQPIP,(void*)qp); CHKERRQ(info);
563:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_BQPIP,TaoSetDown_BQPIP); CHKERRQ(info);
564:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_BQPIP); CHKERRQ(info);
565:   info=TaoSetTaoViewRoutine(tao,TaoView_BQPIP); CHKERRQ(info);
566:   info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_BQPIP); CHKERRQ(info);

568:   info = TaoSetMaximumIterates(tao,100); CHKERRQ(info);
569:   info = TaoSetMaximumFunctionEvaluations(tao,500); CHKERRQ(info);
570:   info = TaoSetTolerances(tao,1e-12,1e-12,1e-12,0); CHKERRQ(info);
571:   /*
572:   tao->defaultmonitor     = TaoDefaultMonitor_BQPIP;
573:   */
574:   /* Initialize pointers and variables */
575:   qp->n              = 0;
576:   qp->m              = 0;
577:   qp->ksp_tol       = 0.1;
578:   qp->dobj           = 0.0;
579:   qp->pobj           = 1.0;
580:   qp->gap            = 10.0;
581:   qp->rgap           = 1.0;
582:   qp->mu             = 1.0;
583:   qp->sigma          = 1.0;
584:   qp->dinfeas        = 1.0;
585:   qp->predcorr       = 1;
586:   qp->psteplength    = 0.0;
587:   qp->dsteplength    = 0.0;

589:   TaoFunctionReturn(0);
590: }