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