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