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 CHKERRQ(VecDestroy(&tron->X_New)); 13 CHKERRQ(VecDestroy(&tron->G_New)); 14 CHKERRQ(VecDestroy(&tron->Work)); 15 CHKERRQ(VecDestroy(&tron->DXFree)); 16 CHKERRQ(VecDestroy(&tron->R)); 17 CHKERRQ(VecDestroy(&tron->diag)); 18 CHKERRQ(VecScatterDestroy(&tron->scatter)); 19 CHKERRQ(ISDestroy(&tron->Free_Local)); 20 CHKERRQ(MatDestroy(&tron->H_sub)); 21 CHKERRQ(MatDestroy(&tron->Hpre_sub)); 22 CHKERRQ(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 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"Newton Trust Region Method for bound constrained optimization")); 34 CHKERRQ(PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg)); 35 CHKERRQ(PetscOptionsTail()); 36 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 48 if (isascii) { 49 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its)); 50 CHKERRQ(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 CHKERRQ(VecDuplicate(tao->solution, &tron->diag)); 64 CHKERRQ(VecDuplicate(tao->solution, &tron->X_New)); 65 CHKERRQ(VecDuplicate(tao->solution, &tron->G_New)); 66 CHKERRQ(VecDuplicate(tao->solution, &tron->Work)); 67 CHKERRQ(VecDuplicate(tao->solution, &tao->gradient)); 68 CHKERRQ(VecDuplicate(tao->solution, &tao->stepdirection)); 69 if (!tao->XL) { 70 CHKERRQ(VecDuplicate(tao->solution, &tao->XL)); 71 CHKERRQ(VecSet(tao->XL, PETSC_NINFINITY)); 72 } 73 if (!tao->XU) { 74 CHKERRQ(VecDuplicate(tao->solution, &tao->XU)); 75 CHKERRQ(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 CHKERRQ(TaoComputeVariableBounds(tao)); 92 CHKERRQ(TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU)); 93 94 /* Project the initial point onto the feasible region */ 95 CHKERRQ(VecMedian(tao->XL,tao->solution,tao->XU,tao->solution)); 96 97 /* Compute the objective function and gradient */ 98 CHKERRQ(TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient)); 99 CHKERRQ(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 CHKERRQ(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient)); 104 CHKERRQ(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 CHKERRQ(TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its)); 118 CHKERRQ(TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,tron->stepsize)); 119 CHKERRQ((*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 CHKERRQ((*tao->ops->update)(tao, tao->niter, tao->user_update)); 124 } 125 126 /* Perform projected gradient iterations */ 127 CHKERRQ(TronGradientProjections(tao,tron)); 128 129 CHKERRQ(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient)); 130 CHKERRQ(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 CHKERRQ(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 136 137 /* Generate index set (IS) of which bound constraints are active */ 138 CHKERRQ(ISDestroy(&tron->Free_Local)); 139 CHKERRQ(VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local)); 140 CHKERRQ(ISGetSize(tron->Free_Local, &tron->n_free)); 141 142 /* If no free variables */ 143 if (tron->n_free == 0) { 144 CHKERRQ(VecNorm(tao->gradient,NORM_2,&tron->gnorm)); 145 CHKERRQ(TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its)); 146 CHKERRQ(TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta)); 147 CHKERRQ((*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 CHKERRQ(TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R)); 155 CHKERRQ(TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree)); 156 CHKERRQ(VecSet(tron->DXFree,0.0)); 157 CHKERRQ(VecScale(tron->R, -1.0)); 158 CHKERRQ(TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub)); 159 if (tao->hessian == tao->hessian_pre) { 160 CHKERRQ(MatDestroy(&tron->Hpre_sub)); 161 CHKERRQ(PetscObjectReference((PetscObject)(tron->H_sub))); 162 tron->Hpre_sub = tron->H_sub; 163 } else { 164 CHKERRQ(TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub)); 165 } 166 CHKERRQ(KSPReset(tao->ksp)); 167 CHKERRQ(KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub)); 168 while (1) { 169 170 /* Approximately solve the reduced linear system */ 171 CHKERRQ(KSPCGSetRadius(tao->ksp,delta)); 172 173 CHKERRQ(KSPSolve(tao->ksp, tron->R, tron->DXFree)); 174 CHKERRQ(KSPGetIterationNumber(tao->ksp,&its)); 175 tao->ksp_its+=its; 176 tao->ksp_tot_its+=its; 177 CHKERRQ(VecSet(tao->stepdirection,0.0)); 178 179 /* Add dxfree matrix to compute step direction vector */ 180 CHKERRQ(VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree)); 181 182 CHKERRQ(VecDot(tao->gradient, tao->stepdirection, &gdx)); 183 CHKERRQ(VecCopy(tao->solution, tron->X_New)); 184 CHKERRQ(VecCopy(tao->gradient, tron->G_New)); 185 186 stepsize=1.0;f_new=f; 187 188 CHKERRQ(TaoLineSearchSetInitialStepLength(tao->linesearch,1.0)); 189 CHKERRQ(TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason)); 190 CHKERRQ(TaoAddLineSearchCounts(tao)); 191 192 CHKERRQ(MatMult(tao->hessian, tao->stepdirection, tron->Work)); 193 CHKERRQ(VecAYPX(tron->Work, 0.5, tao->gradient)); 194 CHKERRQ(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 CHKERRQ(VecCopy(tron->X_New, tao->stepdirection)); 208 CHKERRQ(VecAXPY(tao->stepdirection, -1.0, tao->solution)); 209 210 CHKERRQ(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 CHKERRQ(VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient)); 222 CHKERRQ(ISDestroy(&tron->Free_Local)); 223 CHKERRQ(VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local)); 224 f=f_new; 225 CHKERRQ(VecNorm(tao->gradient,NORM_2,&tron->gnorm)); 226 CHKERRQ(VecCopy(tron->X_New, tao->solution)); 227 CHKERRQ(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 CHKERRQ(TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its)); 241 CHKERRQ(TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,stepsize)); 242 CHKERRQ((*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 CHKERRQ(VecCopy(tao->gradient,tao->stepdirection)); 271 CHKERRQ(VecScale(tao->stepdirection,-1.0)); 272 CHKERRQ(TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize)); 273 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, 274 &tron->pgstepsize, &ls_reason);CHKERRQ(ierr); 275 CHKERRQ(TaoAddLineSearchCounts(tao)); 276 277 CHKERRQ(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient)); 278 CHKERRQ(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 CHKERRQ(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work)); 300 CHKERRQ(VecCopy(tron->Work,DXL)); 301 CHKERRQ(VecAXPY(DXL,-1.0,tao->gradient)); 302 CHKERRQ(VecSet(DXU,0.0)); 303 CHKERRQ(VecPointwiseMax(DXL,DXL,DXU)); 304 305 CHKERRQ(VecCopy(tao->gradient,DXU)); 306 CHKERRQ(VecAXPY(DXU,-1.0,tron->Work)); 307 CHKERRQ(VecSet(tron->Work,0.0)); 308 CHKERRQ(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 CHKERRQ(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 CHKERRQ(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 373 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 374 CHKERRQ(TaoLineSearchSetType(tao->linesearch,morethuente_type)); 375 CHKERRQ(TaoLineSearchUseTaoRoutines(tao->linesearch,tao)); 376 CHKERRQ(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 377 378 CHKERRQ(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 379 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 380 CHKERRQ(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix)); 381 CHKERRQ(KSPSetType(tao->ksp,KSPSTCG)); 382 PetscFunctionReturn(0); 383 } 384