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