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 = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 55 ierr = PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);CHKERRQ(ierr); 56 ierr = PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);CHKERRQ(ierr); 57 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 58 } 59 PetscFunctionReturn(0); 60 } 61 62 /* ---------------------------------------------------------- */ 63 static PetscErrorCode TaoSetup_TRON(Tao tao) 64 { 65 PetscErrorCode ierr; 66 TAO_TRON *tron = (TAO_TRON *)tao->data; 67 68 PetscFunctionBegin; 69 70 /* Allocate some arrays */ 71 ierr = VecDuplicate(tao->solution, &tron->diag);CHKERRQ(ierr); 72 ierr = VecDuplicate(tao->solution, &tron->X_New);CHKERRQ(ierr); 73 ierr = VecDuplicate(tao->solution, &tron->G_New);CHKERRQ(ierr); 74 ierr = VecDuplicate(tao->solution, &tron->Work);CHKERRQ(ierr); 75 ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr); 76 ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr); 77 if (!tao->XL) { 78 ierr = VecDuplicate(tao->solution, &tao->XL);CHKERRQ(ierr); 79 ierr = VecSet(tao->XL, PETSC_NINFINITY);CHKERRQ(ierr); 80 } 81 if (!tao->XU) { 82 ierr = VecDuplicate(tao->solution, &tao->XU);CHKERRQ(ierr); 83 ierr = VecSet(tao->XU, PETSC_INFINITY);CHKERRQ(ierr); 84 } 85 PetscFunctionReturn(0); 86 } 87 88 static PetscErrorCode TaoSolve_TRON(Tao tao) 89 { 90 TAO_TRON *tron = (TAO_TRON *)tao->data; 91 PetscErrorCode ierr; 92 PetscInt its; 93 TaoConvergedReason reason = TAO_CONTINUE_ITERATING; 94 TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING; 95 PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize; 96 97 PetscFunctionBegin; 98 tron->pgstepsize = 1.0; 99 tao->trust = tao->trust0; 100 /* Project the current point onto the feasible set */ 101 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 102 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 103 104 /* Project the initial point onto the feasible region */ 105 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 106 107 /* Compute the objective function and gradient */ 108 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);CHKERRQ(ierr); 109 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 110 if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 111 112 /* Project the gradient and calculate the norm */ 113 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 114 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 115 116 /* Initialize trust region radius */ 117 tao->trust=tao->trust0; 118 if (tao->trust <= 0) { 119 tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0); 120 } 121 122 /* Initialize step sizes for the line searches */ 123 tron->pgstepsize=1.0; 124 tron->stepsize=tao->trust; 125 126 ierr = TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,tron->stepsize,&reason);CHKERRQ(ierr); 127 while (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 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 141 ierr = VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);CHKERRQ(ierr); 142 ierr = ISGetSize(tron->Free_Local, &tron->n_free);CHKERRQ(ierr); 143 144 /* If no free variables */ 145 if (tron->n_free == 0) { 146 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 147 ierr = TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta, &reason);CHKERRQ(ierr); 148 if (!reason) { 149 reason = TAO_CONVERGED_STEPTOL; 150 ierr = TaoSetConvergedReason(tao,reason);CHKERRQ(ierr); 151 } 152 break; 153 } 154 /* use free_local to mask/submat gradient, hessian, stepdirection */ 155 ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);CHKERRQ(ierr); 156 ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);CHKERRQ(ierr); 157 ierr = VecSet(tron->DXFree,0.0);CHKERRQ(ierr); 158 ierr = VecScale(tron->R, -1.0);CHKERRQ(ierr); 159 ierr = TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);CHKERRQ(ierr); 160 if (tao->hessian == tao->hessian_pre) { 161 ierr = MatDestroy(&tron->Hpre_sub);CHKERRQ(ierr); 162 ierr = PetscObjectReference((PetscObject)(tron->H_sub));CHKERRQ(ierr); 163 tron->Hpre_sub = tron->H_sub; 164 } else { 165 ierr = TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);CHKERRQ(ierr); 166 } 167 ierr = KSPReset(tao->ksp);CHKERRQ(ierr); 168 ierr = KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);CHKERRQ(ierr); 169 while (1) { 170 171 /* Approximately solve the reduced linear system */ 172 ierr = KSPCGSetRadius(tao->ksp,delta);CHKERRQ(ierr); 173 174 ierr = KSPSolve(tao->ksp, tron->R, tron->DXFree);CHKERRQ(ierr); 175 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 176 tao->ksp_its+=its; 177 tao->ksp_tot_its+=its; 178 ierr = VecSet(tao->stepdirection,0.0);CHKERRQ(ierr); 179 180 /* Add dxfree matrix to compute step direction vector */ 181 ierr = VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);CHKERRQ(ierr); 182 183 ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr); 184 ierr = VecCopy(tao->solution, tron->X_New);CHKERRQ(ierr); 185 ierr = VecCopy(tao->gradient, tron->G_New);CHKERRQ(ierr); 186 187 stepsize=1.0;f_new=f; 188 189 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr); 190 ierr = TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);CHKERRQ(ierr);CHKERRQ(ierr); 191 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 192 193 ierr = MatMult(tao->hessian, tao->stepdirection, tron->Work);CHKERRQ(ierr); 194 ierr = VecAYPX(tron->Work, 0.5, tao->gradient);CHKERRQ(ierr); 195 ierr = VecDot(tao->stepdirection, tron->Work, &prered);CHKERRQ(ierr); 196 actred = f_new - f; 197 if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) { 198 rhok = 1.0; 199 } else if (actred<0) { 200 rhok=PetscAbs(-actred/prered); 201 } else { 202 rhok=0.0; 203 } 204 205 /* Compare actual improvement to the quadratic model */ 206 if (rhok > tron->eta1) { /* Accept the point */ 207 /* d = x_new - x */ 208 ierr = VecCopy(tron->X_New, tao->stepdirection);CHKERRQ(ierr); 209 ierr = VecAXPY(tao->stepdirection, -1.0, tao->solution);CHKERRQ(ierr); 210 211 ierr = VecNorm(tao->stepdirection, NORM_2, &xdiff);CHKERRQ(ierr); 212 xdiff *= stepsize; 213 214 /* Adjust trust region size */ 215 if (rhok < tron->eta2 ){ 216 delta = PetscMin(xdiff,delta)*tron->sigma1; 217 } else if (rhok > tron->eta4 ){ 218 delta= PetscMin(xdiff,delta)*tron->sigma3; 219 } else if (rhok > tron->eta3 ){ 220 delta=PetscMin(xdiff,delta)*tron->sigma2; 221 } 222 ierr = VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 223 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 224 ierr = VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);CHKERRQ(ierr); 225 f=f_new; 226 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 227 ierr = VecCopy(tron->X_New, tao->solution);CHKERRQ(ierr); 228 ierr = VecCopy(tron->G_New, tao->gradient);CHKERRQ(ierr); 229 break; 230 } 231 else if (delta <= 1e-30) { 232 break; 233 } 234 else { 235 delta /= 4.0; 236 } 237 } /* end linear solve loop */ 238 239 tron->f=f; tron->actred=actred; tao->trust=delta; 240 tao->niter++; 241 ierr = TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta, &reason);CHKERRQ(ierr); 242 } /* END MAIN LOOP */ 243 PetscFunctionReturn(0); 244 } 245 246 static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron) 247 { 248 PetscErrorCode ierr; 249 PetscInt i; 250 TaoLineSearchConvergedReason ls_reason; 251 PetscReal actred=-1.0,actred_max=0.0; 252 PetscReal f_new; 253 /* 254 The gradient and function value passed into and out of this 255 routine should be current and correct. 256 257 The free, active, and binding variables should be already identified 258 */ 259 PetscFunctionBegin; 260 261 for (i=0;i<tron->maxgpits;++i){ 262 263 if (-actred <= (tron->pg_ftol)*actred_max) break; 264 265 ++tron->gp_iterates; 266 ++tron->total_gp_its; 267 f_new=tron->f; 268 269 ierr = VecCopy(tao->gradient,tao->stepdirection);CHKERRQ(ierr); 270 ierr = VecScale(tao->stepdirection,-1.0);CHKERRQ(ierr); 271 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);CHKERRQ(ierr); 272 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, 273 &tron->pgstepsize, &ls_reason);CHKERRQ(ierr); 274 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 275 276 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 277 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 278 279 /* Update the iterate */ 280 actred = f_new - tron->f; 281 actred_max = PetscMax(actred_max,-(f_new - tron->f)); 282 tron->f = f_new; 283 } 284 PetscFunctionReturn(0); 285 } 286 287 static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) { 288 289 TAO_TRON *tron = (TAO_TRON *)tao->data; 290 PetscErrorCode ierr; 291 292 PetscFunctionBegin; 293 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 294 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 295 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 296 if (!tron->Work || !tao->gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n"); 297 298 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);CHKERRQ(ierr); 299 ierr = VecCopy(tron->Work,DXL);CHKERRQ(ierr); 300 ierr = VecAXPY(DXL,-1.0,tao->gradient);CHKERRQ(ierr); 301 ierr = VecSet(DXU,0.0);CHKERRQ(ierr); 302 ierr = VecPointwiseMax(DXL,DXL,DXU);CHKERRQ(ierr); 303 304 ierr = VecCopy(tao->gradient,DXU);CHKERRQ(ierr); 305 ierr = VecAXPY(DXU,-1.0,tron->Work);CHKERRQ(ierr); 306 ierr = VecSet(tron->Work,0.0);CHKERRQ(ierr); 307 ierr = VecPointwiseMin(DXU,tron->Work,DXU);CHKERRQ(ierr); 308 PetscFunctionReturn(0); 309 } 310 311 /*------------------------------------------------------------*/ 312 /*MC 313 TAOTRON - The TRON algorithm is an active-set Newton trust region method 314 for bound-constrained minimization. 315 316 Options Database Keys: 317 + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate 318 - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets 319 320 Level: beginner 321 M*/ 322 PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao) 323 { 324 TAO_TRON *tron; 325 PetscErrorCode ierr; 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 ierr = PetscNewLog(tao,&tron);CHKERRQ(ierr); 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 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);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 = KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);CHKERRQ(ierr); 379 ierr = KSPSetType(tao->ksp,KSPCGSTCG);CHKERRQ(ierr); 380 PetscFunctionReturn(0); 381 } 382