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