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