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