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