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