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