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