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