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