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