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