1 #include "../src/tao/matrix/lmvmmat.h" 2 #include "ntr.h" 3 4 #include "petscksp.h" 5 #include "petscpc.h" 6 #include "petsc-private/kspimpl.h" 7 #include "petsc-private/pcimpl.h" 8 9 #define NTR_KSP_NASH 0 10 #define NTR_KSP_STCG 1 11 #define NTR_KSP_GLTR 2 12 #define NTR_KSP_TYPES 3 13 14 #define NTR_PC_NONE 0 15 #define NTR_PC_AHESS 1 16 #define NTR_PC_BFGS 2 17 #define NTR_PC_PETSC 3 18 #define NTR_PC_TYPES 4 19 20 #define BFGS_SCALE_AHESS 0 21 #define BFGS_SCALE_BFGS 1 22 #define BFGS_SCALE_TYPES 2 23 24 #define NTR_INIT_CONSTANT 0 25 #define NTR_INIT_DIRECTION 1 26 #define NTR_INIT_INTERPOLATION 2 27 #define NTR_INIT_TYPES 3 28 29 #define NTR_UPDATE_REDUCTION 0 30 #define NTR_UPDATE_INTERPOLATION 1 31 #define NTR_UPDATE_TYPES 2 32 33 static const char *NTR_KSP[64] = { "nash", "stcg", "gltr"}; 34 35 static const char *NTR_PC[64] = { "none", "ahess", "bfgs", "petsc"}; 36 37 static const char *BFGS_SCALE[64] = { "ahess", "bfgs"}; 38 39 static const char *NTR_INIT[64] = { "constant", "direction", "interpolation"}; 40 41 static const char *NTR_UPDATE[64] = { "reduction", "interpolation"}; 42 43 /* Routine for BFGS preconditioner */ 44 static PetscErrorCode MatLMVMSolveShell(PC pc, Vec xin, Vec xout); 45 46 /* 47 TaoSolve_NTR - Implements Newton's Method with a trust region approach 48 for solving unconstrained minimization problems. 49 50 The basic algorithm is taken from MINPACK-2 (dstrn). 51 52 TaoSolve_NTR computes a local minimizer of a twice differentiable function 53 f by applying a trust region variant of Newton's method. At each stage 54 of the algorithm, we use the prconditioned conjugate gradient method to 55 determine an approximate minimizer of the quadratic equation 56 57 q(s) = <s, Hs + g> 58 59 subject to the trust region constraint 60 61 || s ||_M <= radius, 62 63 where radius is the trust region radius and M is a symmetric positive 64 definite matrix (the preconditioner). Here g is the gradient and H 65 is the Hessian matrix. 66 67 Note: TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG, 68 or KSPGLTR. Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this 69 routine regardless of what the user may have previously specified. 70 */ 71 #undef __FUNCT__ 72 #define __FUNCT__ "TaoSolve_NTR" 73 static PetscErrorCode TaoSolve_NTR(TaoSolver tao) 74 { 75 TAO_NTR *tr = (TAO_NTR *)tao->data; 76 PC pc; 77 KSPConvergedReason ksp_reason; 78 TaoSolverTerminationReason reason; 79 MatStructure matflag; 80 PetscReal fmin, ftrial, prered, actred, kappa, sigma, beta; 81 PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius; 82 PetscReal f, gnorm; 83 84 PetscReal delta; 85 PetscReal norm_d; 86 PetscErrorCode ierr; 87 88 PetscInt iter = 0; 89 PetscInt bfgsUpdates = 0; 90 PetscInt needH; 91 92 PetscInt i_max = 5; 93 PetscInt j_max = 1; 94 PetscInt i, j, N, n, its; 95 96 PetscFunctionBegin; 97 if (tao->XL || tao->XU || tao->ops->computebounds) { 98 ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");CHKERRQ(ierr); 99 } 100 101 tao->trust = tao->trust0; 102 103 /* Modify the radius if it is too large or small */ 104 tao->trust = PetscMax(tao->trust, tr->min_radius); 105 tao->trust = PetscMin(tao->trust, tr->max_radius); 106 107 108 if (NTR_PC_BFGS == tr->pc_type && !tr->M) { 109 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 110 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 111 ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tr->M);CHKERRQ(ierr); 112 ierr = MatLMVMAllocateVectors(tr->M,tao->solution);CHKERRQ(ierr); 113 } 114 115 /* Check convergence criteria */ 116 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr); 117 ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 118 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 119 needH = 1; 120 121 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr); 122 if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 123 124 /* Create vectors for the limited memory preconditioner */ 125 if ((NTR_PC_BFGS == tr->pc_type) && 126 (BFGS_SCALE_BFGS != tr->bfgs_scale_type)) { 127 if (!tr->Diag) { 128 ierr = VecDuplicate(tao->solution, &tr->Diag);CHKERRQ(ierr); 129 } 130 } 131 132 switch(tr->ksp_type) { 133 case NTR_KSP_NASH: 134 ierr = KSPSetType(tao->ksp, KSPNASH);CHKERRQ(ierr); 135 if (tao->ksp->ops->setfromoptions) { 136 (*tao->ksp->ops->setfromoptions)(tao->ksp); 137 } 138 break; 139 140 case NTR_KSP_STCG: 141 ierr = KSPSetType(tao->ksp, KSPSTCG);CHKERRQ(ierr); 142 if (tao->ksp->ops->setfromoptions) { 143 (*tao->ksp->ops->setfromoptions)(tao->ksp); 144 } 145 break; 146 147 default: 148 ierr = KSPSetType(tao->ksp, KSPGLTR);CHKERRQ(ierr); 149 if (tao->ksp->ops->setfromoptions) { 150 (*tao->ksp->ops->setfromoptions)(tao->ksp); 151 } 152 break; 153 } 154 155 /* Modify the preconditioner to use the bfgs approximation */ 156 ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr); 157 switch(tr->pc_type) { 158 case NTR_PC_NONE: 159 ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr); 160 if (pc->ops->setfromoptions) { 161 (*pc->ops->setfromoptions)(pc); 162 } 163 break; 164 165 case NTR_PC_AHESS: 166 ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr); 167 if (pc->ops->setfromoptions) { 168 (*pc->ops->setfromoptions)(pc); 169 } 170 ierr = PCJacobiSetUseAbs(pc);CHKERRQ(ierr); 171 break; 172 173 case NTR_PC_BFGS: 174 ierr = PCSetType(pc, PCSHELL);CHKERRQ(ierr); 175 if (pc->ops->setfromoptions) { 176 (*pc->ops->setfromoptions)(pc); 177 } 178 ierr = PCShellSetName(pc, "bfgs");CHKERRQ(ierr); 179 ierr = PCShellSetContext(pc, tr->M);CHKERRQ(ierr); 180 ierr = PCShellSetApply(pc, MatLMVMSolveShell);CHKERRQ(ierr); 181 break; 182 183 default: 184 /* Use the pc method set by pc_type */ 185 break; 186 } 187 188 /* Initialize trust-region radius */ 189 switch(tr->init_type) { 190 case NTR_INIT_CONSTANT: 191 /* Use the initial radius specified */ 192 break; 193 194 case NTR_INIT_INTERPOLATION: 195 /* Use the initial radius specified */ 196 max_radius = 0.0; 197 198 for (j = 0; j < j_max; ++j) { 199 fmin = f; 200 sigma = 0.0; 201 202 if (needH) { 203 ierr = TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag);CHKERRQ(ierr); 204 needH = 0; 205 } 206 207 for (i = 0; i < i_max; ++i) { 208 209 ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr); 210 ierr = VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);CHKERRQ(ierr); 211 ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr); 212 213 if (PetscIsInfOrNanReal(ftrial)) { 214 tau = tr->gamma1_i; 215 } 216 else { 217 if (ftrial < fmin) { 218 fmin = ftrial; 219 sigma = -tao->trust / gnorm; 220 } 221 222 ierr = MatMult(tao->hessian, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 223 ierr = VecDot(tao->gradient, tao->stepdirection, &prered);CHKERRQ(ierr); 224 225 prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm)); 226 actred = f - ftrial; 227 if ((PetscAbsScalar(actred) <= tr->epsilon) && 228 (PetscAbsScalar(prered) <= tr->epsilon)) { 229 kappa = 1.0; 230 } 231 else { 232 kappa = actred / prered; 233 } 234 235 tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred); 236 tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred); 237 tau_min = PetscMin(tau_1, tau_2); 238 tau_max = PetscMax(tau_1, tau_2); 239 240 if (PetscAbsScalar(kappa - 1.0) <= tr->mu1_i) { 241 /* Great agreement */ 242 max_radius = PetscMax(max_radius, tao->trust); 243 244 if (tau_max < 1.0) { 245 tau = tr->gamma3_i; 246 } 247 else if (tau_max > tr->gamma4_i) { 248 tau = tr->gamma4_i; 249 } 250 else { 251 tau = tau_max; 252 } 253 } 254 else if (PetscAbsScalar(kappa - 1.0) <= tr->mu2_i) { 255 /* Good agreement */ 256 max_radius = PetscMax(max_radius, tao->trust); 257 258 if (tau_max < tr->gamma2_i) { 259 tau = tr->gamma2_i; 260 } 261 else if (tau_max > tr->gamma3_i) { 262 tau = tr->gamma3_i; 263 } 264 else { 265 tau = tau_max; 266 } 267 } 268 else { 269 /* Not good agreement */ 270 if (tau_min > 1.0) { 271 tau = tr->gamma2_i; 272 } 273 else if (tau_max < tr->gamma1_i) { 274 tau = tr->gamma1_i; 275 } 276 else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) { 277 tau = tr->gamma1_i; 278 } 279 else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) && 280 ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) { 281 tau = tau_1; 282 } 283 else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) && 284 ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) { 285 tau = tau_2; 286 } 287 else { 288 tau = tau_max; 289 } 290 } 291 } 292 tao->trust = tau * tao->trust; 293 } 294 295 if (fmin < f) { 296 f = fmin; 297 ierr = VecAXPY(tao->solution, sigma, tao->gradient);CHKERRQ(ierr); 298 ierr = TaoComputeGradient(tao,tao->solution, tao->gradient);CHKERRQ(ierr); 299 300 ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr); 301 302 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 303 needH = 1; 304 305 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr); 306 if (reason != TAO_CONTINUE_ITERATING) { 307 PetscFunctionReturn(0); 308 } 309 } 310 } 311 tao->trust = PetscMax(tao->trust, max_radius); 312 313 /* Modify the radius if it is too large or small */ 314 tao->trust = PetscMax(tao->trust, tr->min_radius); 315 tao->trust = PetscMin(tao->trust, tr->max_radius); 316 break; 317 318 default: 319 /* Norm of the first direction will initialize radius */ 320 tao->trust = 0.0; 321 break; 322 } 323 324 /* Set initial scaling for the BFGS preconditioner 325 This step is done after computing the initial trust-region radius 326 since the function value may have decreased */ 327 if (NTR_PC_BFGS == tr->pc_type) { 328 if (f != 0.0) { 329 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 330 } 331 else { 332 delta = 2.0 / (gnorm*gnorm); 333 } 334 ierr = MatLMVMSetDelta(tr->M,delta);CHKERRQ(ierr); 335 } 336 337 /* Have not converged; continue with Newton method */ 338 while (reason == TAO_CONTINUE_ITERATING) { 339 ++iter; 340 341 /* Compute the Hessian */ 342 if (needH) { 343 ierr = TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag);CHKERRQ(ierr); 344 needH = 0; 345 } 346 347 if (NTR_PC_BFGS == tr->pc_type) { 348 if (BFGS_SCALE_AHESS == tr->bfgs_scale_type) { 349 /* Obtain diagonal for the bfgs preconditioner */ 350 ierr = MatGetDiagonal(tao->hessian, tr->Diag);CHKERRQ(ierr); 351 ierr = VecAbs(tr->Diag);CHKERRQ(ierr); 352 ierr = VecReciprocal(tr->Diag);CHKERRQ(ierr); 353 ierr = MatLMVMSetScale(tr->M,tr->Diag);CHKERRQ(ierr); 354 } 355 356 /* Update the limited memory preconditioner */ 357 ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr); 358 ++bfgsUpdates; 359 } 360 361 while (reason == TAO_CONTINUE_ITERATING) { 362 ierr = KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre, matflag);CHKERRQ(ierr); 363 364 /* Solve the trust region subproblem */ 365 if (NTR_KSP_NASH == tr->ksp_type) { 366 ierr = KSPNASHSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 367 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 368 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 369 tao->ksp_its+=its; 370 ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 371 } else if (NTR_KSP_STCG == tr->ksp_type) { 372 ierr = KSPSTCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 373 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 374 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 375 tao->ksp_its+=its; 376 ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 377 } else { /* NTR_KSP_GLTR */ 378 ierr = KSPGLTRSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 379 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 380 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 381 tao->ksp_its+=its; 382 ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 383 } 384 385 if (0.0 == tao->trust) { 386 /* Radius was uninitialized; use the norm of the direction */ 387 if (norm_d > 0.0) { 388 tao->trust = norm_d; 389 390 /* Modify the radius if it is too large or small */ 391 tao->trust = PetscMax(tao->trust, tr->min_radius); 392 tao->trust = PetscMin(tao->trust, tr->max_radius); 393 } 394 else { 395 /* The direction was bad; set radius to default value and re-solve 396 the trust-region subproblem to get a direction */ 397 tao->trust = tao->trust0; 398 399 /* Modify the radius if it is too large or small */ 400 tao->trust = PetscMax(tao->trust, tr->min_radius); 401 tao->trust = PetscMin(tao->trust, tr->max_radius); 402 403 if (NTR_KSP_NASH == tr->ksp_type) { 404 ierr = KSPNASHSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 405 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 406 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 407 tao->ksp_its+=its; 408 ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 409 } else if (NTR_KSP_STCG == tr->ksp_type) { 410 ierr = KSPSTCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 411 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 412 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 413 tao->ksp_its+=its; 414 ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 415 } else { /* NTR_KSP_GLTR */ 416 ierr = KSPGLTRSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 417 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 418 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 419 tao->ksp_its+=its; 420 ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 421 } 422 423 if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero"); 424 } 425 } 426 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 427 ierr = KSPGetConvergedReason(tao->ksp, &ksp_reason);CHKERRQ(ierr); 428 if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && 429 (NTR_PC_BFGS == tr->pc_type) && (bfgsUpdates > 1)) { 430 /* Preconditioner is numerically indefinite; reset the 431 approximate if using BFGS preconditioning. */ 432 433 if (f != 0.0) { 434 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 435 } 436 else { 437 delta = 2.0 / (gnorm*gnorm); 438 } 439 ierr = MatLMVMSetDelta(tr->M, delta);CHKERRQ(ierr); 440 ierr = MatLMVMReset(tr->M);CHKERRQ(ierr); 441 ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr); 442 bfgsUpdates = 1; 443 } 444 445 if (NTR_UPDATE_REDUCTION == tr->update_type) { 446 /* Get predicted reduction */ 447 if (NTR_KSP_NASH == tr->ksp_type) { 448 ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 449 } else if (NTR_KSP_STCG == tr->ksp_type) { 450 ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 451 } else { /* gltr */ 452 ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 453 } 454 455 if (prered >= 0.0) { 456 /* The predicted reduction has the wrong sign. This cannot 457 happen in infinite precision arithmetic. Step should 458 be rejected! */ 459 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 460 } 461 else { 462 /* Compute trial step and function value */ 463 ierr = VecCopy(tao->solution,tr->W);CHKERRQ(ierr); 464 ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr); 465 ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr); 466 467 if (PetscIsInfOrNanReal(ftrial)) { 468 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 469 } else { 470 /* Compute and actual reduction */ 471 actred = f - ftrial; 472 prered = -prered; 473 if ((PetscAbsScalar(actred) <= tr->epsilon) && 474 (PetscAbsScalar(prered) <= tr->epsilon)) { 475 kappa = 1.0; 476 } 477 else { 478 kappa = actred / prered; 479 } 480 481 /* Accept or reject the step and update radius */ 482 if (kappa < tr->eta1) { 483 /* Reject the step */ 484 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 485 } 486 else { 487 /* Accept the step */ 488 if (kappa < tr->eta2) { 489 /* Marginal bad step */ 490 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d); 491 } 492 else if (kappa < tr->eta3) { 493 /* Reasonable step */ 494 tao->trust = tr->alpha3 * tao->trust; 495 } 496 else if (kappa < tr->eta4) { 497 /* Good step */ 498 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust); 499 } 500 else { 501 /* Very good step */ 502 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust); 503 } 504 break; 505 } 506 } 507 } 508 } 509 else { 510 /* Get predicted reduction */ 511 if (NTR_KSP_NASH == tr->ksp_type) { 512 ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 513 } else if (NTR_KSP_STCG == tr->ksp_type) { 514 ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 515 } else { /* gltr */ 516 ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 517 } 518 519 if (prered >= 0.0) { 520 /* The predicted reduction has the wrong sign. This cannot 521 happen in infinite precision arithmetic. Step should 522 be rejected! */ 523 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 524 } 525 else { 526 ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr); 527 ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr); 528 ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr); 529 if (PetscIsInfOrNanReal(ftrial)) { 530 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 531 } 532 else { 533 ierr = VecDot(tao->gradient, tao->stepdirection, &beta);CHKERRQ(ierr); 534 actred = f - ftrial; 535 prered = -prered; 536 if ((PetscAbsScalar(actred) <= tr->epsilon) && 537 (PetscAbsScalar(prered) <= tr->epsilon)) { 538 kappa = 1.0; 539 } 540 else { 541 kappa = actred / prered; 542 } 543 544 tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred); 545 tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred); 546 tau_min = PetscMin(tau_1, tau_2); 547 tau_max = PetscMax(tau_1, tau_2); 548 549 if (kappa >= 1.0 - tr->mu1) { 550 /* Great agreement; accept step and update radius */ 551 if (tau_max < 1.0) { 552 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d); 553 } 554 else if (tau_max > tr->gamma4) { 555 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d); 556 } 557 else { 558 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 559 } 560 break; 561 } 562 else if (kappa >= 1.0 - tr->mu2) { 563 /* Good agreement */ 564 565 if (tau_max < tr->gamma2) { 566 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d); 567 } 568 else if (tau_max > tr->gamma3) { 569 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d); 570 } 571 else if (tau_max < 1.0) { 572 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 573 } 574 else { 575 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 576 } 577 break; 578 } 579 else { 580 /* Not good agreement */ 581 if (tau_min > 1.0) { 582 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d); 583 } 584 else if (tau_max < tr->gamma1) { 585 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 586 } 587 else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) { 588 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 589 } 590 else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) && 591 ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) { 592 tao->trust = tau_1 * PetscMin(tao->trust, norm_d); 593 } 594 else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) && 595 ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) { 596 tao->trust = tau_2 * PetscMin(tao->trust, norm_d); 597 } 598 else { 599 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 600 } 601 } 602 } 603 } 604 } 605 606 /* The step computed was not good and the radius was decreased. 607 Monitor the radius to terminate. */ 608 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, tao->trust, &reason);CHKERRQ(ierr); 609 } 610 611 /* The radius may have been increased; modify if it is too large */ 612 tao->trust = PetscMin(tao->trust, tr->max_radius); 613 614 if (reason == TAO_CONTINUE_ITERATING) { 615 ierr = VecCopy(tr->W, tao->solution);CHKERRQ(ierr); 616 f = ftrial; 617 ierr = TaoComputeGradient(tao, tao->solution, tao->gradient); 618 ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr); 619 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 620 needH = 1; 621 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, tao->trust, &reason);CHKERRQ(ierr); 622 } 623 } 624 PetscFunctionReturn(0); 625 } 626 627 /*------------------------------------------------------------*/ 628 #undef __FUNCT__ 629 #define __FUNCT__ "TaoSetUp_NTR" 630 static PetscErrorCode TaoSetUp_NTR(TaoSolver tao) 631 { 632 TAO_NTR *tr = (TAO_NTR *)tao->data; 633 PetscErrorCode ierr; 634 635 PetscFunctionBegin; 636 637 if (!tao->gradient) {ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr);} 638 if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);} 639 if (!tr->W) {ierr = VecDuplicate(tao->solution, &tr->W);CHKERRQ(ierr);} 640 641 tr->Diag = 0; 642 tr->M = 0; 643 644 645 PetscFunctionReturn(0); 646 } 647 648 /*------------------------------------------------------------*/ 649 #undef __FUNCT__ 650 #define __FUNCT__ "TaoDestroy_NTR" 651 static PetscErrorCode TaoDestroy_NTR(TaoSolver tao) 652 { 653 TAO_NTR *tr = (TAO_NTR *)tao->data; 654 PetscErrorCode ierr; 655 656 PetscFunctionBegin; 657 if (tao->setupcalled) { 658 ierr = VecDestroy(&tr->W);CHKERRQ(ierr); 659 } 660 ierr = MatDestroy(&tr->M);CHKERRQ(ierr); 661 ierr = VecDestroy(&tr->Diag);CHKERRQ(ierr); 662 ierr = PetscFree(tao->data);CHKERRQ(ierr); 663 PetscFunctionReturn(0); 664 } 665 666 /*------------------------------------------------------------*/ 667 #undef __FUNCT__ 668 #define __FUNCT__ "TaoSetFromOptions_NTR" 669 static PetscErrorCode TaoSetFromOptions_NTR(TaoSolver tao) 670 { 671 TAO_NTR *tr = (TAO_NTR *)tao->data; 672 PetscErrorCode ierr; 673 674 PetscFunctionBegin; 675 ierr = PetscOptionsHead("Newton trust region method for unconstrained optimization");CHKERRQ(ierr); 676 ierr = PetscOptionsEList("-tao_ntr_ksp_type", "ksp type", "", NTR_KSP, NTR_KSP_TYPES, NTR_KSP[tr->ksp_type], &tr->ksp_type, 0);CHKERRQ(ierr); 677 ierr = PetscOptionsEList("-tao_ntr_pc_type", "pc type", "", NTR_PC, NTR_PC_TYPES, NTR_PC[tr->pc_type], &tr->pc_type, 0);CHKERRQ(ierr); 678 ierr = PetscOptionsEList("-tao_ntr_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tr->bfgs_scale_type], &tr->bfgs_scale_type, 0);CHKERRQ(ierr); 679 ierr = PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type, 0);CHKERRQ(ierr); 680 ierr = PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type, 0);CHKERRQ(ierr); 681 ierr = PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1, 0);CHKERRQ(ierr); 682 ierr = PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2, 0);CHKERRQ(ierr); 683 ierr = PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3, 0);CHKERRQ(ierr); 684 ierr = PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4, 0);CHKERRQ(ierr); 685 ierr = PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1, 0);CHKERRQ(ierr); 686 ierr = PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2, 0);CHKERRQ(ierr); 687 ierr = PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3, 0);CHKERRQ(ierr); 688 ierr = PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4, 0);CHKERRQ(ierr); 689 ierr = PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5, 0);CHKERRQ(ierr); 690 ierr = PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1, 0);CHKERRQ(ierr); 691 ierr = PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2, 0);CHKERRQ(ierr); 692 ierr = PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1, 0);CHKERRQ(ierr); 693 ierr = PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2, 0);CHKERRQ(ierr); 694 ierr = PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3, 0);CHKERRQ(ierr); 695 ierr = PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4, 0);CHKERRQ(ierr); 696 ierr = PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta, 0);CHKERRQ(ierr); 697 ierr = PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i, 0);CHKERRQ(ierr); 698 ierr = PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i, 0);CHKERRQ(ierr); 699 ierr = PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i, 0);CHKERRQ(ierr); 700 ierr = PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i, 0);CHKERRQ(ierr); 701 ierr = PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i, 0);CHKERRQ(ierr); 702 ierr = PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i, 0);CHKERRQ(ierr); 703 ierr = PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i, 0);CHKERRQ(ierr); 704 ierr = PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius, 0);CHKERRQ(ierr); 705 ierr = PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius, 0);CHKERRQ(ierr); 706 ierr = PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon, 0);CHKERRQ(ierr); 707 ierr = PetscOptionsTail();CHKERRQ(ierr); 708 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 709 PetscFunctionReturn(0); 710 } 711 712 /*------------------------------------------------------------*/ 713 #undef __FUNCT__ 714 #define __FUNCT__ "TaoView_NTR" 715 static PetscErrorCode TaoView_NTR(TaoSolver tao, PetscViewer viewer) 716 { 717 TAO_NTR *tr = (TAO_NTR *)tao->data; 718 PetscErrorCode ierr; 719 PetscInt nrejects; 720 PetscBool isascii; 721 722 PetscFunctionBegin; 723 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 724 if (isascii) { 725 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 726 if (NTR_PC_BFGS == tr->pc_type && tr->M) { 727 ierr = MatLMVMGetRejects(tr->M, &nrejects);CHKERRQ(ierr); 728 ierr = PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n", nrejects);CHKERRQ(ierr); 729 } 730 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 731 } 732 PetscFunctionReturn(0); 733 } 734 735 /*------------------------------------------------------------*/ 736 EXTERN_C_BEGIN 737 #undef __FUNCT__ 738 #define __FUNCT__ "TaoCreate_NTR" 739 PetscErrorCode TaoCreate_NTR(TaoSolver tao) 740 { 741 TAO_NTR *tr; 742 PetscErrorCode ierr; 743 744 PetscFunctionBegin; 745 746 ierr = PetscNewLog(tao,&tr);CHKERRQ(ierr); 747 748 tao->ops->setup = TaoSetUp_NTR; 749 tao->ops->solve = TaoSolve_NTR; 750 tao->ops->view = TaoView_NTR; 751 tao->ops->setfromoptions = TaoSetFromOptions_NTR; 752 tao->ops->destroy = TaoDestroy_NTR; 753 754 tao->max_it = 50; 755 tao->fatol = 1e-10; 756 tao->frtol = 1e-10; 757 tao->data = (void*)tr; 758 759 tao->trust0 = 100.0; 760 761 /* Standard trust region update parameters */ 762 tr->eta1 = 1.0e-4; 763 tr->eta2 = 0.25; 764 tr->eta3 = 0.50; 765 tr->eta4 = 0.90; 766 767 tr->alpha1 = 0.25; 768 tr->alpha2 = 0.50; 769 tr->alpha3 = 1.00; 770 tr->alpha4 = 2.00; 771 tr->alpha5 = 4.00; 772 773 /* Interpolation parameters */ 774 tr->mu1_i = 0.35; 775 tr->mu2_i = 0.50; 776 777 tr->gamma1_i = 0.0625; 778 tr->gamma2_i = 0.50; 779 tr->gamma3_i = 2.00; 780 tr->gamma4_i = 5.00; 781 782 tr->theta_i = 0.25; 783 784 /* Interpolation trust region update parameters */ 785 tr->mu1 = 0.10; 786 tr->mu2 = 0.50; 787 788 tr->gamma1 = 0.25; 789 tr->gamma2 = 0.50; 790 tr->gamma3 = 2.00; 791 tr->gamma4 = 4.00; 792 793 tr->theta = 0.05; 794 795 tr->min_radius = 1.0e-10; 796 tr->max_radius = 1.0e10; 797 tr->epsilon = 1.0e-6; 798 799 tr->ksp_type = NTR_KSP_STCG; 800 tr->pc_type = NTR_PC_BFGS; 801 tr->bfgs_scale_type = BFGS_SCALE_AHESS; 802 tr->init_type = NTR_INIT_INTERPOLATION; 803 tr->update_type = NTR_UPDATE_REDUCTION; 804 805 806 /* Set linear solver to default for trust region */ 807 ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr); 808 809 PetscFunctionReturn(0); 810 811 812 } 813 EXTERN_C_END 814 815 816 #undef __FUNCT__ 817 #define __FUNCT__ "MatLMVMSolveShell" 818 static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x) 819 { 820 PetscErrorCode ierr; 821 Mat M; 822 PetscFunctionBegin; 823 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 824 PetscValidHeaderSpecific(b,VEC_CLASSID,2); 825 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 826 ierr = PCShellGetContext(pc,(void**)&M);CHKERRQ(ierr); 827 ierr = MatLMVMSolve(M, b, x);CHKERRQ(ierr); 828 PetscFunctionReturn(0); 829 } 830