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