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