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 19 ierr = VecDestroy(&tron->X_New);CHKERRQ(ierr); 20 ierr = VecDestroy(&tron->G_New);CHKERRQ(ierr); 21 ierr = VecDestroy(&tron->Work);CHKERRQ(ierr); 22 ierr = VecDestroy(&tron->DXFree);CHKERRQ(ierr); 23 ierr = VecDestroy(&tron->R); CHKERRQ(ierr); 24 ierr = VecDestroy(&tron->diag); CHKERRQ(ierr); 25 ierr = VecScatterDestroy(&tron->scatter); CHKERRQ(ierr); 26 ierr = ISDestroy(&tron->Free_Local); CHKERRQ(ierr); 27 ierr = MatDestroy(&tron->H_sub); CHKERRQ(ierr); 28 ierr = MatDestroy(&tron->Hpre_sub); CHKERRQ(ierr); 29 ierr = PetscFree(tao->data); 30 tao->data = PETSC_NULL; 31 32 PetscFunctionReturn(0); 33 } 34 35 /*------------------------------------------------------------*/ 36 #undef __FUNCT__ 37 #define __FUNCT__ "TaoSetFromOptions_TRON" 38 static PetscErrorCode TaoSetFromOptions_TRON(TaoSolver tao) 39 { 40 TAO_TRON *tron = (TAO_TRON *)tao->data; 41 PetscErrorCode ierr; 42 PetscBool flg; 43 44 PetscFunctionBegin; 45 46 ierr = PetscOptionsHead("Newton Trust Region Method for bound constrained optimization");CHKERRQ(ierr); 47 48 ierr = PetscOptionsInt("-tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg); 49 CHKERRQ(ierr); 50 51 ierr = PetscOptionsTail();CHKERRQ(ierr); 52 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 53 ierr = KSPSetFromOptions(tao->ksp); CHKERRQ(ierr); 54 55 PetscFunctionReturn(0); 56 } 57 58 /*------------------------------------------------------------*/ 59 #undef __FUNCT__ 60 #define __FUNCT__ "TaoView_TRON" 61 static PetscErrorCode TaoView_TRON(TaoSolver tao, PetscViewer viewer) 62 { 63 TAO_TRON *tron = (TAO_TRON *)tao->data; 64 PetscBool isascii; 65 PetscErrorCode ierr; 66 67 PetscFunctionBegin; 68 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 69 if (isascii) { 70 ierr = PetscViewerASCIIPushTab(viewer); CHKERRQ(ierr); 71 ierr = PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);CHKERRQ(ierr); 72 ierr = PetscViewerASCIIPrintf(viewer,"PG tolerance: %G \n",tron->pg_ftol);CHKERRQ(ierr); 73 ierr = PetscViewerASCIIPopTab(viewer); CHKERRQ(ierr); 74 } else { 75 SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO TRON",((PetscObject)viewer)->type_name); 76 } 77 PetscFunctionReturn(0); 78 } 79 80 81 /* ---------------------------------------------------------- */ 82 #undef __FUNCT__ 83 #define __FUNCT__ "TaoSetup_TRON" 84 static PetscErrorCode TaoSetup_TRON(TaoSolver tao) 85 { 86 PetscErrorCode ierr; 87 TAO_TRON *tron = (TAO_TRON *)tao->data; 88 89 PetscFunctionBegin; 90 91 /* Allocate some arrays */ 92 ierr = VecDuplicate(tao->solution, &tron->diag); CHKERRQ(ierr); 93 ierr = VecDuplicate(tao->solution, &tron->X_New); CHKERRQ(ierr); 94 ierr = VecDuplicate(tao->solution, &tron->G_New); CHKERRQ(ierr); 95 ierr = VecDuplicate(tao->solution, &tron->Work); CHKERRQ(ierr); 96 ierr = VecDuplicate(tao->solution, &tao->gradient); CHKERRQ(ierr); 97 ierr = VecDuplicate(tao->solution, &tao->stepdirection); CHKERRQ(ierr); 98 if (!tao->XL) { 99 ierr = VecDuplicate(tao->solution, &tao->XL); CHKERRQ(ierr); 100 ierr = VecSet(tao->XL, TAO_NINFINITY); CHKERRQ(ierr); 101 } 102 if (!tao->XU) { 103 ierr = VecDuplicate(tao->solution, &tao->XU); CHKERRQ(ierr); 104 ierr = VecSet(tao->XU, TAO_INFINITY); CHKERRQ(ierr); 105 } 106 107 PetscFunctionReturn(0); 108 } 109 110 111 112 #undef __FUNCT__ 113 #define __FUNCT__ "TaoSolve_TRON" 114 static PetscErrorCode TaoSolve_TRON(TaoSolver tao){ 115 116 TAO_TRON *tron = (TAO_TRON *)tao->data; 117 PetscErrorCode ierr; 118 PetscInt iter=0,its; 119 120 TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING; 121 TaoLineSearchTerminationReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING; 122 PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize; 123 PetscFunctionBegin; 124 125 tron->pgstepsize=1.0; 126 127 tao->trust = tao->trust0; 128 /* Project the current point onto the feasible set */ 129 ierr = TaoComputeVariableBounds(tao); CHKERRQ(ierr); 130 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution); CHKERRQ(ierr); 131 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU); CHKERRQ(ierr); 132 133 134 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);CHKERRQ(ierr); 135 if (tron->Free_Local) { 136 ierr = ISDestroy(&tron->Free_Local); CHKERRQ(ierr); 137 tron->Free_Local = PETSC_NULL; 138 } 139 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local); CHKERRQ(ierr); 140 141 /* Project the gradient and calculate the norm */ 142 ierr = VecBoundGradientProjection(tao->gradient,tao->solution, tao->XL, tao->XU, tao->gradient); CHKERRQ(ierr); 143 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm); CHKERRQ(ierr); 144 145 if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) { 146 SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN"); 147 } 148 149 if (tao->trust <= 0) { 150 tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0); 151 } 152 153 tron->stepsize=tao->trust; 154 ierr = TaoMonitor(tao, iter, tron->f, tron->gnorm, 0.0, tron->stepsize, &reason); CHKERRQ(ierr); 155 while (reason==TAO_CONTINUE_ITERATING){ 156 157 ierr = TronGradientProjections(tao,tron); CHKERRQ(ierr); 158 f=tron->f; delta=tao->trust; 159 tron->n_free_last = tron->n_free; 160 ierr = TaoComputeHessian(tao,tao->solution,&tao->hessian, &tao->hessian_pre, &tron->matflag);CHKERRQ(ierr); 161 162 ierr = ISGetSize(tron->Free_Local, &tron->n_free); CHKERRQ(ierr); 163 164 /* If no free variables */ 165 if (tron->n_free == 0) { 166 actred=0; 167 PetscInfo(tao,"No free variables in tron iteration."); 168 break; 169 170 } 171 /* use free_local to mask/submat gradient, hessian, stepdirection */ 172 ierr = VecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R); 173 CHKERRQ(ierr); 174 ierr = VecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree); 175 CHKERRQ(ierr); 176 ierr = VecSet(tron->DXFree,0.0); CHKERRQ(ierr); 177 ierr = VecScale(tron->R, -1.0); CHKERRQ(ierr); 178 ierr = MatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub); CHKERRQ(ierr); 179 if (tao->hessian == tao->hessian_pre) { 180 ierr = MatDestroy(&tron->Hpre_sub); CHKERRQ(ierr); 181 ierr = PetscObjectReference((PetscObject)(tron->H_sub)); CHKERRQ(ierr); 182 tron->Hpre_sub = tron->H_sub; 183 } else { 184 ierr = MatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub); 185 CHKERRQ(ierr); 186 } 187 ierr = KSPReset(tao->ksp); CHKERRQ(ierr); 188 ierr = KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub, tron->matflag); CHKERRQ(ierr); 189 while (1) { 190 191 /* Approximately solve the reduced linear system */ 192 ierr = KSPSTCGSetRadius(tao->ksp,delta); CHKERRQ(ierr); 193 ierr = KSPSolve(tao->ksp, tron->R, tron->DXFree); CHKERRQ(ierr); 194 ierr = KSPGetIterationNumber(tao->ksp,&its); CHKERRQ(ierr); 195 tao->ksp_its+=its; 196 ierr = VecSet(tao->stepdirection,0.0); CHKERRQ(ierr); 197 198 /* Add dxfree matrix to compute step direction vector */ 199 ierr = VecReducedXPY(tao->stepdirection,tron->DXFree,tron->Free_Local); CHKERRQ(ierr); 200 if (0) { 201 PetscReal rhs,stepnorm; 202 ierr = VecNorm(tron->R,NORM_2,&rhs); CHKERRQ(ierr); 203 ierr = VecNorm(tron->DXFree,NORM_2,&stepnorm); CHKERRQ(ierr); 204 ierr = PetscPrintf(PETSC_COMM_WORLD,"|rhs|=%g\t|s|=%g\n",rhs,stepnorm); CHKERRQ(ierr); 205 } 206 207 208 ierr = VecDot(tao->gradient, tao->stepdirection, &gdx); CHKERRQ(ierr); 209 ierr = PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",gdx); 210 CHKERRQ(ierr); 211 212 ierr = VecCopy(tao->solution, tron->X_New); CHKERRQ(ierr); 213 ierr = VecCopy(tao->gradient, tron->G_New); CHKERRQ(ierr); 214 215 stepsize=1.0;f_new=f; 216 217 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0); CHKERRQ(ierr); 218 ierr = TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection, 219 &stepsize,&ls_reason); CHKERRQ(ierr); CHKERRQ(ierr); 220 ierr = TaoAddLineSearchCounts(tao); CHKERRQ(ierr); 221 222 ierr = MatMult(tao->hessian, tao->stepdirection, tron->Work); CHKERRQ(ierr); 223 ierr = VecAYPX(tron->Work, 0.5, tao->gradient); CHKERRQ(ierr); 224 ierr = VecDot(tao->stepdirection, tron->Work, &prered); CHKERRQ(ierr); 225 actred = f_new - f; 226 if (actred<0) { 227 rhok=PetscAbs(-actred/prered); 228 } else { 229 rhok=0.0; 230 } 231 232 /* Compare actual improvement to the quadratic model */ 233 if (rhok > tron->eta1) { /* Accept the point */ 234 /* d = x_new - x */ 235 ierr = VecCopy(tron->X_New, tao->stepdirection); CHKERRQ(ierr); 236 ierr = VecAXPY(tao->stepdirection, -1.0, tao->solution); CHKERRQ(ierr); 237 238 ierr = VecNorm(tao->stepdirection, NORM_2, &xdiff); CHKERRQ(ierr); 239 xdiff *= stepsize; 240 241 /* Adjust trust region size */ 242 if (rhok < tron->eta2 ){ 243 delta = PetscMin(xdiff,delta)*tron->sigma1; 244 } else if (rhok > tron->eta4 ){ 245 delta= PetscMin(xdiff,delta)*tron->sigma3; 246 } else if (rhok > tron->eta3 ){ 247 delta=PetscMin(xdiff,delta)*tron->sigma2; 248 } 249 ierr = VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient); CHKERRQ(ierr); 250 if (tron->Free_Local) { 251 ierr = ISDestroy(&tron->Free_Local); CHKERRQ(ierr); 252 tron->Free_Local=PETSC_NULL; 253 } 254 ierr = VecWhichBetween(tao->XL, tron->X_New, tao->XU, &tron->Free_Local); CHKERRQ(ierr); 255 f=f_new; 256 ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm); CHKERRQ(ierr); 257 ierr = VecCopy(tron->X_New, tao->solution); CHKERRQ(ierr); 258 ierr = VecCopy(tron->G_New, tao->gradient); CHKERRQ(ierr); 259 break; 260 } 261 else if (delta <= 1e-30) { 262 break; 263 } 264 else { 265 delta /= 4.0; 266 } 267 } /* end linear solve loop */ 268 269 270 tron->f=f; tron->actred=actred; tao->trust=delta; 271 iter++; 272 ierr = TaoMonitor(tao, iter, tron->f, tron->gnorm, 0.0, delta, &reason); CHKERRQ(ierr); 273 } /* END MAIN LOOP */ 274 275 PetscFunctionReturn(0); 276 } 277 278 279 #undef __FUNCT__ 280 #define __FUNCT__ "TronGradientProjections" 281 static PetscErrorCode TronGradientProjections(TaoSolver tao,TAO_TRON *tron) 282 { 283 PetscErrorCode ierr; 284 PetscInt i; 285 TaoLineSearchTerminationReason ls_reason; 286 PetscReal actred=-1.0,actred_max=0.0; 287 PetscReal f_new; 288 /* 289 The gradient and function value passed into and out of this 290 routine should be current and correct. 291 292 The free, active, and binding variables should be already identified 293 */ 294 295 PetscFunctionBegin; 296 if (tron->Free_Local) { 297 ierr = ISDestroy(&tron->Free_Local); CHKERRQ(ierr); 298 tron->Free_Local = PETSC_NULL; 299 } 300 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local); CHKERRQ(ierr); 301 302 for (i=0;i<tron->maxgpits;i++){ 303 304 if ( -actred <= (tron->pg_ftol)*actred_max) break; 305 306 tron->gp_iterates++; tron->total_gp_its++; 307 f_new=tron->f; 308 309 ierr = VecCopy(tao->gradient, tao->stepdirection); CHKERRQ(ierr); 310 ierr = VecScale(tao->stepdirection, -1.0); CHKERRQ(ierr); 311 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize); CHKERRQ(ierr); 312 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, 313 &tron->pgstepsize, &ls_reason); CHKERRQ(ierr); 314 ierr = TaoAddLineSearchCounts(tao); CHKERRQ(ierr); 315 316 317 /* Update the iterate */ 318 actred = f_new - tron->f; 319 actred_max = PetscMax(actred_max,-(f_new - tron->f)); 320 tron->f = f_new; 321 if (tron->Free_Local) { 322 ierr = ISDestroy(&tron->Free_Local); CHKERRQ(ierr); 323 tron->Free_Local = PETSC_NULL; 324 } 325 ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local); CHKERRQ(ierr); 326 } 327 328 PetscFunctionReturn(0); 329 } 330 331 #undef __FUNCT__ 332 #define __FUNCT__ "TaoComputeDual_TRON" 333 static PetscErrorCode TaoComputeDual_TRON(TaoSolver tao, Vec DXL, Vec DXU) { 334 335 TAO_TRON *tron = (TAO_TRON *)tao->data; 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 340 PetscValidHeaderSpecific(tao,TAOSOLVER_CLASSID,1); 341 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 342 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 343 344 if (!tron->Work || !tao->gradient) { 345 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n"); 346 } 347 348 ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work); CHKERRQ(ierr); 349 ierr = VecCopy(tron->Work,DXL); CHKERRQ(ierr); 350 ierr = VecAXPY(DXL,-1.0,tao->gradient); CHKERRQ(ierr); 351 ierr = VecSet(DXU,0.0); CHKERRQ(ierr); 352 ierr = VecPointwiseMax(DXL,DXL,DXU); CHKERRQ(ierr); 353 354 ierr = VecCopy(tao->gradient,DXU); CHKERRQ(ierr); 355 ierr = VecAXPY(DXU,-1.0,tron->Work); CHKERRQ(ierr); 356 ierr = VecSet(tron->Work,0.0); CHKERRQ(ierr); 357 ierr = VecPointwiseMin(DXU,tron->Work,DXU); CHKERRQ(ierr); 358 359 PetscFunctionReturn(0); 360 } 361 362 /*------------------------------------------------------------*/ 363 EXTERN_C_BEGIN 364 #undef __FUNCT__ 365 #define __FUNCT__ "TaoCreate_TRON" 366 PetscErrorCode TaoCreate_TRON(TaoSolver tao) 367 { 368 TAO_TRON *tron; 369 PetscErrorCode ierr; 370 const char *morethuente_type = TAOLINESEARCH_MT; 371 PetscFunctionBegin; 372 373 tao->ops->setup = TaoSetup_TRON; 374 tao->ops->solve = TaoSolve_TRON; 375 tao->ops->view = TaoView_TRON; 376 tao->ops->setfromoptions = TaoSetFromOptions_TRON; 377 tao->ops->destroy = TaoDestroy_TRON; 378 tao->ops->computedual = TaoComputeDual_TRON; 379 380 ierr = PetscNewLog(tao,TAO_TRON,&tron); CHKERRQ(ierr); 381 382 tao->max_it = 50; 383 tao->fatol = 1e-10; 384 tao->frtol = 1e-10; 385 tao->data = (void*)tron; 386 tao->steptol = 1e-12; 387 tao->trust0 = 1.0; 388 389 /* Initialize pointers and variables */ 390 tron->n = 0; 391 tron->maxgpits = 3; 392 tron->pg_ftol = 0.001; 393 394 tron->eta1 = 1.0e-4; 395 tron->eta2 = 0.25; 396 tron->eta3 = 0.50; 397 tron->eta4 = 0.90; 398 399 tron->sigma1 = 0.5; 400 tron->sigma2 = 2.0; 401 tron->sigma3 = 4.0; 402 403 tron->gp_iterates = 0; /* Cumulative number */ 404 tron->total_gp_its = 0; 405 406 tron->n_free = 0; 407 408 tron->DXFree=PETSC_NULL; 409 tron->R=PETSC_NULL; 410 tron->X_New=PETSC_NULL; 411 tron->G_New=PETSC_NULL; 412 tron->Work=PETSC_NULL; 413 tron->Free_Local=PETSC_NULL; 414 tron->H_sub=PETSC_NULL; 415 tron->Hpre_sub=PETSC_NULL; 416 tao->subset_type = TAO_SUBSET_SUBVEC; 417 418 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); CHKERRQ(ierr); 419 ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type); CHKERRQ(ierr); 420 ierr = TaoLineSearchUseTaoSolverRoutines(tao->linesearch,tao); CHKERRQ(ierr); 421 422 ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp); CHKERRQ(ierr); 423 ierr = KSPSetType(tao->ksp,KSPSTCG); CHKERRQ(ierr); 424 425 PetscFunctionReturn(0); 426 } 427 EXTERN_C_END 428