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 213 PetscCall(VecDuplicate(tao->solution,&qp->Work)); 214 PetscCall(VecDuplicate(tao->solution,&qp->XU)); 215 PetscCall(VecDuplicate(tao->solution,&qp->XL)); 216 PetscCall(VecDuplicate(tao->solution,&qp->HDiag)); 217 PetscCall(VecDuplicate(tao->solution,&qp->DiagAxpy)); 218 PetscCall(VecDuplicate(tao->solution,&qp->RHS)); 219 PetscCall(VecDuplicate(tao->solution,&qp->RHS2)); 220 PetscCall(VecDuplicate(tao->solution,&qp->C)); 221 222 PetscCall(VecDuplicate(tao->solution,&qp->G)); 223 PetscCall(VecDuplicate(tao->solution,&qp->DG)); 224 PetscCall(VecDuplicate(tao->solution,&qp->S)); 225 PetscCall(VecDuplicate(tao->solution,&qp->Z)); 226 PetscCall(VecDuplicate(tao->solution,&qp->DZ)); 227 PetscCall(VecDuplicate(tao->solution,&qp->GZwork)); 228 PetscCall(VecDuplicate(tao->solution,&qp->R3)); 229 230 PetscCall(VecDuplicate(tao->solution,&qp->T)); 231 PetscCall(VecDuplicate(tao->solution,&qp->DT)); 232 PetscCall(VecDuplicate(tao->solution,&qp->DS)); 233 PetscCall(VecDuplicate(tao->solution,&qp->TSwork)); 234 PetscCall(VecDuplicate(tao->solution,&qp->R5)); 235 qp->m=2*qp->n; 236 PetscFunctionReturn(0); 237 } 238 239 static PetscErrorCode TaoSolve_BQPIP(Tao tao) 240 { 241 TAO_BQPIP *qp = (TAO_BQPIP*)tao->data; 242 PetscInt its; 243 PetscReal d1,d2,ksptol,sigmamu; 244 PetscReal gnorm,dstep,pstep,step=0; 245 PetscReal gap[4]; 246 PetscBool getdiagop; 247 248 PetscFunctionBegin; 249 qp->dobj = 0.0; 250 qp->pobj = 1.0; 251 qp->gap = 10.0; 252 qp->rgap = 1.0; 253 qp->mu = 1.0; 254 qp->dinfeas = 1.0; 255 qp->psteplength = 0.0; 256 qp->dsteplength = 0.0; 257 258 /* TODO 259 - Remove fixed variables and treat them correctly 260 - Use index sets for the infinite versus finite bounds 261 - Update remaining code for fixed and free variables 262 - Fix inexact solves for predictor and corrector 263 */ 264 265 /* Tighten infinite bounds, things break when we don't do this 266 -- see test_bqpip.c 267 */ 268 PetscCall(TaoComputeVariableBounds(tao)); 269 PetscCall(VecSet(qp->XU,1.0e20)); 270 PetscCall(VecSet(qp->XL,-1.0e20)); 271 if (tao->XL) PetscCall(VecPointwiseMax(qp->XL,qp->XL,tao->XL)); 272 if (tao->XU) PetscCall(VecPointwiseMin(qp->XU,qp->XU,tao->XU)); 273 PetscCall(VecMedian(qp->XL,tao->solution,qp->XU,tao->solution)); 274 275 /* Evaluate gradient and Hessian at zero to get the correct values 276 without contaminating them with numerical artifacts. 277 */ 278 PetscCall(VecSet(qp->Work,0)); 279 PetscCall(TaoComputeObjectiveAndGradient(tao,qp->Work,&qp->d,qp->C)); 280 PetscCall(TaoComputeHessian(tao,qp->Work,tao->hessian,tao->hessian_pre)); 281 PetscCall(MatHasOperation(tao->hessian,MATOP_GET_DIAGONAL,&getdiagop)); 282 if (getdiagop) PetscCall(MatGetDiagonal(tao->hessian,qp->HDiag)); 283 284 /* Initialize starting point and residuals */ 285 PetscCall(QPIPSetInitialPoint(qp,tao)); 286 PetscCall(QPIPComputeResidual(qp,tao)); 287 288 /* Enter main loop */ 289 tao->reason = TAO_CONTINUE_ITERATING; 290 while (1) { 291 292 /* Check Stopping Condition */ 293 gnorm = PetscSqrtScalar(qp->gap + qp->dinfeas); 294 PetscCall(TaoLogConvergenceHistory(tao,qp->pobj,gnorm,qp->pinfeas,tao->ksp_its)); 295 PetscCall(TaoMonitor(tao,tao->niter,qp->pobj,gnorm,qp->pinfeas,step)); 296 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 297 if (tao->reason != TAO_CONTINUE_ITERATING) break; 298 /* Call general purpose update function */ 299 if (tao->ops->update) PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 300 tao->niter++; 301 tao->ksp_its = 0; 302 303 /* 304 Dual Infeasibility Direction should already be in the right 305 hand side from computing the residuals 306 */ 307 308 PetscCall(QPIPComputeNormFromCentralPath(qp,&d1)); 309 310 PetscCall(VecSet(qp->DZ,0.0)); 311 PetscCall(VecSet(qp->DS,0.0)); 312 313 /* 314 Compute the Primal Infeasiblitiy RHS and the 315 Diagonal Matrix to be added to H and store in Work 316 */ 317 PetscCall(VecPointwiseDivide(qp->DiagAxpy,qp->Z,qp->G)); 318 PetscCall(VecPointwiseMult(qp->GZwork,qp->DiagAxpy,qp->R3)); 319 PetscCall(VecAXPY(qp->RHS,-1.0,qp->GZwork)); 320 321 PetscCall(VecPointwiseDivide(qp->TSwork,qp->S,qp->T)); 322 PetscCall(VecAXPY(qp->DiagAxpy,1.0,qp->TSwork)); 323 PetscCall(VecPointwiseMult(qp->TSwork,qp->TSwork,qp->R5)); 324 PetscCall(VecAXPY(qp->RHS,-1.0,qp->TSwork)); 325 326 /* Determine the solving tolerance */ 327 ksptol = qp->mu/10.0; 328 ksptol = PetscMin(ksptol,0.001); 329 PetscCall(KSPSetTolerances(tao->ksp,ksptol,1e-30,1e30,PetscMax(10,qp->n))); 330 331 /* Shift the diagonals of the Hessian matrix */ 332 PetscCall(MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES)); 333 if (!getdiagop) { 334 PetscCall(VecCopy(qp->DiagAxpy,qp->HDiag)); 335 PetscCall(VecScale(qp->HDiag,-1.0)); 336 } 337 PetscCall(MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY)); 338 PetscCall(MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY)); 339 340 PetscCall(KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre)); 341 PetscCall(KSPSolve(tao->ksp,qp->RHS,tao->stepdirection)); 342 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 343 tao->ksp_its += its; 344 tao->ksp_tot_its += its; 345 346 /* Restore the true diagonal of the Hessian matrix */ 347 if (getdiagop) { 348 PetscCall(MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES)); 349 } else { 350 PetscCall(MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES)); 351 } 352 PetscCall(MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY)); 353 PetscCall(MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY)); 354 PetscCall(QPIPComputeStepDirection(qp,tao)); 355 PetscCall(QPIPStepLength(qp)); 356 357 /* Calculate New Residual R1 in Work vector */ 358 PetscCall(MatMult(tao->hessian,tao->stepdirection,qp->RHS2)); 359 PetscCall(VecAXPY(qp->RHS2,1.0,qp->DS)); 360 PetscCall(VecAXPY(qp->RHS2,-1.0,qp->DZ)); 361 PetscCall(VecAYPX(qp->RHS2,qp->dsteplength,tao->gradient)); 362 363 PetscCall(VecNorm(qp->RHS2,NORM_2,&qp->dinfeas)); 364 PetscCall(VecDot(qp->DZ,qp->DG,gap)); 365 PetscCall(VecDot(qp->DS,qp->DT,gap+1)); 366 367 qp->rnorm = (qp->dinfeas+qp->psteplength*qp->pinfeas)/(qp->m+qp->n); 368 pstep = qp->psteplength; 369 step = PetscMin(qp->psteplength,qp->dsteplength); 370 sigmamu = (pstep*pstep*(gap[0]+gap[1]) + (1 - pstep)*qp->gap)/qp->m; 371 372 if (qp->predcorr && step < 0.9) { 373 if (sigmamu < qp->mu) { 374 sigmamu = sigmamu/qp->mu; 375 sigmamu = sigmamu*sigmamu*sigmamu; 376 } else { 377 sigmamu = 1.0; 378 } 379 sigmamu = sigmamu*qp->mu; 380 381 /* Compute Corrector Step */ 382 PetscCall(VecPointwiseMult(qp->DZ,qp->DG,qp->DZ)); 383 PetscCall(VecScale(qp->DZ,-1.0)); 384 PetscCall(VecShift(qp->DZ,sigmamu)); 385 PetscCall(VecPointwiseDivide(qp->DZ,qp->DZ,qp->G)); 386 387 PetscCall(VecPointwiseMult(qp->DS,qp->DS,qp->DT)); 388 PetscCall(VecScale(qp->DS,-1.0)); 389 PetscCall(VecShift(qp->DS,sigmamu)); 390 PetscCall(VecPointwiseDivide(qp->DS,qp->DS,qp->T)); 391 392 PetscCall(VecCopy(qp->DZ,qp->RHS2)); 393 PetscCall(VecAXPY(qp->RHS2,-1.0,qp->DS)); 394 PetscCall(VecAXPY(qp->RHS2,1.0,qp->RHS)); 395 396 /* Approximately solve the linear system */ 397 PetscCall(MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES)); 398 if (!getdiagop) { 399 PetscCall(VecCopy(qp->DiagAxpy,qp->HDiag)); 400 PetscCall(VecScale(qp->HDiag,-1.0)); 401 } 402 PetscCall(MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY)); 403 PetscCall(MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY)); 404 405 /* Solve using the previous tolerances that were set */ 406 PetscCall(KSPSolve(tao->ksp,qp->RHS2,tao->stepdirection)); 407 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 408 tao->ksp_its += its; 409 tao->ksp_tot_its += its; 410 411 if (getdiagop) { 412 PetscCall(MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES)); 413 } else { 414 PetscCall(MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES)); 415 } 416 PetscCall(MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY)); 417 PetscCall(MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY)); 418 PetscCall(QPIPComputeStepDirection(qp,tao)); 419 PetscCall(QPIPStepLength(qp)); 420 } /* End Corrector step */ 421 422 /* Take the step */ 423 dstep = qp->dsteplength; 424 425 PetscCall(VecAXPY(qp->Z,dstep,qp->DZ)); 426 PetscCall(VecAXPY(qp->S,dstep,qp->DS)); 427 PetscCall(VecAXPY(tao->solution,dstep,tao->stepdirection)); 428 PetscCall(VecAXPY(qp->G,dstep,qp->DG)); 429 PetscCall(VecAXPY(qp->T,dstep,qp->DT)); 430 431 /* Compute Residuals */ 432 PetscCall(QPIPComputeResidual(qp,tao)); 433 434 /* Evaluate quadratic function */ 435 PetscCall(MatMult(tao->hessian,tao->solution,qp->Work)); 436 437 PetscCall(VecDot(tao->solution,qp->Work,&d1)); 438 PetscCall(VecDot(tao->solution,qp->C,&d2)); 439 PetscCall(VecDot(qp->G,qp->Z,gap)); 440 PetscCall(VecDot(qp->T,qp->S,gap+1)); 441 442 /* Compute the duality gap */ 443 qp->pobj = d1/2.0 + d2+qp->d; 444 qp->gap = gap[0]+gap[1]; 445 qp->dobj = qp->pobj - qp->gap; 446 if (qp->m > 0) { 447 qp->mu = qp->gap/(qp->m); 448 } 449 qp->rgap = qp->gap/(PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0); 450 } /* END MAIN LOOP */ 451 PetscFunctionReturn(0); 452 } 453 454 static PetscErrorCode TaoView_BQPIP(Tao tao,PetscViewer viewer) 455 { 456 PetscFunctionBegin; 457 PetscFunctionReturn(0); 458 } 459 460 static PetscErrorCode TaoSetFromOptions_BQPIP(PetscOptionItems *PetscOptionsObject,Tao tao) 461 { 462 TAO_BQPIP *qp = (TAO_BQPIP*)tao->data; 463 464 PetscFunctionBegin; 465 PetscOptionsHeadBegin(PetscOptionsObject,"Interior point method for bound constrained quadratic optimization"); 466 PetscCall(PetscOptionsInt("-tao_bqpip_predcorr","Use a predictor-corrector method","",qp->predcorr,&qp->predcorr,NULL)); 467 PetscOptionsHeadEnd(); 468 PetscCall(KSPSetFromOptions(tao->ksp)); 469 PetscFunctionReturn(0); 470 } 471 472 static PetscErrorCode TaoDestroy_BQPIP(Tao tao) 473 { 474 TAO_BQPIP *qp = (TAO_BQPIP*)tao->data; 475 476 PetscFunctionBegin; 477 if (tao->setupcalled) { 478 PetscCall(VecDestroy(&qp->G)); 479 PetscCall(VecDestroy(&qp->DG)); 480 PetscCall(VecDestroy(&qp->Z)); 481 PetscCall(VecDestroy(&qp->DZ)); 482 PetscCall(VecDestroy(&qp->GZwork)); 483 PetscCall(VecDestroy(&qp->R3)); 484 PetscCall(VecDestroy(&qp->S)); 485 PetscCall(VecDestroy(&qp->DS)); 486 PetscCall(VecDestroy(&qp->T)); 487 488 PetscCall(VecDestroy(&qp->DT)); 489 PetscCall(VecDestroy(&qp->TSwork)); 490 PetscCall(VecDestroy(&qp->R5)); 491 PetscCall(VecDestroy(&qp->HDiag)); 492 PetscCall(VecDestroy(&qp->Work)); 493 PetscCall(VecDestroy(&qp->XL)); 494 PetscCall(VecDestroy(&qp->XU)); 495 PetscCall(VecDestroy(&qp->DiagAxpy)); 496 PetscCall(VecDestroy(&qp->RHS)); 497 PetscCall(VecDestroy(&qp->RHS2)); 498 PetscCall(VecDestroy(&qp->C)); 499 } 500 PetscCall(KSPDestroy(&tao->ksp)); 501 PetscCall(PetscFree(tao->data)); 502 PetscFunctionReturn(0); 503 } 504 505 static PetscErrorCode TaoComputeDual_BQPIP(Tao tao,Vec DXL,Vec DXU) 506 { 507 TAO_BQPIP *qp = (TAO_BQPIP*)tao->data; 508 509 PetscFunctionBegin; 510 PetscCall(VecCopy(qp->Z,DXL)); 511 PetscCall(VecCopy(qp->S,DXU)); 512 PetscCall(VecScale(DXU,-1.0)); 513 PetscFunctionReturn(0); 514 } 515 516 /*MC 517 TAOBQPIP - interior-point method for quadratic programs with 518 box constraints. 519 520 Notes: 521 This algorithms solves quadratic problems only, the Hessian will 522 only be computed once. 523 524 Options Database Keys: 525 . -tao_bqpip_predcorr - use a predictor/corrector method 526 527 Level: beginner 528 M*/ 529 530 PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao) 531 { 532 TAO_BQPIP *qp; 533 534 PetscFunctionBegin; 535 PetscCall(PetscNewLog(tao,&qp)); 536 537 tao->ops->setup = TaoSetup_BQPIP; 538 tao->ops->solve = TaoSolve_BQPIP; 539 tao->ops->view = TaoView_BQPIP; 540 tao->ops->setfromoptions = TaoSetFromOptions_BQPIP; 541 tao->ops->destroy = TaoDestroy_BQPIP; 542 tao->ops->computedual = TaoComputeDual_BQPIP; 543 544 /* Override default settings (unless already changed) */ 545 if (!tao->max_it_changed) tao->max_it=100; 546 if (!tao->max_funcs_changed) tao->max_funcs = 500; 547 #if defined(PETSC_USE_REAL_SINGLE) 548 if (!tao->catol_changed) tao->catol=1e-6; 549 #else 550 if (!tao->catol_changed) tao->catol=1e-12; 551 #endif 552 553 /* Initialize pointers and variables */ 554 qp->n = 0; 555 qp->m = 0; 556 557 qp->predcorr = 1; 558 tao->data = (void*)qp; 559 560 PetscCall(KSPCreate(((PetscObject)tao)->comm,&tao->ksp)); 561 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 562 PetscCall(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 563 PetscCall(KSPSetType(tao->ksp,KSPCG)); 564 PetscCall(KSPSetTolerances(tao->ksp,1e-14,1e-30,1e30,PetscMax(10,qp->n))); 565 PetscFunctionReturn(0); 566 } 567