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 129 /* Perform projected gradient iterations */ 130 ierr = TronGradientProjections(tao,tron);CHKERRQ(ierr); 131 132 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 133 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 134 135 tao->ksp_its=0; 136 f=tron->f; delta=tao->trust; 137 tron->n_free_last = tron->n_free; 138 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 139 140 /* Generate index set (IS) of which bound constraints are active */ 141 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 142 ierr = VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);CHKERRQ(ierr); 143 ierr = ISGetSize(tron->Free_Local, &tron->n_free);CHKERRQ(ierr); 144 145 /* If no free variables */ 146 if (tron->n_free == 0) { 147 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 148 ierr = TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 149 ierr = TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta);CHKERRQ(ierr); 150 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 151 if (!tao->reason) { 152 tao->reason = TAO_CONVERGED_STEPTOL; 153 } 154 break; 155 } 156 /* use free_local to mask/submat gradient, hessian, stepdirection */ 157 ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);CHKERRQ(ierr); 158 ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);CHKERRQ(ierr); 159 ierr = VecSet(tron->DXFree,0.0);CHKERRQ(ierr); 160 ierr = VecScale(tron->R, -1.0);CHKERRQ(ierr); 161 ierr = TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);CHKERRQ(ierr); 162 if (tao->hessian == tao->hessian_pre) { 163 ierr = MatDestroy(&tron->Hpre_sub);CHKERRQ(ierr); 164 ierr = PetscObjectReference((PetscObject)(tron->H_sub));CHKERRQ(ierr); 165 tron->Hpre_sub = tron->H_sub; 166 } else { 167 ierr = TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);CHKERRQ(ierr); 168 } 169 ierr = KSPReset(tao->ksp);CHKERRQ(ierr); 170 ierr = KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);CHKERRQ(ierr); 171 while (1) { 172 173 /* Approximately solve the reduced linear system */ 174 ierr = KSPCGSetRadius(tao->ksp,delta);CHKERRQ(ierr); 175 176 ierr = KSPSolve(tao->ksp, tron->R, tron->DXFree);CHKERRQ(ierr); 177 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 178 tao->ksp_its+=its; 179 tao->ksp_tot_its+=its; 180 ierr = VecSet(tao->stepdirection,0.0);CHKERRQ(ierr); 181 182 /* Add dxfree matrix to compute step direction vector */ 183 ierr = VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);CHKERRQ(ierr); 184 185 ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr); 186 ierr = VecCopy(tao->solution, tron->X_New);CHKERRQ(ierr); 187 ierr = VecCopy(tao->gradient, tron->G_New);CHKERRQ(ierr); 188 189 stepsize=1.0;f_new=f; 190 191 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr); 192 ierr = TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);CHKERRQ(ierr);CHKERRQ(ierr); 193 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 194 195 ierr = MatMult(tao->hessian, tao->stepdirection, tron->Work);CHKERRQ(ierr); 196 ierr = VecAYPX(tron->Work, 0.5, tao->gradient);CHKERRQ(ierr); 197 ierr = VecDot(tao->stepdirection, tron->Work, &prered);CHKERRQ(ierr); 198 actred = f_new - f; 199 if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) { 200 rhok = 1.0; 201 } else if (actred<0) { 202 rhok=PetscAbs(-actred/prered); 203 } else { 204 rhok=0.0; 205 } 206 207 /* Compare actual improvement to the quadratic model */ 208 if (rhok > tron->eta1) { /* Accept the point */ 209 /* d = x_new - x */ 210 ierr = VecCopy(tron->X_New, tao->stepdirection);CHKERRQ(ierr); 211 ierr = VecAXPY(tao->stepdirection, -1.0, tao->solution);CHKERRQ(ierr); 212 213 ierr = VecNorm(tao->stepdirection, NORM_2, &xdiff);CHKERRQ(ierr); 214 xdiff *= stepsize; 215 216 /* Adjust trust region size */ 217 if (rhok < tron->eta2 ){ 218 delta = PetscMin(xdiff,delta)*tron->sigma1; 219 } else if (rhok > tron->eta4 ){ 220 delta= PetscMin(xdiff,delta)*tron->sigma3; 221 } else if (rhok > tron->eta3 ){ 222 delta=PetscMin(xdiff,delta)*tron->sigma2; 223 } 224 ierr = VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 225 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 226 ierr = VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);CHKERRQ(ierr); 227 f=f_new; 228 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 229 ierr = VecCopy(tron->X_New, tao->solution);CHKERRQ(ierr); 230 ierr = VecCopy(tron->G_New, tao->gradient);CHKERRQ(ierr); 231 break; 232 } 233 else if (delta <= 1e-30) { 234 break; 235 } 236 else { 237 delta /= 4.0; 238 } 239 } /* end linear solve loop */ 240 241 tron->f=f; tron->actred=actred; tao->trust=delta; 242 tao->niter++; 243 ierr = TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 244 ierr = TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta);CHKERRQ(ierr); 245 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 246 } /* END MAIN LOOP */ 247 PetscFunctionReturn(0); 248 } 249 250 static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron) 251 { 252 PetscErrorCode ierr; 253 PetscInt i; 254 TaoLineSearchConvergedReason ls_reason; 255 PetscReal actred=-1.0,actred_max=0.0; 256 PetscReal f_new; 257 /* 258 The gradient and function value passed into and out of this 259 routine should be current and correct. 260 261 The free, active, and binding variables should be already identified 262 */ 263 PetscFunctionBegin; 264 265 for (i=0;i<tron->maxgpits;++i){ 266 267 if (-actred <= (tron->pg_ftol)*actred_max) break; 268 269 ++tron->gp_iterates; 270 ++tron->total_gp_its; 271 f_new=tron->f; 272 273 ierr = VecCopy(tao->gradient,tao->stepdirection);CHKERRQ(ierr); 274 ierr = VecScale(tao->stepdirection,-1.0);CHKERRQ(ierr); 275 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);CHKERRQ(ierr); 276 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, 277 &tron->pgstepsize, &ls_reason);CHKERRQ(ierr); 278 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 279 280 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 281 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 282 283 /* Update the iterate */ 284 actred = f_new - tron->f; 285 actred_max = PetscMax(actred_max,-(f_new - tron->f)); 286 tron->f = f_new; 287 } 288 PetscFunctionReturn(0); 289 } 290 291 static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) { 292 293 TAO_TRON *tron = (TAO_TRON *)tao->data; 294 PetscErrorCode ierr; 295 296 PetscFunctionBegin; 297 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 298 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 299 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 300 if (!tron->Work || !tao->gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n"); 301 302 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);CHKERRQ(ierr); 303 ierr = VecCopy(tron->Work,DXL);CHKERRQ(ierr); 304 ierr = VecAXPY(DXL,-1.0,tao->gradient);CHKERRQ(ierr); 305 ierr = VecSet(DXU,0.0);CHKERRQ(ierr); 306 ierr = VecPointwiseMax(DXL,DXL,DXU);CHKERRQ(ierr); 307 308 ierr = VecCopy(tao->gradient,DXU);CHKERRQ(ierr); 309 ierr = VecAXPY(DXU,-1.0,tron->Work);CHKERRQ(ierr); 310 ierr = VecSet(tron->Work,0.0);CHKERRQ(ierr); 311 ierr = VecPointwiseMin(DXU,tron->Work,DXU);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 /*------------------------------------------------------------*/ 316 /*MC 317 TAOTRON - The TRON algorithm is an active-set Newton trust region method 318 for bound-constrained minimization. 319 320 Options Database Keys: 321 + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate 322 - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets 323 324 Level: beginner 325 M*/ 326 PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao) 327 { 328 TAO_TRON *tron; 329 PetscErrorCode ierr; 330 const char *morethuente_type = TAOLINESEARCHMT; 331 332 PetscFunctionBegin; 333 tao->ops->setup = TaoSetup_TRON; 334 tao->ops->solve = TaoSolve_TRON; 335 tao->ops->view = TaoView_TRON; 336 tao->ops->setfromoptions = TaoSetFromOptions_TRON; 337 tao->ops->destroy = TaoDestroy_TRON; 338 tao->ops->computedual = TaoComputeDual_TRON; 339 340 ierr = PetscNewLog(tao,&tron);CHKERRQ(ierr); 341 tao->data = (void*)tron; 342 343 /* Override default settings (unless already changed) */ 344 if (!tao->max_it_changed) tao->max_it = 50; 345 if (!tao->trust0_changed) tao->trust0 = 1.0; 346 if (!tao->steptol_changed) tao->steptol = 0.0; 347 348 /* Initialize pointers and variables */ 349 tron->n = 0; 350 tron->maxgpits = 3; 351 tron->pg_ftol = 0.001; 352 353 tron->eta1 = 1.0e-4; 354 tron->eta2 = 0.25; 355 tron->eta3 = 0.50; 356 tron->eta4 = 0.90; 357 358 tron->sigma1 = 0.5; 359 tron->sigma2 = 2.0; 360 tron->sigma3 = 4.0; 361 362 tron->gp_iterates = 0; /* Cumulative number */ 363 tron->total_gp_its = 0; 364 tron->n_free = 0; 365 366 tron->DXFree=NULL; 367 tron->R=NULL; 368 tron->X_New=NULL; 369 tron->G_New=NULL; 370 tron->Work=NULL; 371 tron->Free_Local=NULL; 372 tron->H_sub=NULL; 373 tron->Hpre_sub=NULL; 374 tao->subset_type = TAO_SUBSET_SUBVEC; 375 376 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr); 377 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr); 378 ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr); 379 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr); 380 ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr); 381 382 ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr); 383 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);CHKERRQ(ierr); 384 ierr = KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);CHKERRQ(ierr); 385 ierr = KSPSetType(tao->ksp,KSPCGSTCG);CHKERRQ(ierr); 386 PetscFunctionReturn(0); 387 } 388