1 #include <../src/tao/unconstrained/impls/ntr/ntrimpl.h> 2 3 #include <petscksp.h> 4 5 #define NTR_INIT_CONSTANT 0 6 #define NTR_INIT_DIRECTION 1 7 #define NTR_INIT_INTERPOLATION 2 8 #define NTR_INIT_TYPES 3 9 10 #define NTR_UPDATE_REDUCTION 0 11 #define NTR_UPDATE_INTERPOLATION 1 12 #define NTR_UPDATE_TYPES 2 13 14 static const char *NTR_INIT[64] = {"constant","direction","interpolation"}; 15 16 static const char *NTR_UPDATE[64] = {"reduction","interpolation"}; 17 18 /* 19 TaoSolve_NTR - Implements Newton's Method with a trust region approach 20 for solving unconstrained minimization problems. 21 22 The basic algorithm is taken from MINPACK-2 (dstrn). 23 24 TaoSolve_NTR computes a local minimizer of a twice differentiable function 25 f by applying a trust region variant of Newton's method. At each stage 26 of the algorithm, we use the prconditioned conjugate gradient method to 27 determine an approximate minimizer of the quadratic equation 28 29 q(s) = <s, Hs + g> 30 31 subject to the trust region constraint 32 33 || s ||_M <= radius, 34 35 where radius is the trust region radius and M is a symmetric positive 36 definite matrix (the preconditioner). Here g is the gradient and H 37 is the Hessian matrix. 38 39 Note: TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG, 40 or KSPGLTR. Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this 41 routine regardless of what the user may have previously specified. 42 */ 43 static PetscErrorCode TaoSolve_NTR(Tao tao) 44 { 45 TAO_NTR *tr = (TAO_NTR *)tao->data; 46 KSPType ksp_type; 47 PetscBool is_nash,is_stcg,is_gltr,is_bfgs,is_jacobi,is_symmetric,sym_set; 48 KSPConvergedReason ksp_reason; 49 PC pc; 50 PetscReal fmin, ftrial, prered, actred, kappa, sigma, beta; 51 PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius; 52 PetscReal f, gnorm; 53 54 PetscReal norm_d; 55 PetscInt bfgsUpdates = 0; 56 PetscInt needH; 57 58 PetscInt i_max = 5; 59 PetscInt j_max = 1; 60 PetscInt i, j, N, n, its; 61 62 PetscFunctionBegin; 63 if (tao->XL || tao->XU || tao->ops->computebounds) { 64 PetscCall(PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n")); 65 } 66 67 PetscCall(KSPGetType(tao->ksp,&ksp_type)); 68 PetscCall(PetscStrcmp(ksp_type,KSPNASH,&is_nash)); 69 PetscCall(PetscStrcmp(ksp_type,KSPSTCG,&is_stcg)); 70 PetscCall(PetscStrcmp(ksp_type,KSPGLTR,&is_gltr)); 71 PetscCheck(is_nash || is_stcg || is_gltr,PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"TAO_NTR requires nash, stcg, or gltr for the KSP"); 72 73 /* Initialize the radius and modify if it is too large or small */ 74 tao->trust = tao->trust0; 75 tao->trust = PetscMax(tao->trust, tr->min_radius); 76 tao->trust = PetscMin(tao->trust, tr->max_radius); 77 78 /* Allocate the vectors needed for the BFGS approximation */ 79 PetscCall(KSPGetPC(tao->ksp, &pc)); 80 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs)); 81 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi)); 82 if (is_bfgs) { 83 tr->bfgs_pre = pc; 84 PetscCall(PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M)); 85 PetscCall(VecGetLocalSize(tao->solution, &n)); 86 PetscCall(VecGetSize(tao->solution, &N)); 87 PetscCall(MatSetSizes(tr->M, n, n, N, N)); 88 PetscCall(MatLMVMAllocate(tr->M, tao->solution, tao->gradient)); 89 PetscCall(MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric)); 90 PetscCheck(sym_set && is_symmetric,PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric."); 91 } else if (is_jacobi) { 92 PetscCall(PCJacobiSetUseAbs(pc,PETSC_TRUE)); 93 } 94 95 /* Check convergence criteria */ 96 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient)); 97 PetscCall(TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm)); 98 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER,"User provided compute function generated Inf or NaN"); 99 needH = 1; 100 101 tao->reason = TAO_CONTINUE_ITERATING; 102 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 103 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0)); 104 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 105 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 106 107 /* Initialize trust-region radius */ 108 switch (tr->init_type) { 109 case NTR_INIT_CONSTANT: 110 /* Use the initial radius specified */ 111 break; 112 113 case NTR_INIT_INTERPOLATION: 114 /* Use the initial radius specified */ 115 max_radius = 0.0; 116 117 for (j = 0; j < j_max; ++j) { 118 fmin = f; 119 sigma = 0.0; 120 121 if (needH) { 122 PetscCall(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 123 needH = 0; 124 } 125 126 for (i = 0; i < i_max; ++i) { 127 128 PetscCall(VecCopy(tao->solution, tr->W)); 129 PetscCall(VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient)); 130 PetscCall(TaoComputeObjective(tao, tr->W, &ftrial)); 131 132 if (PetscIsInfOrNanReal(ftrial)) { 133 tau = tr->gamma1_i; 134 } 135 else { 136 if (ftrial < fmin) { 137 fmin = ftrial; 138 sigma = -tao->trust / gnorm; 139 } 140 141 PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection)); 142 PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered)); 143 144 prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm)); 145 actred = f - ftrial; 146 if ((PetscAbsScalar(actred) <= tr->epsilon) && 147 (PetscAbsScalar(prered) <= tr->epsilon)) { 148 kappa = 1.0; 149 } 150 else { 151 kappa = actred / prered; 152 } 153 154 tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred); 155 tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred); 156 tau_min = PetscMin(tau_1, tau_2); 157 tau_max = PetscMax(tau_1, tau_2); 158 159 if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu1_i) { 160 /* Great agreement */ 161 max_radius = PetscMax(max_radius, tao->trust); 162 163 if (tau_max < 1.0) { 164 tau = tr->gamma3_i; 165 } 166 else if (tau_max > tr->gamma4_i) { 167 tau = tr->gamma4_i; 168 } 169 else { 170 tau = tau_max; 171 } 172 } 173 else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu2_i) { 174 /* Good agreement */ 175 max_radius = PetscMax(max_radius, tao->trust); 176 177 if (tau_max < tr->gamma2_i) { 178 tau = tr->gamma2_i; 179 } 180 else if (tau_max > tr->gamma3_i) { 181 tau = tr->gamma3_i; 182 } 183 else { 184 tau = tau_max; 185 } 186 } 187 else { 188 /* Not good agreement */ 189 if (tau_min > 1.0) { 190 tau = tr->gamma2_i; 191 } 192 else if (tau_max < tr->gamma1_i) { 193 tau = tr->gamma1_i; 194 } 195 else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) { 196 tau = tr->gamma1_i; 197 } 198 else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) && 199 ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) { 200 tau = tau_1; 201 } 202 else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) && 203 ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) { 204 tau = tau_2; 205 } 206 else { 207 tau = tau_max; 208 } 209 } 210 } 211 tao->trust = tau * tao->trust; 212 } 213 214 if (fmin < f) { 215 f = fmin; 216 PetscCall(VecAXPY(tao->solution, sigma, tao->gradient)); 217 PetscCall(TaoComputeGradient(tao,tao->solution, tao->gradient)); 218 219 PetscCall(TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm)); 220 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 221 needH = 1; 222 223 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 224 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0)); 225 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 226 if (tao->reason != TAO_CONTINUE_ITERATING) { 227 PetscFunctionReturn(0); 228 } 229 } 230 } 231 tao->trust = PetscMax(tao->trust, max_radius); 232 233 /* Modify the radius if it is too large or small */ 234 tao->trust = PetscMax(tao->trust, tr->min_radius); 235 tao->trust = PetscMin(tao->trust, tr->max_radius); 236 break; 237 238 default: 239 /* Norm of the first direction will initialize radius */ 240 tao->trust = 0.0; 241 break; 242 } 243 244 /* Have not converged; continue with Newton method */ 245 while (tao->reason == TAO_CONTINUE_ITERATING) { 246 /* Call general purpose update function */ 247 if (tao->ops->update) { 248 PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 249 } 250 ++tao->niter; 251 tao->ksp_its=0; 252 /* Compute the Hessian */ 253 if (needH) { 254 PetscCall(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 255 needH = 0; 256 } 257 258 if (tr->bfgs_pre) { 259 /* Update the limited memory preconditioner */ 260 PetscCall(MatLMVMUpdate(tr->M, tao->solution, tao->gradient)); 261 ++bfgsUpdates; 262 } 263 264 while (tao->reason == TAO_CONTINUE_ITERATING) { 265 PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre)); 266 267 /* Solve the trust region subproblem */ 268 PetscCall(KSPCGSetRadius(tao->ksp,tao->trust)); 269 PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection)); 270 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 271 tao->ksp_its+=its; 272 tao->ksp_tot_its+=its; 273 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 274 275 if (0.0 == tao->trust) { 276 /* Radius was uninitialized; use the norm of the direction */ 277 if (norm_d > 0.0) { 278 tao->trust = norm_d; 279 280 /* Modify the radius if it is too large or small */ 281 tao->trust = PetscMax(tao->trust, tr->min_radius); 282 tao->trust = PetscMin(tao->trust, tr->max_radius); 283 } 284 else { 285 /* The direction was bad; set radius to default value and re-solve 286 the trust-region subproblem to get a direction */ 287 tao->trust = tao->trust0; 288 289 /* Modify the radius if it is too large or small */ 290 tao->trust = PetscMax(tao->trust, tr->min_radius); 291 tao->trust = PetscMin(tao->trust, tr->max_radius); 292 293 PetscCall(KSPCGSetRadius(tao->ksp,tao->trust)); 294 PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection)); 295 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 296 tao->ksp_its+=its; 297 tao->ksp_tot_its+=its; 298 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 299 300 PetscCheck(norm_d != 0.0,PetscObjectComm((PetscObject)tao),PETSC_ERR_PLIB, "Initial direction zero"); 301 } 302 } 303 PetscCall(VecScale(tao->stepdirection, -1.0)); 304 PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason)); 305 if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) { 306 /* Preconditioner is numerically indefinite; reset the 307 approximate if using BFGS preconditioning. */ 308 PetscCall(MatLMVMReset(tr->M, PETSC_FALSE)); 309 PetscCall(MatLMVMUpdate(tr->M, tao->solution, tao->gradient)); 310 bfgsUpdates = 1; 311 } 312 313 if (NTR_UPDATE_REDUCTION == tr->update_type) { 314 /* Get predicted reduction */ 315 PetscCall(KSPCGGetObjFcn(tao->ksp,&prered)); 316 if (prered >= 0.0) { 317 /* The predicted reduction has the wrong sign. This cannot 318 happen in infinite precision arithmetic. Step should 319 be rejected! */ 320 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 321 } 322 else { 323 /* Compute trial step and function value */ 324 PetscCall(VecCopy(tao->solution,tr->W)); 325 PetscCall(VecAXPY(tr->W, 1.0, tao->stepdirection)); 326 PetscCall(TaoComputeObjective(tao, tr->W, &ftrial)); 327 328 if (PetscIsInfOrNanReal(ftrial)) { 329 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 330 } else { 331 /* Compute and actual reduction */ 332 actred = f - ftrial; 333 prered = -prered; 334 if ((PetscAbsScalar(actred) <= tr->epsilon) && 335 (PetscAbsScalar(prered) <= tr->epsilon)) { 336 kappa = 1.0; 337 } 338 else { 339 kappa = actred / prered; 340 } 341 342 /* Accept or reject the step and update radius */ 343 if (kappa < tr->eta1) { 344 /* Reject the step */ 345 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 346 } 347 else { 348 /* Accept the step */ 349 if (kappa < tr->eta2) { 350 /* Marginal bad step */ 351 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d); 352 } 353 else if (kappa < tr->eta3) { 354 /* Reasonable step */ 355 tao->trust = tr->alpha3 * tao->trust; 356 } 357 else if (kappa < tr->eta4) { 358 /* Good step */ 359 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust); 360 } 361 else { 362 /* Very good step */ 363 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust); 364 } 365 break; 366 } 367 } 368 } 369 } 370 else { 371 /* Get predicted reduction */ 372 PetscCall(KSPCGGetObjFcn(tao->ksp,&prered)); 373 if (prered >= 0.0) { 374 /* The predicted reduction has the wrong sign. This cannot 375 happen in infinite precision arithmetic. Step should 376 be rejected! */ 377 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 378 } 379 else { 380 PetscCall(VecCopy(tao->solution, tr->W)); 381 PetscCall(VecAXPY(tr->W, 1.0, tao->stepdirection)); 382 PetscCall(TaoComputeObjective(tao, tr->W, &ftrial)); 383 if (PetscIsInfOrNanReal(ftrial)) { 384 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 385 } 386 else { 387 PetscCall(VecDot(tao->gradient, tao->stepdirection, &beta)); 388 actred = f - ftrial; 389 prered = -prered; 390 if ((PetscAbsScalar(actred) <= tr->epsilon) && 391 (PetscAbsScalar(prered) <= tr->epsilon)) { 392 kappa = 1.0; 393 } 394 else { 395 kappa = actred / prered; 396 } 397 398 tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred); 399 tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred); 400 tau_min = PetscMin(tau_1, tau_2); 401 tau_max = PetscMax(tau_1, tau_2); 402 403 if (kappa >= 1.0 - tr->mu1) { 404 /* Great agreement; accept step and update radius */ 405 if (tau_max < 1.0) { 406 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d); 407 } 408 else if (tau_max > tr->gamma4) { 409 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d); 410 } 411 else { 412 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 413 } 414 break; 415 } 416 else if (kappa >= 1.0 - tr->mu2) { 417 /* Good agreement */ 418 419 if (tau_max < tr->gamma2) { 420 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d); 421 } 422 else if (tau_max > tr->gamma3) { 423 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d); 424 } 425 else if (tau_max < 1.0) { 426 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 427 } 428 else { 429 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 430 } 431 break; 432 } 433 else { 434 /* Not good agreement */ 435 if (tau_min > 1.0) { 436 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d); 437 } 438 else if (tau_max < tr->gamma1) { 439 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 440 } 441 else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) { 442 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 443 } 444 else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) && 445 ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) { 446 tao->trust = tau_1 * PetscMin(tao->trust, norm_d); 447 } 448 else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) && 449 ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) { 450 tao->trust = tau_2 * PetscMin(tao->trust, norm_d); 451 } 452 else { 453 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 454 } 455 } 456 } 457 } 458 } 459 460 /* The step computed was not good and the radius was decreased. 461 Monitor the radius to terminate. */ 462 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 463 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust)); 464 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 465 } 466 467 /* The radius may have been increased; modify if it is too large */ 468 tao->trust = PetscMin(tao->trust, tr->max_radius); 469 470 if (tao->reason == TAO_CONTINUE_ITERATING) { 471 PetscCall(VecCopy(tr->W, tao->solution)); 472 f = ftrial; 473 PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient)); 474 PetscCall(TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm)); 475 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 476 needH = 1; 477 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 478 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust)); 479 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 480 } 481 } 482 PetscFunctionReturn(0); 483 } 484 485 /*------------------------------------------------------------*/ 486 static PetscErrorCode TaoSetUp_NTR(Tao tao) 487 { 488 TAO_NTR *tr = (TAO_NTR *)tao->data; 489 490 PetscFunctionBegin; 491 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 492 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 493 if (!tr->W) PetscCall(VecDuplicate(tao->solution, &tr->W)); 494 495 tr->bfgs_pre = NULL; 496 tr->M = NULL; 497 PetscFunctionReturn(0); 498 } 499 500 /*------------------------------------------------------------*/ 501 static PetscErrorCode TaoDestroy_NTR(Tao tao) 502 { 503 TAO_NTR *tr = (TAO_NTR *)tao->data; 504 505 PetscFunctionBegin; 506 if (tao->setupcalled) { 507 PetscCall(VecDestroy(&tr->W)); 508 } 509 PetscCall(KSPDestroy(&tao->ksp)); 510 PetscCall(PetscFree(tao->data)); 511 PetscFunctionReturn(0); 512 } 513 514 /*------------------------------------------------------------*/ 515 static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao) 516 { 517 TAO_NTR *tr = (TAO_NTR *)tao->data; 518 519 PetscFunctionBegin; 520 PetscOptionsHeadBegin(PetscOptionsObject,"Newton trust region method for unconstrained optimization"); 521 PetscCall(PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type,NULL)); 522 PetscCall(PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type,NULL)); 523 PetscCall(PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL)); 524 PetscCall(PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL)); 525 PetscCall(PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL)); 526 PetscCall(PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL)); 527 PetscCall(PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL)); 528 PetscCall(PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL)); 529 PetscCall(PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL)); 530 PetscCall(PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL)); 531 PetscCall(PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL)); 532 PetscCall(PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL)); 533 PetscCall(PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL)); 534 PetscCall(PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL)); 535 PetscCall(PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL)); 536 PetscCall(PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL)); 537 PetscCall(PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL)); 538 PetscCall(PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL)); 539 PetscCall(PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL)); 540 PetscCall(PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL)); 541 PetscCall(PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL)); 542 PetscCall(PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL)); 543 PetscCall(PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL)); 544 PetscCall(PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL)); 545 PetscCall(PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL)); 546 PetscCall(PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL)); 547 PetscCall(PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL)); 548 PetscCall(PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL)); 549 PetscOptionsHeadEnd(); 550 PetscCall(KSPSetFromOptions(tao->ksp)); 551 PetscFunctionReturn(0); 552 } 553 554 /*------------------------------------------------------------*/ 555 /*MC 556 TAONTR - Newton's method with trust region for unconstrained minimization. 557 At each iteration, the Newton trust region method solves the system. 558 NTR expects a KSP solver with a trust region radius. 559 min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k 560 561 Options Database Keys: 562 + -tao_ntr_init_type - "constant","direction","interpolation" 563 . -tao_ntr_update_type - "reduction","interpolation" 564 . -tao_ntr_min_radius - lower bound on trust region radius 565 . -tao_ntr_max_radius - upper bound on trust region radius 566 . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction 567 . -tao_ntr_mu1_i - mu1 interpolation init factor 568 . -tao_ntr_mu2_i - mu2 interpolation init factor 569 . -tao_ntr_gamma1_i - gamma1 interpolation init factor 570 . -tao_ntr_gamma2_i - gamma2 interpolation init factor 571 . -tao_ntr_gamma3_i - gamma3 interpolation init factor 572 . -tao_ntr_gamma4_i - gamma4 interpolation init factor 573 . -tao_ntr_theta_i - theta1 interpolation init factor 574 . -tao_ntr_eta1 - eta1 reduction update factor 575 . -tao_ntr_eta2 - eta2 reduction update factor 576 . -tao_ntr_eta3 - eta3 reduction update factor 577 . -tao_ntr_eta4 - eta4 reduction update factor 578 . -tao_ntr_alpha1 - alpha1 reduction update factor 579 . -tao_ntr_alpha2 - alpha2 reduction update factor 580 . -tao_ntr_alpha3 - alpha3 reduction update factor 581 . -tao_ntr_alpha4 - alpha4 reduction update factor 582 . -tao_ntr_alpha4 - alpha4 reduction update factor 583 . -tao_ntr_mu1 - mu1 interpolation update 584 . -tao_ntr_mu2 - mu2 interpolation update 585 . -tao_ntr_gamma1 - gamma1 interpolcation update 586 . -tao_ntr_gamma2 - gamma2 interpolcation update 587 . -tao_ntr_gamma3 - gamma3 interpolcation update 588 . -tao_ntr_gamma4 - gamma4 interpolation update 589 - -tao_ntr_theta - theta interpolation update 590 591 Level: beginner 592 M*/ 593 594 PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao) 595 { 596 TAO_NTR *tr; 597 598 PetscFunctionBegin; 599 600 PetscCall(PetscNewLog(tao,&tr)); 601 602 tao->ops->setup = TaoSetUp_NTR; 603 tao->ops->solve = TaoSolve_NTR; 604 tao->ops->setfromoptions = TaoSetFromOptions_NTR; 605 tao->ops->destroy = TaoDestroy_NTR; 606 607 /* Override default settings (unless already changed) */ 608 if (!tao->max_it_changed) tao->max_it = 50; 609 if (!tao->trust0_changed) tao->trust0 = 100.0; 610 tao->data = (void*)tr; 611 612 /* Standard trust region update parameters */ 613 tr->eta1 = 1.0e-4; 614 tr->eta2 = 0.25; 615 tr->eta3 = 0.50; 616 tr->eta4 = 0.90; 617 618 tr->alpha1 = 0.25; 619 tr->alpha2 = 0.50; 620 tr->alpha3 = 1.00; 621 tr->alpha4 = 2.00; 622 tr->alpha5 = 4.00; 623 624 /* Interpolation trust region update parameters */ 625 tr->mu1 = 0.10; 626 tr->mu2 = 0.50; 627 628 tr->gamma1 = 0.25; 629 tr->gamma2 = 0.50; 630 tr->gamma3 = 2.00; 631 tr->gamma4 = 4.00; 632 633 tr->theta = 0.05; 634 635 /* Interpolation parameters for initialization */ 636 tr->mu1_i = 0.35; 637 tr->mu2_i = 0.50; 638 639 tr->gamma1_i = 0.0625; 640 tr->gamma2_i = 0.50; 641 tr->gamma3_i = 2.00; 642 tr->gamma4_i = 5.00; 643 644 tr->theta_i = 0.25; 645 646 tr->min_radius = 1.0e-10; 647 tr->max_radius = 1.0e10; 648 tr->epsilon = 1.0e-6; 649 650 tr->init_type = NTR_INIT_INTERPOLATION; 651 tr->update_type = NTR_UPDATE_REDUCTION; 652 653 /* Set linear solver to default for trust region */ 654 PetscCall(KSPCreate(((PetscObject)tao)->comm,&tao->ksp)); 655 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1)); 656 PetscCall(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 657 PetscCall(KSPAppendOptionsPrefix(tao->ksp,"tao_ntr_")); 658 PetscCall(KSPSetType(tao->ksp,KSPSTCG)); 659 PetscFunctionReturn(0); 660 } 661