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 PetscBool is_stcg,is_nash,is_gltr; 109 TaoConvergedReason reason = TAO_CONTINUE_ITERATING; 110 TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING; 111 PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize; 112 113 PetscFunctionBegin; 114 tron->pgstepsize=1.0; 115 tao->trust = tao->trust0; 116 /* Project the current point onto the feasible set */ 117 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 118 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 119 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 120 121 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);CHKERRQ(ierr); 122 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 123 124 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 125 126 /* Project the gradient and calculate the norm */ 127 ierr = VecBoundGradientProjection(tao->gradient,tao->solution, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 128 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 129 130 if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN"); 131 if (tao->trust <= 0) { 132 tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0); 133 } 134 135 tron->stepsize=tao->trust; 136 ierr = TaoMonitor(tao, iter, tron->f, tron->gnorm, 0.0, tron->stepsize, &reason);CHKERRQ(ierr); 137 while (reason==TAO_CONTINUE_ITERATING){ 138 139 ierr = TronGradientProjections(tao,tron);CHKERRQ(ierr); 140 f=tron->f; delta=tao->trust; 141 tron->n_free_last = tron->n_free; 142 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 143 144 ierr = ISGetSize(tron->Free_Local, &tron->n_free);CHKERRQ(ierr); 145 146 /* If no free variables */ 147 if (tron->n_free == 0) { 148 actred=0; 149 PetscInfo(tao,"No free variables in tron iteration."); 150 break; 151 152 } 153 /* use free_local to mask/submat gradient, hessian, stepdirection */ 154 ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);CHKERRQ(ierr); 155 ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);CHKERRQ(ierr); 156 ierr = VecSet(tron->DXFree,0.0);CHKERRQ(ierr); 157 ierr = VecScale(tron->R, -1.0);CHKERRQ(ierr); 158 ierr = TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);CHKERRQ(ierr); 159 if (tao->hessian == tao->hessian_pre) { 160 ierr = MatDestroy(&tron->Hpre_sub);CHKERRQ(ierr); 161 ierr = PetscObjectReference((PetscObject)(tron->H_sub));CHKERRQ(ierr); 162 tron->Hpre_sub = tron->H_sub; 163 } else { 164 ierr = TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);CHKERRQ(ierr); 165 } 166 ierr = KSPReset(tao->ksp);CHKERRQ(ierr); 167 ierr = KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);CHKERRQ(ierr); 168 ierr = PetscObjectTypeCompare((PetscObject)(tao->ksp),KSPSTCG,&is_stcg);CHKERRQ(ierr); 169 ierr = PetscObjectTypeCompare((PetscObject)(tao->ksp),KSPNASH,&is_nash);CHKERRQ(ierr); 170 ierr = PetscObjectTypeCompare((PetscObject)(tao->ksp),KSPGLTR,&is_gltr);CHKERRQ(ierr); 171 while (1) { 172 173 /* Approximately solve the reduced linear system */ 174 if (is_stcg) { 175 ierr = KSPSTCGSetRadius(tao->ksp,delta);CHKERRQ(ierr); 176 } else if (is_nash) { 177 ierr = KSPNASHSetRadius(tao->ksp,delta);CHKERRQ(ierr); 178 } else if (is_gltr) { 179 ierr = KSPGLTRSetRadius(tao->ksp,delta);CHKERRQ(ierr); 180 } 181 182 ierr = KSPSolve(tao->ksp, tron->R, tron->DXFree);CHKERRQ(ierr); 183 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 184 tao->ksp_its+=its; 185 ierr = VecSet(tao->stepdirection,0.0);CHKERRQ(ierr); 186 187 /* Add dxfree matrix to compute step direction vector */ 188 ierr = VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);CHKERRQ(ierr); 189 if (0) { 190 PetscReal rhs,stepnorm; 191 ierr = VecNorm(tron->R,NORM_2,&rhs);CHKERRQ(ierr); 192 ierr = VecNorm(tron->DXFree,NORM_2,&stepnorm);CHKERRQ(ierr); 193 ierr = PetscPrintf(PETSC_COMM_WORLD,"|rhs|=%g\t|s|=%g\n",(double)rhs,(double)stepnorm);CHKERRQ(ierr); 194 } 195 196 197 ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr); 198 ierr = PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",(double)gdx);CHKERRQ(ierr); 199 200 ierr = VecCopy(tao->solution, tron->X_New);CHKERRQ(ierr); 201 ierr = VecCopy(tao->gradient, tron->G_New);CHKERRQ(ierr); 202 203 stepsize=1.0;f_new=f; 204 205 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr); 206 ierr = TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);CHKERRQ(ierr);CHKERRQ(ierr); 207 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 208 209 ierr = MatMult(tao->hessian, tao->stepdirection, tron->Work);CHKERRQ(ierr); 210 ierr = VecAYPX(tron->Work, 0.5, tao->gradient);CHKERRQ(ierr); 211 ierr = VecDot(tao->stepdirection, tron->Work, &prered);CHKERRQ(ierr); 212 actred = f_new - f; 213 if (actred<0) { 214 rhok=PetscAbs(-actred/prered); 215 } else { 216 rhok=0.0; 217 } 218 219 /* Compare actual improvement to the quadratic model */ 220 if (rhok > tron->eta1) { /* Accept the point */ 221 /* d = x_new - x */ 222 ierr = VecCopy(tron->X_New, tao->stepdirection);CHKERRQ(ierr); 223 ierr = VecAXPY(tao->stepdirection, -1.0, tao->solution);CHKERRQ(ierr); 224 225 ierr = VecNorm(tao->stepdirection, NORM_2, &xdiff);CHKERRQ(ierr); 226 xdiff *= stepsize; 227 228 /* Adjust trust region size */ 229 if (rhok < tron->eta2 ){ 230 delta = PetscMin(xdiff,delta)*tron->sigma1; 231 } else if (rhok > tron->eta4 ){ 232 delta= PetscMin(xdiff,delta)*tron->sigma3; 233 } else if (rhok > tron->eta3 ){ 234 delta=PetscMin(xdiff,delta)*tron->sigma2; 235 } 236 ierr = VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 237 if (tron->Free_Local) { 238 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 239 } 240 ierr = VecWhichBetween(tao->XL, tron->X_New, tao->XU, &tron->Free_Local);CHKERRQ(ierr); 241 f=f_new; 242 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 243 ierr = VecCopy(tron->X_New, tao->solution);CHKERRQ(ierr); 244 ierr = VecCopy(tron->G_New, tao->gradient);CHKERRQ(ierr); 245 break; 246 } 247 else if (delta <= 1e-30) { 248 break; 249 } 250 else { 251 delta /= 4.0; 252 } 253 } /* end linear solve loop */ 254 255 256 tron->f=f; tron->actred=actred; tao->trust=delta; 257 iter++; 258 ierr = TaoMonitor(tao, iter, tron->f, tron->gnorm, 0.0, delta, &reason);CHKERRQ(ierr); 259 } /* END MAIN LOOP */ 260 261 PetscFunctionReturn(0); 262 } 263 264 265 #undef __FUNCT__ 266 #define __FUNCT__ "TronGradientProjections" 267 static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron) 268 { 269 PetscErrorCode ierr; 270 PetscInt i; 271 TaoLineSearchConvergedReason ls_reason; 272 PetscReal actred=-1.0,actred_max=0.0; 273 PetscReal f_new; 274 /* 275 The gradient and function value passed into and out of this 276 routine should be current and correct. 277 278 The free, active, and binding variables should be already identified 279 */ 280 PetscFunctionBegin; 281 if (tron->Free_Local) { 282 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 283 } 284 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 285 286 for (i=0;i<tron->maxgpits;i++){ 287 288 if ( -actred <= (tron->pg_ftol)*actred_max) break; 289 290 tron->gp_iterates++; tron->total_gp_its++; 291 f_new=tron->f; 292 293 ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr); 294 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 295 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);CHKERRQ(ierr); 296 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, 297 &tron->pgstepsize, &ls_reason);CHKERRQ(ierr); 298 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 299 300 301 /* Update the iterate */ 302 actred = f_new - tron->f; 303 actred_max = PetscMax(actred_max,-(f_new - tron->f)); 304 tron->f = f_new; 305 if (tron->Free_Local) { 306 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 307 } 308 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 309 } 310 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "TaoComputeDual_TRON" 316 static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) { 317 318 TAO_TRON *tron = (TAO_TRON *)tao->data; 319 PetscErrorCode ierr; 320 321 PetscFunctionBegin; 322 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 323 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 324 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 325 if (!tron->Work || !tao->gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n"); 326 327 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);CHKERRQ(ierr); 328 ierr = VecCopy(tron->Work,DXL);CHKERRQ(ierr); 329 ierr = VecAXPY(DXL,-1.0,tao->gradient);CHKERRQ(ierr); 330 ierr = VecSet(DXU,0.0);CHKERRQ(ierr); 331 ierr = VecPointwiseMax(DXL,DXL,DXU);CHKERRQ(ierr); 332 333 ierr = VecCopy(tao->gradient,DXU);CHKERRQ(ierr); 334 ierr = VecAXPY(DXU,-1.0,tron->Work);CHKERRQ(ierr); 335 ierr = VecSet(tron->Work,0.0);CHKERRQ(ierr); 336 ierr = VecPointwiseMin(DXU,tron->Work,DXU);CHKERRQ(ierr); 337 PetscFunctionReturn(0); 338 } 339 340 /*------------------------------------------------------------*/ 341 EXTERN_C_BEGIN 342 #undef __FUNCT__ 343 #define __FUNCT__ "TaoCreate_TRON" 344 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 EXTERN_C_END 410