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