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 PetscCall(VecScale(qp->R3,dtmp)); 18 PetscCall(VecScale(qp->R5,dtmp)); 19 qp->pinfeas=dtmp*qp->pinfeas; 20 21 PetscCall(VecCopy(qp->S,tao->gradient)); 22 PetscCall(VecAXPY(tao->gradient,-1.0,qp->Z)); 23 24 PetscCall(MatMult(tao->hessian,tao->solution,qp->RHS)); 25 PetscCall(VecScale(qp->RHS,-1.0)); 26 PetscCall(VecAXPY(qp->RHS,-1.0,qp->C)); 27 PetscCall(VecAXPY(tao->gradient,-1.0,qp->RHS)); 28 29 PetscCall(VecNorm(tao->gradient,NORM_1,&qp->dinfeas)); 30 qp->rnorm=(qp->dinfeas+qp->pinfeas)/(qp->m+qp->n); 31 PetscFunctionReturn(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 PetscCall(MatMult(tao->hessian,tao->solution,tao->gradient)); 42 PetscCall(VecCopy(qp->C,qp->Work)); 43 PetscCall(VecAXPY(qp->Work,0.5,tao->gradient)); 44 PetscCall(VecAXPY(tao->gradient,1.0,qp->C)); 45 PetscCall(VecDot(tao->solution,qp->Work,&fff)); 46 qp->pobj = fff + qp->d; 47 48 PetscCheck(!PetscIsInfOrNanReal(qp->pobj),PETSC_COMM_SELF,PETSC_ERR_USER, "User provided data contains Inf or NaN"); 49 50 /* Initialize slack vectors */ 51 /* T = XU - X; G = X - XL */ 52 PetscCall(VecCopy(qp->XU,qp->T)); 53 PetscCall(VecAXPY(qp->T,-1.0,tao->solution)); 54 PetscCall(VecCopy(tao->solution,qp->G)); 55 PetscCall(VecAXPY(qp->G,-1.0,qp->XL)); 56 57 PetscCall(VecSet(qp->GZwork,p01)); 58 PetscCall(VecSet(qp->TSwork,p01)); 59 60 PetscCall(VecPointwiseMax(qp->G,qp->G,qp->GZwork)); 61 PetscCall(VecPointwiseMax(qp->T,qp->T,qp->TSwork)); 62 63 /* Initialize Dual Variable Vectors */ 64 PetscCall(VecCopy(qp->G,qp->Z)); 65 PetscCall(VecReciprocal(qp->Z)); 66 67 PetscCall(VecCopy(qp->T,qp->S)); 68 PetscCall(VecReciprocal(qp->S)); 69 70 PetscCall(MatMult(tao->hessian,qp->Work,qp->RHS)); 71 PetscCall(VecAbs(qp->RHS)); 72 PetscCall(VecSet(qp->Work,p01)); 73 PetscCall(VecPointwiseMax(qp->RHS,qp->RHS,qp->Work)); 74 75 PetscCall(VecPointwiseDivide(qp->RHS,tao->gradient,qp->RHS)); 76 PetscCall(VecNorm(qp->RHS,NORM_1,&gap1)); 77 mu = PetscMin(10.0,(gap1+10.0)/qp->m); 78 79 PetscCall(VecScale(qp->S,mu)); 80 PetscCall(VecScale(qp->Z,mu)); 81 82 PetscCall(VecSet(qp->TSwork,p01)); 83 PetscCall(VecSet(qp->GZwork,p01)); 84 PetscCall(VecPointwiseMax(qp->S,qp->S,qp->TSwork)); 85 PetscCall(VecPointwiseMax(qp->Z,qp->Z,qp->GZwork)); 86 87 qp->mu=0;qp->dinfeas=1.0;qp->pinfeas=1.0; 88 while ((qp->dinfeas+qp->pinfeas)/(qp->m+qp->n) >= qp->mu) { 89 PetscCall(VecScale(qp->G,two)); 90 PetscCall(VecScale(qp->Z,two)); 91 PetscCall(VecScale(qp->S,two)); 92 PetscCall(VecScale(qp->T,two)); 93 94 PetscCall(QPIPComputeResidual(qp,tao)); 95 96 PetscCall(VecCopy(tao->solution,qp->R3)); 97 PetscCall(VecAXPY(qp->R3,-1.0,qp->G)); 98 PetscCall(VecAXPY(qp->R3,-1.0,qp->XL)); 99 100 PetscCall(VecCopy(tao->solution,qp->R5)); 101 PetscCall(VecAXPY(qp->R5,1.0,qp->T)); 102 PetscCall(VecAXPY(qp->R5,-1.0,qp->XU)); 103 104 PetscCall(VecNorm(qp->R3,NORM_INFINITY,&gap1)); 105 PetscCall(VecNorm(qp->R5,NORM_INFINITY,&gap2)); 106 qp->pinfeas=PetscMax(gap1,gap2); 107 108 /* Compute the duality gap */ 109 PetscCall(VecDot(qp->G,qp->Z,&gap1)); 110 PetscCall(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 PetscCall(VecStepMax(qp->G,qp->DG,&tstep1)); 131 PetscCall(VecStepMax(qp->T,qp->DT,&tstep2)); 132 PetscCall(VecStepMax(qp->S,qp->DS,&tstep3)); 133 PetscCall(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 PetscCall(VecPointwiseMult(qp->GZwork,qp->G,qp->Z)); 152 PetscCall(VecPointwiseMult(qp->TSwork,qp->T,qp->S)); 153 PetscCall(VecNorm(qp->TSwork,NORM_1,&mu[0])); 154 PetscCall(VecNorm(qp->GZwork,NORM_1,&mu[1])); 155 156 nmu=-(mu[0]+mu[1])/qp->m; 157 158 PetscCall(VecShift(qp->GZwork,nmu)); 159 PetscCall(VecShift(qp->TSwork,nmu)); 160 161 PetscCall(VecNorm(qp->GZwork,NORM_2,&gap[0])); 162 PetscCall(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 PetscCall(VecCopy(tao->stepdirection,qp->DG)); 176 PetscCall(VecAXPY(qp->DG,1.0,qp->R3)); 177 178 /* Calculate DT */ 179 PetscCall(VecCopy(tao->stepdirection,qp->DT)); 180 PetscCall(VecScale(qp->DT,-1.0)); 181 PetscCall(VecAXPY(qp->DT,-1.0,qp->R5)); 182 183 /* Calculate DZ */ 184 PetscCall(VecAXPY(qp->DZ,-1.0,qp->Z)); 185 PetscCall(VecPointwiseDivide(qp->GZwork,qp->DG,qp->G)); 186 PetscCall(VecPointwiseMult(qp->GZwork,qp->GZwork,qp->Z)); 187 PetscCall(VecAXPY(qp->DZ,-1.0,qp->GZwork)); 188 189 /* Calculate DS */ 190 PetscCall(VecAXPY(qp->DS,-1.0,qp->S)); 191 PetscCall(VecPointwiseDivide(qp->TSwork,qp->DT,qp->T)); 192 PetscCall(VecPointwiseMult(qp->TSwork,qp->TSwork,qp->S)); 193 PetscCall(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 PetscCall(VecGetSize(tao->solution,&qp->n)); 204 205 /* Allocate some arrays */ 206 if (!tao->gradient) { 207 PetscCall(VecDuplicate(tao->solution,&tao->gradient)); 208 } 209 if (!tao->stepdirection) { 210 PetscCall(VecDuplicate(tao->solution,&tao->stepdirection)); 211 } 212 if (!tao->XL) { 213 PetscCall(VecDuplicate(tao->solution,&tao->XL)); 214 PetscCall(VecSet(tao->XL,PETSC_NINFINITY)); 215 } 216 if (!tao->XU) { 217 PetscCall(VecDuplicate(tao->solution,&tao->XU)); 218 PetscCall(VecSet(tao->XU,PETSC_INFINITY)); 219 } 220 221 PetscCall(VecDuplicate(tao->solution,&qp->Work)); 222 PetscCall(VecDuplicate(tao->solution,&qp->XU)); 223 PetscCall(VecDuplicate(tao->solution,&qp->XL)); 224 PetscCall(VecDuplicate(tao->solution,&qp->HDiag)); 225 PetscCall(VecDuplicate(tao->solution,&qp->DiagAxpy)); 226 PetscCall(VecDuplicate(tao->solution,&qp->RHS)); 227 PetscCall(VecDuplicate(tao->solution,&qp->RHS2)); 228 PetscCall(VecDuplicate(tao->solution,&qp->C)); 229 230 PetscCall(VecDuplicate(tao->solution,&qp->G)); 231 PetscCall(VecDuplicate(tao->solution,&qp->DG)); 232 PetscCall(VecDuplicate(tao->solution,&qp->S)); 233 PetscCall(VecDuplicate(tao->solution,&qp->Z)); 234 PetscCall(VecDuplicate(tao->solution,&qp->DZ)); 235 PetscCall(VecDuplicate(tao->solution,&qp->GZwork)); 236 PetscCall(VecDuplicate(tao->solution,&qp->R3)); 237 238 PetscCall(VecDuplicate(tao->solution,&qp->T)); 239 PetscCall(VecDuplicate(tao->solution,&qp->DT)); 240 PetscCall(VecDuplicate(tao->solution,&qp->DS)); 241 PetscCall(VecDuplicate(tao->solution,&qp->TSwork)); 242 PetscCall(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 PetscCall(TaoComputeVariableBounds(tao)); 277 PetscCall(VecSet(qp->XU,1.0e20)); 278 PetscCall(VecSet(qp->XL,-1.0e20)); 279 PetscCall(VecPointwiseMax(qp->XL,qp->XL,tao->XL)); 280 PetscCall(VecPointwiseMin(qp->XU,qp->XU,tao->XU)); 281 PetscCall(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 PetscCall(VecSet(qp->Work,0)); 287 PetscCall(TaoComputeObjectiveAndGradient(tao,qp->Work,&qp->d,qp->C)); 288 PetscCall(TaoComputeHessian(tao,qp->Work,tao->hessian,tao->hessian_pre)); 289 PetscCall(MatHasOperation(tao->hessian,MATOP_GET_DIAGONAL,&getdiagop)); 290 if (getdiagop) { 291 PetscCall(MatGetDiagonal(tao->hessian,qp->HDiag)); 292 } 293 294 /* Initialize starting point and residuals */ 295 PetscCall(QPIPSetInitialPoint(qp,tao)); 296 PetscCall(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 PetscCall(TaoLogConvergenceHistory(tao,qp->pobj,gnorm,qp->pinfeas,tao->ksp_its)); 305 PetscCall(TaoMonitor(tao,tao->niter,qp->pobj,gnorm,qp->pinfeas,step)); 306 PetscCall((*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 PetscCall((*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 PetscCall(QPIPComputeNormFromCentralPath(qp,&d1)); 321 322 PetscCall(VecSet(qp->DZ,0.0)); 323 PetscCall(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 PetscCall(VecPointwiseDivide(qp->DiagAxpy,qp->Z,qp->G)); 330 PetscCall(VecPointwiseMult(qp->GZwork,qp->DiagAxpy,qp->R3)); 331 PetscCall(VecAXPY(qp->RHS,-1.0,qp->GZwork)); 332 333 PetscCall(VecPointwiseDivide(qp->TSwork,qp->S,qp->T)); 334 PetscCall(VecAXPY(qp->DiagAxpy,1.0,qp->TSwork)); 335 PetscCall(VecPointwiseMult(qp->TSwork,qp->TSwork,qp->R5)); 336 PetscCall(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 PetscCall(KSPSetTolerances(tao->ksp,ksptol,1e-30,1e30,PetscMax(10,qp->n))); 342 343 /* Shift the diagonals of the Hessian matrix */ 344 PetscCall(MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES)); 345 if (!getdiagop) { 346 PetscCall(VecCopy(qp->DiagAxpy,qp->HDiag)); 347 PetscCall(VecScale(qp->HDiag,-1.0)); 348 } 349 PetscCall(MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY)); 350 PetscCall(MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY)); 351 352 PetscCall(KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre)); 353 PetscCall(KSPSolve(tao->ksp,qp->RHS,tao->stepdirection)); 354 PetscCall(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 PetscCall(MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES)); 361 } else { 362 PetscCall(MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES)); 363 } 364 PetscCall(MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY)); 365 PetscCall(MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY)); 366 PetscCall(QPIPComputeStepDirection(qp,tao)); 367 PetscCall(QPIPStepLength(qp)); 368 369 /* Calculate New Residual R1 in Work vector */ 370 PetscCall(MatMult(tao->hessian,tao->stepdirection,qp->RHS2)); 371 PetscCall(VecAXPY(qp->RHS2,1.0,qp->DS)); 372 PetscCall(VecAXPY(qp->RHS2,-1.0,qp->DZ)); 373 PetscCall(VecAYPX(qp->RHS2,qp->dsteplength,tao->gradient)); 374 375 PetscCall(VecNorm(qp->RHS2,NORM_2,&qp->dinfeas)); 376 PetscCall(VecDot(qp->DZ,qp->DG,gap)); 377 PetscCall(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 PetscCall(VecPointwiseMult(qp->DZ,qp->DG,qp->DZ)); 395 PetscCall(VecScale(qp->DZ,-1.0)); 396 PetscCall(VecShift(qp->DZ,sigmamu)); 397 PetscCall(VecPointwiseDivide(qp->DZ,qp->DZ,qp->G)); 398 399 PetscCall(VecPointwiseMult(qp->DS,qp->DS,qp->DT)); 400 PetscCall(VecScale(qp->DS,-1.0)); 401 PetscCall(VecShift(qp->DS,sigmamu)); 402 PetscCall(VecPointwiseDivide(qp->DS,qp->DS,qp->T)); 403 404 PetscCall(VecCopy(qp->DZ,qp->RHS2)); 405 PetscCall(VecAXPY(qp->RHS2,-1.0,qp->DS)); 406 PetscCall(VecAXPY(qp->RHS2,1.0,qp->RHS)); 407 408 /* Approximately solve the linear system */ 409 PetscCall(MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES)); 410 if (!getdiagop) { 411 PetscCall(VecCopy(qp->DiagAxpy,qp->HDiag)); 412 PetscCall(VecScale(qp->HDiag,-1.0)); 413 } 414 PetscCall(MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY)); 415 PetscCall(MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY)); 416 417 /* Solve using the previous tolerances that were set */ 418 PetscCall(KSPSolve(tao->ksp,qp->RHS2,tao->stepdirection)); 419 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 420 tao->ksp_its += its; 421 tao->ksp_tot_its += its; 422 423 if (getdiagop) { 424 PetscCall(MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES)); 425 } else { 426 PetscCall(MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES)); 427 } 428 PetscCall(MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY)); 429 PetscCall(MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY)); 430 PetscCall(QPIPComputeStepDirection(qp,tao)); 431 PetscCall(QPIPStepLength(qp)); 432 } /* End Corrector step */ 433 434 /* Take the step */ 435 dstep = qp->dsteplength; 436 437 PetscCall(VecAXPY(qp->Z,dstep,qp->DZ)); 438 PetscCall(VecAXPY(qp->S,dstep,qp->DS)); 439 PetscCall(VecAXPY(tao->solution,dstep,tao->stepdirection)); 440 PetscCall(VecAXPY(qp->G,dstep,qp->DG)); 441 PetscCall(VecAXPY(qp->T,dstep,qp->DT)); 442 443 /* Compute Residuals */ 444 PetscCall(QPIPComputeResidual(qp,tao)); 445 446 /* Evaluate quadratic function */ 447 PetscCall(MatMult(tao->hessian,tao->solution,qp->Work)); 448 449 PetscCall(VecDot(tao->solution,qp->Work,&d1)); 450 PetscCall(VecDot(tao->solution,qp->C,&d2)); 451 PetscCall(VecDot(qp->G,qp->Z,gap)); 452 PetscCall(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 PetscCall(PetscOptionsHead(PetscOptionsObject,"Interior point method for bound constrained quadratic optimization")); 478 PetscCall(PetscOptionsInt("-tao_bqpip_predcorr","Use a predictor-corrector method","",qp->predcorr,&qp->predcorr,NULL)); 479 PetscCall(PetscOptionsTail()); 480 PetscCall(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 PetscCall(VecDestroy(&qp->G)); 491 PetscCall(VecDestroy(&qp->DG)); 492 PetscCall(VecDestroy(&qp->Z)); 493 PetscCall(VecDestroy(&qp->DZ)); 494 PetscCall(VecDestroy(&qp->GZwork)); 495 PetscCall(VecDestroy(&qp->R3)); 496 PetscCall(VecDestroy(&qp->S)); 497 PetscCall(VecDestroy(&qp->DS)); 498 PetscCall(VecDestroy(&qp->T)); 499 500 PetscCall(VecDestroy(&qp->DT)); 501 PetscCall(VecDestroy(&qp->TSwork)); 502 PetscCall(VecDestroy(&qp->R5)); 503 PetscCall(VecDestroy(&qp->HDiag)); 504 PetscCall(VecDestroy(&qp->Work)); 505 PetscCall(VecDestroy(&qp->XL)); 506 PetscCall(VecDestroy(&qp->XU)); 507 PetscCall(VecDestroy(&qp->DiagAxpy)); 508 PetscCall(VecDestroy(&qp->RHS)); 509 PetscCall(VecDestroy(&qp->RHS2)); 510 PetscCall(VecDestroy(&qp->C)); 511 } 512 PetscCall(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 PetscCall(VecCopy(qp->Z,DXL)); 522 PetscCall(VecCopy(qp->S,DXU)); 523 PetscCall(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 PetscCall(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 PetscCall(KSPCreate(((PetscObject)tao)->comm,&tao->ksp)); 572 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 573 PetscCall(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 574 PetscCall(KSPSetType(tao->ksp,KSPCG)); 575 PetscCall(KSPSetTolerances(tao->ksp,1e-14,1e-30,1e30,PetscMax(10,qp->n))); 576 PetscFunctionReturn(0); 577 } 578