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