1 #include "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(TaoSolver,TAO_TRON*); 9 /*------------------------------------------------------------*/ 10 #undef __FUNCT__ 11 #define __FUNCT__ "TaoDestroy_TRON" 12 static PetscErrorCode TaoDestroy_TRON(TaoSolver 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 tao->data = PETSC_NULL; 30 PetscFunctionReturn(0); 31 } 32 33 /*------------------------------------------------------------*/ 34 #undef __FUNCT__ 35 #define __FUNCT__ "TaoSetFromOptions_TRON" 36 static PetscErrorCode TaoSetFromOptions_TRON(TaoSolver tao) 37 { 38 TAO_TRON *tron = (TAO_TRON *)tao->data; 39 PetscErrorCode ierr; 40 PetscBool flg; 41 42 PetscFunctionBegin; 43 ierr = PetscOptionsHead("Newton Trust Region Method for bound constrained optimization");CHKERRQ(ierr); 44 ierr = PetscOptionsInt("-tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);CHKERRQ(ierr); 45 ierr = PetscOptionsTail();CHKERRQ(ierr); 46 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 47 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 48 PetscFunctionReturn(0); 49 } 50 51 /*------------------------------------------------------------*/ 52 #undef __FUNCT__ 53 #define __FUNCT__ "TaoView_TRON" 54 static PetscErrorCode TaoView_TRON(TaoSolver tao, PetscViewer viewer) 55 { 56 TAO_TRON *tron = (TAO_TRON *)tao->data; 57 PetscBool isascii; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 62 if (isascii) { 63 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 64 ierr = PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);CHKERRQ(ierr); 65 ierr = PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);CHKERRQ(ierr); 66 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 67 } 68 PetscFunctionReturn(0); 69 } 70 71 72 /* ---------------------------------------------------------- */ 73 #undef __FUNCT__ 74 #define __FUNCT__ "TaoSetup_TRON" 75 static PetscErrorCode TaoSetup_TRON(TaoSolver tao) 76 { 77 PetscErrorCode ierr; 78 TAO_TRON *tron = (TAO_TRON *)tao->data; 79 80 PetscFunctionBegin; 81 82 /* Allocate some arrays */ 83 ierr = VecDuplicate(tao->solution, &tron->diag);CHKERRQ(ierr); 84 ierr = VecDuplicate(tao->solution, &tron->X_New);CHKERRQ(ierr); 85 ierr = VecDuplicate(tao->solution, &tron->G_New);CHKERRQ(ierr); 86 ierr = VecDuplicate(tao->solution, &tron->Work);CHKERRQ(ierr); 87 ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr); 88 ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr); 89 if (!tao->XL) { 90 ierr = VecDuplicate(tao->solution, &tao->XL);CHKERRQ(ierr); 91 ierr = VecSet(tao->XL, TAO_NINFINITY);CHKERRQ(ierr); 92 } 93 if (!tao->XU) { 94 ierr = VecDuplicate(tao->solution, &tao->XU);CHKERRQ(ierr); 95 ierr = VecSet(tao->XU, TAO_INFINITY);CHKERRQ(ierr); 96 } 97 98 PetscFunctionReturn(0); 99 } 100 101 102 103 #undef __FUNCT__ 104 #define __FUNCT__ "TaoSolve_TRON" 105 static PetscErrorCode TaoSolve_TRON(TaoSolver tao) 106 { 107 TAO_TRON *tron = (TAO_TRON *)tao->data; 108 PetscErrorCode ierr; 109 PetscInt iter=0,its; 110 TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING; 111 TaoLineSearchTerminationReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING; 112 PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize; 113 114 PetscFunctionBegin; 115 tron->pgstepsize=1.0; 116 tao->trust = tao->trust0; 117 /* Project the current point onto the feasible set */ 118 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 119 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 120 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 121 122 123 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);CHKERRQ(ierr); 124 if (tron->Free_Local) { 125 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 126 tron->Free_Local = PETSC_NULL; 127 } 128 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 129 130 /* Project the gradient and calculate the norm */ 131 ierr = VecBoundGradientProjection(tao->gradient,tao->solution, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 132 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr); 133 134 if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) { 135 SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN"); 136 } 137 138 if (tao->trust <= 0) { 139 tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0); 140 } 141 142 tron->stepsize=tao->trust; 143 ierr = TaoMonitor(tao, iter, tron->f, tron->gnorm, 0.0, tron->stepsize, &reason);CHKERRQ(ierr); 144 while (reason==TAO_CONTINUE_ITERATING){ 145 146 ierr = TronGradientProjections(tao,tron);CHKERRQ(ierr); 147 f=tron->f; delta=tao->trust; 148 tron->n_free_last = tron->n_free; 149 ierr = TaoComputeHessian(tao,tao->solution,&tao->hessian, &tao->hessian_pre, &tron->matflag);CHKERRQ(ierr); 150 151 ierr = ISGetSize(tron->Free_Local, &tron->n_free);CHKERRQ(ierr); 152 153 /* If no free variables */ 154 if (tron->n_free == 0) { 155 actred=0; 156 PetscInfo(tao,"No free variables in tron iteration."); 157 break; 158 159 } 160 /* use free_local to mask/submat gradient, hessian, stepdirection */ 161 ierr = VecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);CHKERRQ(ierr); 162 ierr = VecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);CHKERRQ(ierr); 163 ierr = VecSet(tron->DXFree,0.0);CHKERRQ(ierr); 164 ierr = VecScale(tron->R, -1.0);CHKERRQ(ierr); 165 ierr = MatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);CHKERRQ(ierr); 166 if (tao->hessian == tao->hessian_pre) { 167 ierr = MatDestroy(&tron->Hpre_sub);CHKERRQ(ierr); 168 ierr = PetscObjectReference((PetscObject)(tron->H_sub));CHKERRQ(ierr); 169 tron->Hpre_sub = tron->H_sub; 170 } else { 171 ierr = MatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);CHKERRQ(ierr); 172 } 173 ierr = KSPReset(tao->ksp);CHKERRQ(ierr); 174 ierr = KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub, tron->matflag);CHKERRQ(ierr); 175 while (1) { 176 177 /* Approximately solve the reduced linear system */ 178 ierr = KSPSTCGSetRadius(tao->ksp,delta);CHKERRQ(ierr); 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 ierr = VecSet(tao->stepdirection,0.0);CHKERRQ(ierr); 183 184 /* Add dxfree matrix to compute step direction vector */ 185 ierr = VecReducedXPY(tao->stepdirection,tron->DXFree,tron->Free_Local);CHKERRQ(ierr); 186 if (0) { 187 PetscReal rhs,stepnorm; 188 ierr = VecNorm(tron->R,NORM_2,&rhs);CHKERRQ(ierr); 189 ierr = VecNorm(tron->DXFree,NORM_2,&stepnorm);CHKERRQ(ierr); 190 ierr = PetscPrintf(PETSC_COMM_WORLD,"|rhs|=%g\t|s|=%g\n",rhs,stepnorm);CHKERRQ(ierr); 191 } 192 193 194 ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr); 195 ierr = PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",gdx);CHKERRQ(ierr); 196 197 ierr = VecCopy(tao->solution, tron->X_New);CHKERRQ(ierr); 198 ierr = VecCopy(tao->gradient, tron->G_New);CHKERRQ(ierr); 199 200 stepsize=1.0;f_new=f; 201 202 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr); 203 ierr = TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);CHKERRQ(ierr);CHKERRQ(ierr); 204 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 205 206 ierr = MatMult(tao->hessian, tao->stepdirection, tron->Work);CHKERRQ(ierr); 207 ierr = VecAYPX(tron->Work, 0.5, tao->gradient);CHKERRQ(ierr); 208 ierr = VecDot(tao->stepdirection, tron->Work, &prered);CHKERRQ(ierr); 209 actred = f_new - f; 210 if (actred<0) { 211 rhok=PetscAbs(-actred/prered); 212 } else { 213 rhok=0.0; 214 } 215 216 /* Compare actual improvement to the quadratic model */ 217 if (rhok > tron->eta1) { /* Accept the point */ 218 /* d = x_new - x */ 219 ierr = VecCopy(tron->X_New, tao->stepdirection);CHKERRQ(ierr); 220 ierr = VecAXPY(tao->stepdirection, -1.0, tao->solution);CHKERRQ(ierr); 221 222 ierr = VecNorm(tao->stepdirection, NORM_2, &xdiff);CHKERRQ(ierr); 223 xdiff *= stepsize; 224 225 /* Adjust trust region size */ 226 if (rhok < tron->eta2 ){ 227 delta = PetscMin(xdiff,delta)*tron->sigma1; 228 } else if (rhok > tron->eta4 ){ 229 delta= PetscMin(xdiff,delta)*tron->sigma3; 230 } else if (rhok > tron->eta3 ){ 231 delta=PetscMin(xdiff,delta)*tron->sigma2; 232 } 233 ierr = VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 234 if (tron->Free_Local) { 235 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 236 tron->Free_Local=PETSC_NULL; 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 iter++; 256 ierr = TaoMonitor(tao, iter, 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(TaoSolver tao,TAO_TRON *tron) 266 { 267 PetscErrorCode ierr; 268 PetscInt i; 269 TaoLineSearchTerminationReason 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 tron->Free_Local = PETSC_NULL; 282 } 283 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr); 284 285 for (i=0;i<tron->maxgpits;i++){ 286 287 if ( -actred <= (tron->pg_ftol)*actred_max) break; 288 289 tron->gp_iterates++; tron->total_gp_its++; 290 f_new=tron->f; 291 292 ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr); 293 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 294 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);CHKERRQ(ierr); 295 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, 296 &tron->pgstepsize, &ls_reason);CHKERRQ(ierr); 297 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 298 299 300 /* Update the iterate */ 301 actred = f_new - tron->f; 302 actred_max = PetscMax(actred_max,-(f_new - tron->f)); 303 tron->f = f_new; 304 if (tron->Free_Local) { 305 ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr); 306 tron->Free_Local = PETSC_NULL; 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(TaoSolver tao, Vec DXL, Vec DXU) { 317 318 TAO_TRON *tron = (TAO_TRON *)tao->data; 319 PetscErrorCode ierr; 320 321 PetscFunctionBegin; 322 PetscValidHeaderSpecific(tao,TAOSOLVER_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(TaoSolver tao) 345 { 346 TAO_TRON *tron; 347 PetscErrorCode ierr; 348 const char *morethuente_type = TAOLINESEARCH_MT; 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 tao->fatol = 1e-10; 362 tao->frtol = 1e-10; 363 tao->data = (void*)tron; 364 tao->steptol = 1e-12; 365 tao->trust0 = 1.0; 366 367 /* Initialize pointers and variables */ 368 tron->n = 0; 369 tron->maxgpits = 3; 370 tron->pg_ftol = 0.001; 371 372 tron->eta1 = 1.0e-4; 373 tron->eta2 = 0.25; 374 tron->eta3 = 0.50; 375 tron->eta4 = 0.90; 376 377 tron->sigma1 = 0.5; 378 tron->sigma2 = 2.0; 379 tron->sigma3 = 4.0; 380 381 tron->gp_iterates = 0; /* Cumulative number */ 382 tron->total_gp_its = 0; 383 tron->n_free = 0; 384 385 tron->DXFree=PETSC_NULL; 386 tron->R=PETSC_NULL; 387 tron->X_New=PETSC_NULL; 388 tron->G_New=PETSC_NULL; 389 tron->Work=PETSC_NULL; 390 tron->Free_Local=PETSC_NULL; 391 tron->H_sub=PETSC_NULL; 392 tron->Hpre_sub=PETSC_NULL; 393 tao->subset_type = TAO_SUBSET_SUBVEC; 394 395 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr); 396 ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr); 397 ierr = TaoLineSearchUseTaoSolverRoutines(tao->linesearch,tao);CHKERRQ(ierr); 398 399 ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr); 400 ierr = KSPSetType(tao->ksp,KSPSTCG);CHKERRQ(ierr); 401 PetscFunctionReturn(0); 402 } 403 EXTERN_C_END 404