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 if (tron->Free_Local) { 214 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 215 } 216 ierr = VecWhichBetween(tao->XL, tron->X_New, tao->XU, &tron->Free_Local);CHKERRQ(ierr); 217 f=f_new; 218 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 219 ierr = VecCopy(tron->X_New, tao->solution);CHKERRQ(ierr); 220 ierr = VecCopy(tron->G_New, tao->gradient);CHKERRQ(ierr); 221 break; 222 } 223 else if (delta <= 1e-30) { 224 break; 225 } 226 else { 227 delta /= 4.0; 228 } 229 } /* end linear solve loop */ 230 231 tron->f=f; tron->actred=actred; tao->trust=delta; 232 tao->niter++; 233 ierr = TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta, &reason);CHKERRQ(ierr); 234 } /* END MAIN LOOP */ 235 236 PetscFunctionReturn(0); 237 } 238 239 240 static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron) 241 { 242 PetscErrorCode ierr; 243 PetscInt i; 244 TaoLineSearchConvergedReason ls_reason; 245 PetscReal actred=-1.0,actred_max=0.0; 246 PetscReal f_new; 247 /* 248 The gradient and function value passed into and out of this 249 routine should be current and correct. 250 251 The free, active, and binding variables should be already identified 252 */ 253 PetscFunctionBegin; 254 if (tron->Free_Local) { 255 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 256 } 257 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 258 259 for (i=0;i<tron->maxgpits;i++){ 260 261 if ( -actred <= (tron->pg_ftol)*actred_max) break; 262 263 tron->gp_iterates++; tron->total_gp_its++; 264 f_new=tron->f; 265 266 ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr); 267 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 268 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);CHKERRQ(ierr); 269 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, 270 &tron->pgstepsize, &ls_reason);CHKERRQ(ierr); 271 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 272 273 274 /* Update the iterate */ 275 actred = f_new - tron->f; 276 actred_max = PetscMax(actred_max,-(f_new - tron->f)); 277 tron->f = f_new; 278 if (tron->Free_Local) { 279 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 280 } 281 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 282 } 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 383