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