1 #include <../src/tao/bound/impls/tron/tron.h> 2 #include <../src/tao/matrix/submatfree.h> 3 4 /* TRON Routines */ 5 static PetscErrorCode TronGradientProjections(Tao,TAO_TRON*); 6 /*------------------------------------------------------------*/ 7 static PetscErrorCode TaoDestroy_TRON(Tao tao) 8 { 9 TAO_TRON *tron = (TAO_TRON *)tao->data; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 ierr = VecDestroy(&tron->X_New);CHKERRQ(ierr); 14 ierr = VecDestroy(&tron->G_New);CHKERRQ(ierr); 15 ierr = VecDestroy(&tron->Work);CHKERRQ(ierr); 16 ierr = VecDestroy(&tron->DXFree);CHKERRQ(ierr); 17 ierr = VecDestroy(&tron->R);CHKERRQ(ierr); 18 ierr = VecDestroy(&tron->diag);CHKERRQ(ierr); 19 ierr = VecScatterDestroy(&tron->scatter);CHKERRQ(ierr); 20 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 21 ierr = MatDestroy(&tron->H_sub);CHKERRQ(ierr); 22 ierr = MatDestroy(&tron->Hpre_sub);CHKERRQ(ierr); 23 ierr = PetscFree(tao->data);CHKERRQ(ierr); 24 PetscFunctionReturn(0); 25 } 26 27 /*------------------------------------------------------------*/ 28 static PetscErrorCode TaoSetFromOptions_TRON(PetscOptionItems *PetscOptionsObject,Tao tao) 29 { 30 TAO_TRON *tron = (TAO_TRON *)tao->data; 31 PetscErrorCode ierr; 32 PetscBool flg; 33 34 PetscFunctionBegin; 35 ierr = PetscOptionsHead(PetscOptionsObject,"Newton Trust Region Method for bound constrained optimization");CHKERRQ(ierr); 36 ierr = PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);CHKERRQ(ierr); 37 ierr = PetscOptionsTail();CHKERRQ(ierr); 38 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 39 PetscFunctionReturn(0); 40 } 41 42 /*------------------------------------------------------------*/ 43 static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer) 44 { 45 TAO_TRON *tron = (TAO_TRON *)tao->data; 46 PetscBool isascii; 47 PetscErrorCode ierr; 48 49 PetscFunctionBegin; 50 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 51 if (isascii) { 52 ierr = PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);CHKERRQ(ierr); 53 ierr = PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);CHKERRQ(ierr); 54 } 55 PetscFunctionReturn(0); 56 } 57 58 /* ---------------------------------------------------------- */ 59 static PetscErrorCode TaoSetup_TRON(Tao tao) 60 { 61 PetscErrorCode ierr; 62 TAO_TRON *tron = (TAO_TRON *)tao->data; 63 64 PetscFunctionBegin; 65 66 /* Allocate some arrays */ 67 ierr = VecDuplicate(tao->solution, &tron->diag);CHKERRQ(ierr); 68 ierr = VecDuplicate(tao->solution, &tron->X_New);CHKERRQ(ierr); 69 ierr = VecDuplicate(tao->solution, &tron->G_New);CHKERRQ(ierr); 70 ierr = VecDuplicate(tao->solution, &tron->Work);CHKERRQ(ierr); 71 ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr); 72 ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr); 73 if (!tao->XL) { 74 ierr = VecDuplicate(tao->solution, &tao->XL);CHKERRQ(ierr); 75 ierr = VecSet(tao->XL, PETSC_NINFINITY);CHKERRQ(ierr); 76 } 77 if (!tao->XU) { 78 ierr = VecDuplicate(tao->solution, &tao->XU);CHKERRQ(ierr); 79 ierr = VecSet(tao->XU, PETSC_INFINITY);CHKERRQ(ierr); 80 } 81 PetscFunctionReturn(0); 82 } 83 84 static PetscErrorCode TaoSolve_TRON(Tao tao) 85 { 86 TAO_TRON *tron = (TAO_TRON *)tao->data; 87 PetscErrorCode ierr; 88 PetscInt its; 89 TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING; 90 PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize; 91 92 PetscFunctionBegin; 93 tron->pgstepsize = 1.0; 94 tao->trust = tao->trust0; 95 /* Project the current point onto the feasible set */ 96 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 97 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 98 99 /* Project the initial point onto the feasible region */ 100 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 101 102 /* Compute the objective function and gradient */ 103 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);CHKERRQ(ierr); 104 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 105 PetscCheck(!PetscIsInfOrNanReal(tron->f) && !PetscIsInfOrNanReal(tron->gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 106 107 /* Project the gradient and calculate the norm */ 108 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 109 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 110 111 /* Initialize trust region radius */ 112 tao->trust=tao->trust0; 113 if (tao->trust <= 0) { 114 tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0); 115 } 116 117 /* Initialize step sizes for the line searches */ 118 tron->pgstepsize=1.0; 119 tron->stepsize=tao->trust; 120 121 tao->reason = TAO_CONTINUE_ITERATING; 122 ierr = TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 123 ierr = TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,tron->stepsize);CHKERRQ(ierr); 124 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 125 while (tao->reason==TAO_CONTINUE_ITERATING) { 126 /* Call general purpose update function */ 127 if (tao->ops->update) { 128 ierr = (*tao->ops->update)(tao, tao->niter, tao->user_update);CHKERRQ(ierr); 129 } 130 131 /* Perform projected gradient iterations */ 132 ierr = TronGradientProjections(tao,tron);CHKERRQ(ierr); 133 134 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 135 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 136 137 tao->ksp_its=0; 138 f=tron->f; delta=tao->trust; 139 tron->n_free_last = tron->n_free; 140 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 141 142 /* Generate index set (IS) of which bound constraints are active */ 143 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 144 ierr = VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);CHKERRQ(ierr); 145 ierr = ISGetSize(tron->Free_Local, &tron->n_free);CHKERRQ(ierr); 146 147 /* If no free variables */ 148 if (tron->n_free == 0) { 149 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 150 ierr = TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 151 ierr = TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta);CHKERRQ(ierr); 152 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 153 if (!tao->reason) { 154 tao->reason = TAO_CONVERGED_STEPTOL; 155 } 156 break; 157 } 158 /* use free_local to mask/submat gradient, hessian, stepdirection */ 159 ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);CHKERRQ(ierr); 160 ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);CHKERRQ(ierr); 161 ierr = VecSet(tron->DXFree,0.0);CHKERRQ(ierr); 162 ierr = VecScale(tron->R, -1.0);CHKERRQ(ierr); 163 ierr = TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);CHKERRQ(ierr); 164 if (tao->hessian == tao->hessian_pre) { 165 ierr = MatDestroy(&tron->Hpre_sub);CHKERRQ(ierr); 166 ierr = PetscObjectReference((PetscObject)(tron->H_sub));CHKERRQ(ierr); 167 tron->Hpre_sub = tron->H_sub; 168 } else { 169 ierr = TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);CHKERRQ(ierr); 170 } 171 ierr = KSPReset(tao->ksp);CHKERRQ(ierr); 172 ierr = KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);CHKERRQ(ierr); 173 while (1) { 174 175 /* Approximately solve the reduced linear system */ 176 ierr = KSPCGSetRadius(tao->ksp,delta);CHKERRQ(ierr); 177 178 ierr = KSPSolve(tao->ksp, tron->R, tron->DXFree);CHKERRQ(ierr); 179 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 180 tao->ksp_its+=its; 181 tao->ksp_tot_its+=its; 182 ierr = VecSet(tao->stepdirection,0.0);CHKERRQ(ierr); 183 184 /* Add dxfree matrix to compute step direction vector */ 185 ierr = VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);CHKERRQ(ierr); 186 187 ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr); 188 ierr = VecCopy(tao->solution, tron->X_New);CHKERRQ(ierr); 189 ierr = VecCopy(tao->gradient, tron->G_New);CHKERRQ(ierr); 190 191 stepsize=1.0;f_new=f; 192 193 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr); 194 ierr = TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);CHKERRQ(ierr); 195 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 196 197 ierr = MatMult(tao->hessian, tao->stepdirection, tron->Work);CHKERRQ(ierr); 198 ierr = VecAYPX(tron->Work, 0.5, tao->gradient);CHKERRQ(ierr); 199 ierr = VecDot(tao->stepdirection, tron->Work, &prered);CHKERRQ(ierr); 200 actred = f_new - f; 201 if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) { 202 rhok = 1.0; 203 } else if (actred<0) { 204 rhok=PetscAbs(-actred/prered); 205 } else { 206 rhok=0.0; 207 } 208 209 /* Compare actual improvement to the quadratic model */ 210 if (rhok > tron->eta1) { /* Accept the point */ 211 /* d = x_new - x */ 212 ierr = VecCopy(tron->X_New, tao->stepdirection);CHKERRQ(ierr); 213 ierr = VecAXPY(tao->stepdirection, -1.0, tao->solution);CHKERRQ(ierr); 214 215 ierr = VecNorm(tao->stepdirection, NORM_2, &xdiff);CHKERRQ(ierr); 216 xdiff *= stepsize; 217 218 /* Adjust trust region size */ 219 if (rhok < tron->eta2) { 220 delta = PetscMin(xdiff,delta)*tron->sigma1; 221 } else if (rhok > tron->eta4) { 222 delta= PetscMin(xdiff,delta)*tron->sigma3; 223 } else if (rhok > tron->eta3) { 224 delta=PetscMin(xdiff,delta)*tron->sigma2; 225 } 226 ierr = VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 227 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 228 ierr = VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);CHKERRQ(ierr); 229 f=f_new; 230 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 231 ierr = VecCopy(tron->X_New, tao->solution);CHKERRQ(ierr); 232 ierr = VecCopy(tron->G_New, tao->gradient);CHKERRQ(ierr); 233 break; 234 } 235 else if (delta <= 1e-30) { 236 break; 237 } 238 else { 239 delta /= 4.0; 240 } 241 } /* end linear solve loop */ 242 243 tron->f=f; tron->actred=actred; tao->trust=delta; 244 tao->niter++; 245 ierr = TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 246 ierr = TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,stepsize);CHKERRQ(ierr); 247 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 248 } /* END MAIN LOOP */ 249 PetscFunctionReturn(0); 250 } 251 252 static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron) 253 { 254 PetscErrorCode ierr; 255 PetscInt i; 256 TaoLineSearchConvergedReason ls_reason; 257 PetscReal actred=-1.0,actred_max=0.0; 258 PetscReal f_new; 259 /* 260 The gradient and function value passed into and out of this 261 routine should be current and correct. 262 263 The free, active, and binding variables should be already identified 264 */ 265 PetscFunctionBegin; 266 267 for (i=0;i<tron->maxgpits;++i) { 268 269 if (-actred <= (tron->pg_ftol)*actred_max) break; 270 271 ++tron->gp_iterates; 272 ++tron->total_gp_its; 273 f_new=tron->f; 274 275 ierr = VecCopy(tao->gradient,tao->stepdirection);CHKERRQ(ierr); 276 ierr = VecScale(tao->stepdirection,-1.0);CHKERRQ(ierr); 277 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);CHKERRQ(ierr); 278 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, 279 &tron->pgstepsize, &ls_reason);CHKERRQ(ierr); 280 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 281 282 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 283 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 284 285 /* Update the iterate */ 286 actred = f_new - tron->f; 287 actred_max = PetscMax(actred_max,-(f_new - tron->f)); 288 tron->f = f_new; 289 } 290 PetscFunctionReturn(0); 291 } 292 293 static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) 294 { 295 296 TAO_TRON *tron = (TAO_TRON *)tao->data; 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 301 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 302 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 303 PetscCheck(tron->Work && tao->gradient,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist."); 304 305 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);CHKERRQ(ierr); 306 ierr = VecCopy(tron->Work,DXL);CHKERRQ(ierr); 307 ierr = VecAXPY(DXL,-1.0,tao->gradient);CHKERRQ(ierr); 308 ierr = VecSet(DXU,0.0);CHKERRQ(ierr); 309 ierr = VecPointwiseMax(DXL,DXL,DXU);CHKERRQ(ierr); 310 311 ierr = VecCopy(tao->gradient,DXU);CHKERRQ(ierr); 312 ierr = VecAXPY(DXU,-1.0,tron->Work);CHKERRQ(ierr); 313 ierr = VecSet(tron->Work,0.0);CHKERRQ(ierr); 314 ierr = VecPointwiseMin(DXU,tron->Work,DXU);CHKERRQ(ierr); 315 PetscFunctionReturn(0); 316 } 317 318 /*------------------------------------------------------------*/ 319 /*MC 320 TAOTRON - The TRON algorithm is an active-set Newton trust region method 321 for bound-constrained minimization. 322 323 Options Database Keys: 324 + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate 325 - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets 326 327 Level: beginner 328 M*/ 329 PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao) 330 { 331 TAO_TRON *tron; 332 PetscErrorCode ierr; 333 const char *morethuente_type = TAOLINESEARCHMT; 334 335 PetscFunctionBegin; 336 tao->ops->setup = TaoSetup_TRON; 337 tao->ops->solve = TaoSolve_TRON; 338 tao->ops->view = TaoView_TRON; 339 tao->ops->setfromoptions = TaoSetFromOptions_TRON; 340 tao->ops->destroy = TaoDestroy_TRON; 341 tao->ops->computedual = TaoComputeDual_TRON; 342 343 ierr = PetscNewLog(tao,&tron);CHKERRQ(ierr); 344 tao->data = (void*)tron; 345 346 /* Override default settings (unless already changed) */ 347 if (!tao->max_it_changed) tao->max_it = 50; 348 if (!tao->trust0_changed) tao->trust0 = 1.0; 349 if (!tao->steptol_changed) tao->steptol = 0.0; 350 351 /* Initialize pointers and variables */ 352 tron->n = 0; 353 tron->maxgpits = 3; 354 tron->pg_ftol = 0.001; 355 356 tron->eta1 = 1.0e-4; 357 tron->eta2 = 0.25; 358 tron->eta3 = 0.50; 359 tron->eta4 = 0.90; 360 361 tron->sigma1 = 0.5; 362 tron->sigma2 = 2.0; 363 tron->sigma3 = 4.0; 364 365 tron->gp_iterates = 0; /* Cumulative number */ 366 tron->total_gp_its = 0; 367 tron->n_free = 0; 368 369 tron->DXFree=NULL; 370 tron->R=NULL; 371 tron->X_New=NULL; 372 tron->G_New=NULL; 373 tron->Work=NULL; 374 tron->Free_Local=NULL; 375 tron->H_sub=NULL; 376 tron->Hpre_sub=NULL; 377 tao->subset_type = TAO_SUBSET_SUBVEC; 378 379 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr); 380 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr); 381 ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr); 382 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr); 383 ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr); 384 385 ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr); 386 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);CHKERRQ(ierr); 387 ierr = KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);CHKERRQ(ierr); 388 ierr = KSPSetType(tao->ksp,KSPSTCG);CHKERRQ(ierr); 389 PetscFunctionReturn(0); 390 } 391