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 #undef __FUNCT__ 9 #define __FUNCT__ "TaoDestroy_TRON" 10 static PetscErrorCode TaoDestroy_TRON(Tao tao) 11 { 12 TAO_TRON *tron = (TAO_TRON *)tao->data; 13 PetscErrorCode ierr; 14 15 PetscFunctionBegin; 16 ierr = VecDestroy(&tron->X_New);CHKERRQ(ierr); 17 ierr = VecDestroy(&tron->G_New);CHKERRQ(ierr); 18 ierr = VecDestroy(&tron->Work);CHKERRQ(ierr); 19 ierr = VecDestroy(&tron->DXFree);CHKERRQ(ierr); 20 ierr = VecDestroy(&tron->R);CHKERRQ(ierr); 21 ierr = VecDestroy(&tron->diag);CHKERRQ(ierr); 22 ierr = VecScatterDestroy(&tron->scatter);CHKERRQ(ierr); 23 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 24 ierr = MatDestroy(&tron->H_sub);CHKERRQ(ierr); 25 ierr = MatDestroy(&tron->Hpre_sub);CHKERRQ(ierr); 26 ierr = PetscFree(tao->data);CHKERRQ(ierr); 27 PetscFunctionReturn(0); 28 } 29 30 /*------------------------------------------------------------*/ 31 #undef __FUNCT__ 32 #define __FUNCT__ "TaoSetFromOptions_TRON" 33 static PetscErrorCode TaoSetFromOptions_TRON(PetscOptionItems *PetscOptionsObject,Tao tao) 34 { 35 TAO_TRON *tron = (TAO_TRON *)tao->data; 36 PetscErrorCode ierr; 37 PetscBool flg; 38 39 PetscFunctionBegin; 40 ierr = PetscOptionsHead(PetscOptionsObject,"Newton Trust Region Method for bound constrained optimization");CHKERRQ(ierr); 41 ierr = PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);CHKERRQ(ierr); 42 ierr = PetscOptionsTail();CHKERRQ(ierr); 43 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 44 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 /*------------------------------------------------------------*/ 49 #undef __FUNCT__ 50 #define __FUNCT__ "TaoView_TRON" 51 static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer) 52 { 53 TAO_TRON *tron = (TAO_TRON *)tao->data; 54 PetscBool isascii; 55 PetscErrorCode ierr; 56 57 PetscFunctionBegin; 58 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 59 if (isascii) { 60 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 61 ierr = PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);CHKERRQ(ierr); 62 ierr = PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);CHKERRQ(ierr); 63 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 64 } 65 PetscFunctionReturn(0); 66 } 67 68 69 /* ---------------------------------------------------------- */ 70 #undef __FUNCT__ 71 #define __FUNCT__ "TaoSetup_TRON" 72 static PetscErrorCode TaoSetup_TRON(Tao tao) 73 { 74 PetscErrorCode ierr; 75 TAO_TRON *tron = (TAO_TRON *)tao->data; 76 77 PetscFunctionBegin; 78 79 /* Allocate some arrays */ 80 ierr = VecDuplicate(tao->solution, &tron->diag);CHKERRQ(ierr); 81 ierr = VecDuplicate(tao->solution, &tron->X_New);CHKERRQ(ierr); 82 ierr = VecDuplicate(tao->solution, &tron->G_New);CHKERRQ(ierr); 83 ierr = VecDuplicate(tao->solution, &tron->Work);CHKERRQ(ierr); 84 ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr); 85 ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr); 86 if (!tao->XL) { 87 ierr = VecDuplicate(tao->solution, &tao->XL);CHKERRQ(ierr); 88 ierr = VecSet(tao->XL, PETSC_NINFINITY);CHKERRQ(ierr); 89 } 90 if (!tao->XU) { 91 ierr = VecDuplicate(tao->solution, &tao->XU);CHKERRQ(ierr); 92 ierr = VecSet(tao->XU, PETSC_INFINITY);CHKERRQ(ierr); 93 } 94 PetscFunctionReturn(0); 95 } 96 97 98 99 #undef __FUNCT__ 100 #define __FUNCT__ "TaoSolve_TRON" 101 static PetscErrorCode TaoSolve_TRON(Tao tao) 102 { 103 TAO_TRON *tron = (TAO_TRON *)tao->data; 104 PetscErrorCode ierr; 105 PetscInt its; 106 TaoConvergedReason reason = TAO_CONTINUE_ITERATING; 107 TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING; 108 PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize; 109 110 PetscFunctionBegin; 111 tron->pgstepsize=1.0; 112 tao->trust = tao->trust0; 113 /* Project the current point onto the feasible set */ 114 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 115 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 116 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 117 118 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);CHKERRQ(ierr); 119 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 120 121 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 122 123 /* Project the gradient and calculate the norm */ 124 ierr = VecBoundGradientProjection(tao->gradient,tao->solution, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 125 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 126 127 if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN"); 128 if (tao->trust <= 0) { 129 tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0); 130 } 131 132 tron->stepsize=tao->trust; 133 ierr = TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, tron->stepsize, &reason);CHKERRQ(ierr); 134 while (reason==TAO_CONTINUE_ITERATING){ 135 tao->ksp_its=0; 136 ierr = TronGradientProjections(tao,tron);CHKERRQ(ierr); 137 f=tron->f; delta=tao->trust; 138 tron->n_free_last = tron->n_free; 139 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 140 141 ierr = ISGetSize(tron->Free_Local, &tron->n_free);CHKERRQ(ierr); 142 143 /* If no free variables */ 144 if (tron->n_free == 0) { 145 ierr = PetscInfo(tao,"No free variables in tron iteration.\n");CHKERRQ(ierr); 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 = PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",(double)gdx);CHKERRQ(ierr); 185 186 ierr = VecCopy(tao->solution, tron->X_New);CHKERRQ(ierr); 187 ierr = VecCopy(tao->gradient, tron->G_New);CHKERRQ(ierr); 188 189 stepsize=1.0;f_new=f; 190 191 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr); 192 ierr = TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);CHKERRQ(ierr);CHKERRQ(ierr); 193 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 194 195 ierr = MatMult(tao->hessian, tao->stepdirection, tron->Work);CHKERRQ(ierr); 196 ierr = VecAYPX(tron->Work, 0.5, tao->gradient);CHKERRQ(ierr); 197 ierr = VecDot(tao->stepdirection, tron->Work, &prered);CHKERRQ(ierr); 198 actred = f_new - f; 199 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 if (tron->Free_Local) { 224 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 225 } 226 ierr = VecWhichBetween(tao->XL, tron->X_New, tao->XU, &tron->Free_Local);CHKERRQ(ierr); 227 f=f_new; 228 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 229 ierr = VecCopy(tron->X_New, tao->solution);CHKERRQ(ierr); 230 ierr = VecCopy(tron->G_New, tao->gradient);CHKERRQ(ierr); 231 break; 232 } 233 else if (delta <= 1e-30) { 234 break; 235 } 236 else { 237 delta /= 4.0; 238 } 239 } /* end linear solve loop */ 240 241 tron->f=f; tron->actred=actred; tao->trust=delta; 242 tao->niter++; 243 ierr = TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta, &reason);CHKERRQ(ierr); 244 } /* END MAIN LOOP */ 245 246 PetscFunctionReturn(0); 247 } 248 249 250 #undef __FUNCT__ 251 #define __FUNCT__ "TronGradientProjections" 252 static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron) 253 { 254 PetscErrorCode ierr; 255 PetscInt i; 256 TaoLineSearchConvergedReason ls_reason; 257 PetscReal actred=-1.0,actred_max=0.0; 258 PetscReal f_new; 259 /* 260 The gradient and function value passed into and out of this 261 routine should be current and correct. 262 263 The free, active, and binding variables should be already identified 264 */ 265 PetscFunctionBegin; 266 if (tron->Free_Local) { 267 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 268 } 269 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 270 271 for (i=0;i<tron->maxgpits;i++){ 272 273 if ( -actred <= (tron->pg_ftol)*actred_max) break; 274 275 tron->gp_iterates++; tron->total_gp_its++; 276 f_new=tron->f; 277 278 ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr); 279 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 280 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);CHKERRQ(ierr); 281 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, 282 &tron->pgstepsize, &ls_reason);CHKERRQ(ierr); 283 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 284 285 286 /* Update the iterate */ 287 actred = f_new - tron->f; 288 actred_max = PetscMax(actred_max,-(f_new - tron->f)); 289 tron->f = f_new; 290 if (tron->Free_Local) { 291 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 292 } 293 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 294 } 295 296 PetscFunctionReturn(0); 297 } 298 299 #undef __FUNCT__ 300 #define __FUNCT__ "TaoComputeDual_TRON" 301 static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) { 302 303 TAO_TRON *tron = (TAO_TRON *)tao->data; 304 PetscErrorCode ierr; 305 306 PetscFunctionBegin; 307 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 308 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 309 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 310 if (!tron->Work || !tao->gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n"); 311 312 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);CHKERRQ(ierr); 313 ierr = VecCopy(tron->Work,DXL);CHKERRQ(ierr); 314 ierr = VecAXPY(DXL,-1.0,tao->gradient);CHKERRQ(ierr); 315 ierr = VecSet(DXU,0.0);CHKERRQ(ierr); 316 ierr = VecPointwiseMax(DXL,DXL,DXU);CHKERRQ(ierr); 317 318 ierr = VecCopy(tao->gradient,DXU);CHKERRQ(ierr); 319 ierr = VecAXPY(DXU,-1.0,tron->Work);CHKERRQ(ierr); 320 ierr = VecSet(tron->Work,0.0);CHKERRQ(ierr); 321 ierr = VecPointwiseMin(DXU,tron->Work,DXU);CHKERRQ(ierr); 322 PetscFunctionReturn(0); 323 } 324 325 /*------------------------------------------------------------*/ 326 /*MC 327 TAOTRON - The TRON algorithm is an active-set Newton trust region method 328 for bound-constrained minimization. 329 330 Options Database Keys: 331 + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate 332 - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets 333 334 Level: beginner 335 M*/ 336 #undef __FUNCT__ 337 #define __FUNCT__ "TaoCreate_TRON" 338 PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao) 339 { 340 TAO_TRON *tron; 341 PetscErrorCode ierr; 342 const char *morethuente_type = TAOLINESEARCHMT; 343 344 PetscFunctionBegin; 345 tao->ops->setup = TaoSetup_TRON; 346 tao->ops->solve = TaoSolve_TRON; 347 tao->ops->view = TaoView_TRON; 348 tao->ops->setfromoptions = TaoSetFromOptions_TRON; 349 tao->ops->destroy = TaoDestroy_TRON; 350 tao->ops->computedual = TaoComputeDual_TRON; 351 352 ierr = PetscNewLog(tao,&tron);CHKERRQ(ierr); 353 tao->data = (void*)tron; 354 355 /* Override default settings (unless already changed) */ 356 if (!tao->max_it_changed) tao->max_it = 50; 357 358 #if defined(PETSC_USE_REAL_SINGLE) 359 if (!tao->steptol_changed) tao->steptol = 1.0e-6; 360 #else 361 if (!tao->steptol_changed) tao->steptol = 1.0e-12; 362 #endif 363 364 if (!tao->trust0_changed) tao->trust0 = 1.0; 365 366 /* Initialize pointers and variables */ 367 tron->n = 0; 368 tron->maxgpits = 3; 369 tron->pg_ftol = 0.001; 370 371 tron->eta1 = 1.0e-4; 372 tron->eta2 = 0.25; 373 tron->eta3 = 0.50; 374 tron->eta4 = 0.90; 375 376 tron->sigma1 = 0.5; 377 tron->sigma2 = 2.0; 378 tron->sigma3 = 4.0; 379 380 tron->gp_iterates = 0; /* Cumulative number */ 381 tron->total_gp_its = 0; 382 tron->n_free = 0; 383 384 tron->DXFree=NULL; 385 tron->R=NULL; 386 tron->X_New=NULL; 387 tron->G_New=NULL; 388 tron->Work=NULL; 389 tron->Free_Local=NULL; 390 tron->H_sub=NULL; 391 tron->Hpre_sub=NULL; 392 tao->subset_type = TAO_SUBSET_SUBVEC; 393 394 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr); 395 ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr); 396 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr); 397 ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr); 398 399 ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr); 400 ierr = KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);CHKERRQ(ierr); 401 ierr = KSPSetType(tao->ksp,KSPCGSTCG);CHKERRQ(ierr); 402 PetscFunctionReturn(0); 403 } 404 405