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