1 #include <petsctaolinesearch.h> 2 #include <../src/tao/unconstrained/impls/nls/nlsimpl.h> 3 4 #include <petscksp.h> 5 6 #define NLS_INIT_CONSTANT 0 7 #define NLS_INIT_DIRECTION 1 8 #define NLS_INIT_INTERPOLATION 2 9 #define NLS_INIT_TYPES 3 10 11 #define NLS_UPDATE_STEP 0 12 #define NLS_UPDATE_REDUCTION 1 13 #define NLS_UPDATE_INTERPOLATION 2 14 #define NLS_UPDATE_TYPES 3 15 16 static const char *NLS_INIT[64] = {"constant", "direction", "interpolation"}; 17 18 static const char *NLS_UPDATE[64] = {"step", "reduction", "interpolation"}; 19 20 /* 21 Implements Newton's Method with a line search approach for solving 22 unconstrained minimization problems. A More'-Thuente line search 23 is used to guarantee that the bfgs preconditioner remains positive 24 definite. 25 26 The method can shift the Hessian matrix. The shifting procedure is 27 adapted from the PATH algorithm for solving complementarity 28 problems. 29 30 The linear system solve should be done with a conjugate gradient 31 method, although any method can be used. 32 */ 33 34 #define NLS_NEWTON 0 35 #define NLS_BFGS 1 36 #define NLS_GRADIENT 2 37 38 static PetscErrorCode TaoSolve_NLS(Tao tao) 39 { 40 TAO_NLS *nlsP = (TAO_NLS *)tao->data; 41 KSPType ksp_type; 42 PetscBool is_nash, is_stcg, is_gltr, is_bfgs, is_jacobi, is_symmetric, sym_set; 43 KSPConvergedReason ksp_reason; 44 PC pc; 45 TaoLineSearchConvergedReason ls_reason; 46 47 PetscReal fmin, ftrial, f_full, prered, actred, kappa, sigma; 48 PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius; 49 PetscReal f, fold, gdx, gnorm, pert; 50 PetscReal step = 1.0; 51 PetscReal norm_d = 0.0, e_min; 52 53 PetscInt stepType; 54 PetscInt bfgsUpdates = 0; 55 PetscInt n, N, kspits; 56 PetscInt needH = 1; 57 58 PetscInt i_max = 5; 59 PetscInt j_max = 1; 60 PetscInt i, j; 61 62 PetscFunctionBegin; 63 if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by nls algorithm\n")); 64 65 /* Initialized variables */ 66 pert = nlsP->sval; 67 68 /* Number of times ksp stopped because of these reasons */ 69 nlsP->ksp_atol = 0; 70 nlsP->ksp_rtol = 0; 71 nlsP->ksp_dtol = 0; 72 nlsP->ksp_ctol = 0; 73 nlsP->ksp_negc = 0; 74 nlsP->ksp_iter = 0; 75 nlsP->ksp_othr = 0; 76 77 /* Initialize trust-region radius when using nash, stcg, or gltr 78 Command automatically ignored for other methods 79 Will be reset during the first iteration 80 */ 81 PetscCall(KSPGetType(tao->ksp, &ksp_type)); 82 PetscCall(PetscStrcmp(ksp_type, KSPNASH, &is_nash)); 83 PetscCall(PetscStrcmp(ksp_type, KSPSTCG, &is_stcg)); 84 PetscCall(PetscStrcmp(ksp_type, KSPGLTR, &is_gltr)); 85 86 PetscCall(KSPCGSetRadius(tao->ksp, nlsP->max_radius)); 87 88 if (is_nash || is_stcg || is_gltr) { 89 PetscCheck(tao->trust0 >= 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_OUTOFRANGE, "Initial radius negative"); 90 tao->trust = tao->trust0; 91 tao->trust = PetscMax(tao->trust, nlsP->min_radius); 92 tao->trust = PetscMin(tao->trust, nlsP->max_radius); 93 } 94 95 /* Check convergence criteria */ 96 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient)); 97 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm)); 98 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 99 100 tao->reason = TAO_CONTINUE_ITERATING; 101 PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its)); 102 PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step)); 103 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 104 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS); 105 106 /* Allocate the vectors needed for the BFGS approximation */ 107 PetscCall(KSPGetPC(tao->ksp, &pc)); 108 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs)); 109 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi)); 110 if (is_bfgs) { 111 nlsP->bfgs_pre = pc; 112 PetscCall(PCLMVMGetMatLMVM(nlsP->bfgs_pre, &nlsP->M)); 113 PetscCall(VecGetLocalSize(tao->solution, &n)); 114 PetscCall(VecGetSize(tao->solution, &N)); 115 PetscCall(MatSetSizes(nlsP->M, n, n, N, N)); 116 PetscCall(MatLMVMAllocate(nlsP->M, tao->solution, tao->gradient)); 117 PetscCall(MatIsSymmetricKnown(nlsP->M, &sym_set, &is_symmetric)); 118 PetscCheck(sym_set && is_symmetric, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric."); 119 } else if (is_jacobi) PetscCall(PCJacobiSetUseAbs(pc, PETSC_TRUE)); 120 121 /* Initialize trust-region radius. The initialization is only performed 122 when we are using Nash, Steihaug-Toint or the Generalized Lanczos method. */ 123 if (is_nash || is_stcg || is_gltr) { 124 switch (nlsP->init_type) { 125 case NLS_INIT_CONSTANT: 126 /* Use the initial radius specified */ 127 break; 128 129 case NLS_INIT_INTERPOLATION: 130 /* Use the initial radius specified */ 131 max_radius = 0.0; 132 133 for (j = 0; j < j_max; ++j) { 134 fmin = f; 135 sigma = 0.0; 136 137 if (needH) { 138 PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre)); 139 needH = 0; 140 } 141 142 for (i = 0; i < i_max; ++i) { 143 PetscCall(VecCopy(tao->solution, nlsP->W)); 144 PetscCall(VecAXPY(nlsP->W, -tao->trust / gnorm, tao->gradient)); 145 PetscCall(TaoComputeObjective(tao, nlsP->W, &ftrial)); 146 if (PetscIsInfOrNanReal(ftrial)) { 147 tau = nlsP->gamma1_i; 148 } else { 149 if (ftrial < fmin) { 150 fmin = ftrial; 151 sigma = -tao->trust / gnorm; 152 } 153 154 PetscCall(MatMult(tao->hessian, tao->gradient, nlsP->D)); 155 PetscCall(VecDot(tao->gradient, nlsP->D, &prered)); 156 157 prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm)); 158 actred = f - ftrial; 159 if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) { 160 kappa = 1.0; 161 } else { 162 kappa = actred / prered; 163 } 164 165 tau_1 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust + (1.0 - nlsP->theta_i) * prered - actred); 166 tau_2 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust - (1.0 + nlsP->theta_i) * prered + actred); 167 tau_min = PetscMin(tau_1, tau_2); 168 tau_max = PetscMax(tau_1, tau_2); 169 170 if (PetscAbsScalar(kappa - (PetscReal)1.0) <= nlsP->mu1_i) { 171 /* Great agreement */ 172 max_radius = PetscMax(max_radius, tao->trust); 173 174 if (tau_max < 1.0) { 175 tau = nlsP->gamma3_i; 176 } else if (tau_max > nlsP->gamma4_i) { 177 tau = nlsP->gamma4_i; 178 } else if (tau_1 >= 1.0 && tau_1 <= nlsP->gamma4_i && tau_2 < 1.0) { 179 tau = tau_1; 180 } else if (tau_2 >= 1.0 && tau_2 <= nlsP->gamma4_i && tau_1 < 1.0) { 181 tau = tau_2; 182 } else { 183 tau = tau_max; 184 } 185 } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= nlsP->mu2_i) { 186 /* Good agreement */ 187 max_radius = PetscMax(max_radius, tao->trust); 188 189 if (tau_max < nlsP->gamma2_i) { 190 tau = nlsP->gamma2_i; 191 } else if (tau_max > nlsP->gamma3_i) { 192 tau = nlsP->gamma3_i; 193 } else { 194 tau = tau_max; 195 } 196 } else { 197 /* Not good agreement */ 198 if (tau_min > 1.0) { 199 tau = nlsP->gamma2_i; 200 } else if (tau_max < nlsP->gamma1_i) { 201 tau = nlsP->gamma1_i; 202 } else if ((tau_min < nlsP->gamma1_i) && (tau_max >= 1.0)) { 203 tau = nlsP->gamma1_i; 204 } else if ((tau_1 >= nlsP->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < nlsP->gamma1_i) || (tau_2 >= 1.0))) { 205 tau = tau_1; 206 } else if ((tau_2 >= nlsP->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < nlsP->gamma1_i) || (tau_2 >= 1.0))) { 207 tau = tau_2; 208 } else { 209 tau = tau_max; 210 } 211 } 212 } 213 tao->trust = tau * tao->trust; 214 } 215 216 if (fmin < f) { 217 f = fmin; 218 PetscCall(VecAXPY(tao->solution, sigma, tao->gradient)); 219 PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient)); 220 221 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm)); 222 PetscCheck(!PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute gradient generated Inf or NaN"); 223 needH = 1; 224 225 PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its)); 226 PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step)); 227 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 228 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS); 229 } 230 } 231 tao->trust = PetscMax(tao->trust, max_radius); 232 233 /* Modify the radius if it is too large or small */ 234 tao->trust = PetscMax(tao->trust, nlsP->min_radius); 235 tao->trust = PetscMin(tao->trust, nlsP->max_radius); 236 break; 237 238 default: 239 /* Norm of the first direction will initialize radius */ 240 tao->trust = 0.0; 241 break; 242 } 243 } 244 245 /* Set counter for gradient/reset steps*/ 246 nlsP->newt = 0; 247 nlsP->bfgs = 0; 248 nlsP->grad = 0; 249 250 /* Have not converged; continue with Newton method */ 251 while (tao->reason == TAO_CONTINUE_ITERATING) { 252 /* Call general purpose update function */ 253 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update); 254 ++tao->niter; 255 tao->ksp_its = 0; 256 257 /* Compute the Hessian */ 258 if (needH) PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre)); 259 260 /* Shift the Hessian matrix */ 261 if (pert > 0) { 262 PetscCall(MatShift(tao->hessian, pert)); 263 if (tao->hessian != tao->hessian_pre) PetscCall(MatShift(tao->hessian_pre, pert)); 264 } 265 266 if (nlsP->bfgs_pre) { 267 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 268 ++bfgsUpdates; 269 } 270 271 /* Solve the Newton system of equations */ 272 PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre)); 273 if (is_nash || is_stcg || is_gltr) { 274 PetscCall(KSPCGSetRadius(tao->ksp, nlsP->max_radius)); 275 PetscCall(KSPSolve(tao->ksp, tao->gradient, nlsP->D)); 276 PetscCall(KSPGetIterationNumber(tao->ksp, &kspits)); 277 tao->ksp_its += kspits; 278 tao->ksp_tot_its += kspits; 279 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 280 281 if (0.0 == tao->trust) { 282 /* Radius was uninitialized; use the norm of the direction */ 283 if (norm_d > 0.0) { 284 tao->trust = norm_d; 285 286 /* Modify the radius if it is too large or small */ 287 tao->trust = PetscMax(tao->trust, nlsP->min_radius); 288 tao->trust = PetscMin(tao->trust, nlsP->max_radius); 289 } else { 290 /* The direction was bad; set radius to default value and re-solve 291 the trust-region subproblem to get a direction */ 292 tao->trust = tao->trust0; 293 294 /* Modify the radius if it is too large or small */ 295 tao->trust = PetscMax(tao->trust, nlsP->min_radius); 296 tao->trust = PetscMin(tao->trust, nlsP->max_radius); 297 298 PetscCall(KSPCGSetRadius(tao->ksp, nlsP->max_radius)); 299 PetscCall(KSPSolve(tao->ksp, tao->gradient, nlsP->D)); 300 PetscCall(KSPGetIterationNumber(tao->ksp, &kspits)); 301 tao->ksp_its += kspits; 302 tao->ksp_tot_its += kspits; 303 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 304 305 PetscCheck(norm_d != 0.0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Initial direction zero"); 306 } 307 } 308 } else { 309 PetscCall(KSPSolve(tao->ksp, tao->gradient, nlsP->D)); 310 PetscCall(KSPGetIterationNumber(tao->ksp, &kspits)); 311 tao->ksp_its += kspits; 312 tao->ksp_tot_its += kspits; 313 } 314 PetscCall(VecScale(nlsP->D, -1.0)); 315 PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason)); 316 if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (nlsP->bfgs_pre)) { 317 /* Preconditioner is numerically indefinite; reset the 318 approximate if using BFGS preconditioning. */ 319 PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE)); 320 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 321 bfgsUpdates = 1; 322 } 323 324 if (KSP_CONVERGED_ATOL == ksp_reason) { 325 ++nlsP->ksp_atol; 326 } else if (KSP_CONVERGED_RTOL == ksp_reason) { 327 ++nlsP->ksp_rtol; 328 } else if (KSP_CONVERGED_STEP_LENGTH == ksp_reason) { 329 ++nlsP->ksp_ctol; 330 } else if (KSP_CONVERGED_NEG_CURVE == ksp_reason) { 331 ++nlsP->ksp_negc; 332 } else if (KSP_DIVERGED_DTOL == ksp_reason) { 333 ++nlsP->ksp_dtol; 334 } else if (KSP_DIVERGED_ITS == ksp_reason) { 335 ++nlsP->ksp_iter; 336 } else { 337 ++nlsP->ksp_othr; 338 } 339 340 /* Check for success (descent direction) */ 341 PetscCall(VecDot(nlsP->D, tao->gradient, &gdx)); 342 if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) { 343 /* Newton step is not descent or direction produced Inf or NaN 344 Update the perturbation for next time */ 345 if (pert <= 0.0) { 346 /* Initialize the perturbation */ 347 pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm)); 348 if (is_gltr) { 349 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 350 pert = PetscMax(pert, -e_min); 351 } 352 } else { 353 /* Increase the perturbation */ 354 pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm)); 355 } 356 357 if (!nlsP->bfgs_pre) { 358 /* We don't have the bfgs matrix around and updated 359 Must use gradient direction in this case */ 360 PetscCall(VecCopy(tao->gradient, nlsP->D)); 361 PetscCall(VecScale(nlsP->D, -1.0)); 362 ++nlsP->grad; 363 stepType = NLS_GRADIENT; 364 } else { 365 /* Attempt to use the BFGS direction */ 366 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 367 PetscCall(VecScale(nlsP->D, -1.0)); 368 369 /* Check for success (descent direction) */ 370 PetscCall(VecDot(tao->gradient, nlsP->D, &gdx)); 371 if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) { 372 /* BFGS direction is not descent or direction produced not a number 373 We can assert bfgsUpdates > 1 in this case because 374 the first solve produces the scaled gradient direction, 375 which is guaranteed to be descent */ 376 377 /* Use steepest descent direction (scaled) */ 378 PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE)); 379 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 380 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 381 PetscCall(VecScale(nlsP->D, -1.0)); 382 383 bfgsUpdates = 1; 384 ++nlsP->grad; 385 stepType = NLS_GRADIENT; 386 } else { 387 PetscCall(MatLMVMGetUpdateCount(nlsP->M, &bfgsUpdates)); 388 if (1 == bfgsUpdates) { 389 /* The first BFGS direction is always the scaled gradient */ 390 ++nlsP->grad; 391 stepType = NLS_GRADIENT; 392 } else { 393 ++nlsP->bfgs; 394 stepType = NLS_BFGS; 395 } 396 } 397 } 398 } else { 399 /* Computed Newton step is descent */ 400 switch (ksp_reason) { 401 case KSP_DIVERGED_NANORINF: 402 case KSP_DIVERGED_BREAKDOWN: 403 case KSP_DIVERGED_INDEFINITE_MAT: 404 case KSP_DIVERGED_INDEFINITE_PC: 405 case KSP_CONVERGED_NEG_CURVE: 406 /* Matrix or preconditioner is indefinite; increase perturbation */ 407 if (pert <= 0.0) { 408 /* Initialize the perturbation */ 409 pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm)); 410 if (is_gltr) { 411 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 412 pert = PetscMax(pert, -e_min); 413 } 414 } else { 415 /* Increase the perturbation */ 416 pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm)); 417 } 418 break; 419 420 default: 421 /* Newton step computation is good; decrease perturbation */ 422 pert = PetscMin(nlsP->psfac * pert, nlsP->pmsfac * gnorm); 423 if (pert < nlsP->pmin) pert = 0.0; 424 break; 425 } 426 427 ++nlsP->newt; 428 stepType = NLS_NEWTON; 429 } 430 431 /* Perform the linesearch */ 432 fold = f; 433 PetscCall(VecCopy(tao->solution, nlsP->Xold)); 434 PetscCall(VecCopy(tao->gradient, nlsP->Gold)); 435 436 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason)); 437 PetscCall(TaoAddLineSearchCounts(tao)); 438 439 while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NLS_GRADIENT) { /* Linesearch failed */ 440 f = fold; 441 PetscCall(VecCopy(nlsP->Xold, tao->solution)); 442 PetscCall(VecCopy(nlsP->Gold, tao->gradient)); 443 444 switch (stepType) { 445 case NLS_NEWTON: 446 /* Failed to obtain acceptable iterate with Newton 1step 447 Update the perturbation for next time */ 448 if (pert <= 0.0) { 449 /* Initialize the perturbation */ 450 pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm)); 451 if (is_gltr) { 452 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 453 pert = PetscMax(pert, -e_min); 454 } 455 } else { 456 /* Increase the perturbation */ 457 pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm)); 458 } 459 460 if (!nlsP->bfgs_pre) { 461 /* We don't have the bfgs matrix around and being updated 462 Must use gradient direction in this case */ 463 PetscCall(VecCopy(tao->gradient, nlsP->D)); 464 ++nlsP->grad; 465 stepType = NLS_GRADIENT; 466 } else { 467 /* Attempt to use the BFGS direction */ 468 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 469 /* Check for success (descent direction) */ 470 PetscCall(VecDot(tao->solution, nlsP->D, &gdx)); 471 if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) { 472 /* BFGS direction is not descent or direction produced not a number 473 We can assert bfgsUpdates > 1 in this case 474 Use steepest descent direction (scaled) */ 475 PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE)); 476 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 477 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 478 479 bfgsUpdates = 1; 480 ++nlsP->grad; 481 stepType = NLS_GRADIENT; 482 } else { 483 if (1 == bfgsUpdates) { 484 /* The first BFGS direction is always the scaled gradient */ 485 ++nlsP->grad; 486 stepType = NLS_GRADIENT; 487 } else { 488 ++nlsP->bfgs; 489 stepType = NLS_BFGS; 490 } 491 } 492 } 493 break; 494 495 case NLS_BFGS: 496 /* Can only enter if pc_type == NLS_PC_BFGS 497 Failed to obtain acceptable iterate with BFGS step 498 Attempt to use the scaled gradient direction */ 499 PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE)); 500 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 501 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 502 503 bfgsUpdates = 1; 504 ++nlsP->grad; 505 stepType = NLS_GRADIENT; 506 break; 507 } 508 PetscCall(VecScale(nlsP->D, -1.0)); 509 510 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason)); 511 PetscCall(TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full)); 512 PetscCall(TaoAddLineSearchCounts(tao)); 513 } 514 515 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 516 /* Failed to find an improving point */ 517 f = fold; 518 PetscCall(VecCopy(nlsP->Xold, tao->solution)); 519 PetscCall(VecCopy(nlsP->Gold, tao->gradient)); 520 step = 0.0; 521 tao->reason = TAO_DIVERGED_LS_FAILURE; 522 break; 523 } 524 525 /* Update trust region radius */ 526 if (is_nash || is_stcg || is_gltr) { 527 switch (nlsP->update_type) { 528 case NLS_UPDATE_STEP: 529 if (stepType == NLS_NEWTON) { 530 if (step < nlsP->nu1) { 531 /* Very bad step taken; reduce radius */ 532 tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust); 533 } else if (step < nlsP->nu2) { 534 /* Reasonably bad step taken; reduce radius */ 535 tao->trust = nlsP->omega2 * PetscMin(norm_d, tao->trust); 536 } else if (step < nlsP->nu3) { 537 /* Reasonable step was taken; leave radius alone */ 538 if (nlsP->omega3 < 1.0) { 539 tao->trust = nlsP->omega3 * PetscMin(norm_d, tao->trust); 540 } else if (nlsP->omega3 > 1.0) { 541 tao->trust = PetscMax(nlsP->omega3 * norm_d, tao->trust); 542 } 543 } else if (step < nlsP->nu4) { 544 /* Full step taken; increase the radius */ 545 tao->trust = PetscMax(nlsP->omega4 * norm_d, tao->trust); 546 } else { 547 /* More than full step taken; increase the radius */ 548 tao->trust = PetscMax(nlsP->omega5 * norm_d, tao->trust); 549 } 550 } else { 551 /* Newton step was not good; reduce the radius */ 552 tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust); 553 } 554 break; 555 556 case NLS_UPDATE_REDUCTION: 557 if (stepType == NLS_NEWTON) { 558 /* Get predicted reduction */ 559 560 PetscCall(KSPCGGetObjFcn(tao->ksp, &prered)); 561 if (prered >= 0.0) { 562 /* The predicted reduction has the wrong sign. This cannot */ 563 /* happen in infinite precision arithmetic. Step should */ 564 /* be rejected! */ 565 tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d); 566 } else { 567 if (PetscIsInfOrNanReal(f_full)) { 568 tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d); 569 } else { 570 /* Compute and actual reduction */ 571 actred = fold - f_full; 572 prered = -prered; 573 if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) { 574 kappa = 1.0; 575 } else { 576 kappa = actred / prered; 577 } 578 579 /* Accept of reject the step and update radius */ 580 if (kappa < nlsP->eta1) { 581 /* Very bad step */ 582 tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d); 583 } else if (kappa < nlsP->eta2) { 584 /* Marginal bad step */ 585 tao->trust = nlsP->alpha2 * PetscMin(tao->trust, norm_d); 586 } else if (kappa < nlsP->eta3) { 587 /* Reasonable step */ 588 if (nlsP->alpha3 < 1.0) { 589 tao->trust = nlsP->alpha3 * PetscMin(norm_d, tao->trust); 590 } else if (nlsP->alpha3 > 1.0) { 591 tao->trust = PetscMax(nlsP->alpha3 * norm_d, tao->trust); 592 } 593 } else if (kappa < nlsP->eta4) { 594 /* Good step */ 595 tao->trust = PetscMax(nlsP->alpha4 * norm_d, tao->trust); 596 } else { 597 /* Very good step */ 598 tao->trust = PetscMax(nlsP->alpha5 * norm_d, tao->trust); 599 } 600 } 601 } 602 } else { 603 /* Newton step was not good; reduce the radius */ 604 tao->trust = nlsP->alpha1 * PetscMin(norm_d, tao->trust); 605 } 606 break; 607 608 default: 609 if (stepType == NLS_NEWTON) { 610 PetscCall(KSPCGGetObjFcn(tao->ksp, &prered)); 611 if (prered >= 0.0) { 612 /* The predicted reduction has the wrong sign. This cannot */ 613 /* happen in infinite precision arithmetic. Step should */ 614 /* be rejected! */ 615 tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d); 616 } else { 617 if (PetscIsInfOrNanReal(f_full)) { 618 tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d); 619 } else { 620 actred = fold - f_full; 621 prered = -prered; 622 if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) { 623 kappa = 1.0; 624 } else { 625 kappa = actred / prered; 626 } 627 628 tau_1 = nlsP->theta * gdx / (nlsP->theta * gdx - (1.0 - nlsP->theta) * prered + actred); 629 tau_2 = nlsP->theta * gdx / (nlsP->theta * gdx + (1.0 + nlsP->theta) * prered - actred); 630 tau_min = PetscMin(tau_1, tau_2); 631 tau_max = PetscMax(tau_1, tau_2); 632 633 if (kappa >= 1.0 - nlsP->mu1) { 634 /* Great agreement */ 635 if (tau_max < 1.0) { 636 tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d); 637 } else if (tau_max > nlsP->gamma4) { 638 tao->trust = PetscMax(tao->trust, nlsP->gamma4 * norm_d); 639 } else { 640 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 641 } 642 } else if (kappa >= 1.0 - nlsP->mu2) { 643 /* Good agreement */ 644 645 if (tau_max < nlsP->gamma2) { 646 tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d); 647 } else if (tau_max > nlsP->gamma3) { 648 tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d); 649 } else if (tau_max < 1.0) { 650 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 651 } else { 652 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 653 } 654 } else { 655 /* Not good agreement */ 656 if (tau_min > 1.0) { 657 tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d); 658 } else if (tau_max < nlsP->gamma1) { 659 tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d); 660 } else if ((tau_min < nlsP->gamma1) && (tau_max >= 1.0)) { 661 tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d); 662 } else if ((tau_1 >= nlsP->gamma1) && (tau_1 < 1.0) && ((tau_2 < nlsP->gamma1) || (tau_2 >= 1.0))) { 663 tao->trust = tau_1 * PetscMin(tao->trust, norm_d); 664 } else if ((tau_2 >= nlsP->gamma1) && (tau_2 < 1.0) && ((tau_1 < nlsP->gamma1) || (tau_2 >= 1.0))) { 665 tao->trust = tau_2 * PetscMin(tao->trust, norm_d); 666 } else { 667 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 668 } 669 } 670 } 671 } 672 } else { 673 /* Newton step was not good; reduce the radius */ 674 tao->trust = nlsP->gamma1 * PetscMin(norm_d, tao->trust); 675 } 676 break; 677 } 678 679 /* The radius may have been increased; modify if it is too large */ 680 tao->trust = PetscMin(tao->trust, nlsP->max_radius); 681 } 682 683 /* Check for termination */ 684 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm)); 685 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number"); 686 needH = 1; 687 PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its)); 688 PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step)); 689 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 690 } 691 PetscFunctionReturn(PETSC_SUCCESS); 692 } 693 694 /* ---------------------------------------------------------- */ 695 static PetscErrorCode TaoSetUp_NLS(Tao tao) 696 { 697 TAO_NLS *nlsP = (TAO_NLS *)tao->data; 698 699 PetscFunctionBegin; 700 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 701 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 702 if (!nlsP->W) PetscCall(VecDuplicate(tao->solution, &nlsP->W)); 703 if (!nlsP->D) PetscCall(VecDuplicate(tao->solution, &nlsP->D)); 704 if (!nlsP->Xold) PetscCall(VecDuplicate(tao->solution, &nlsP->Xold)); 705 if (!nlsP->Gold) PetscCall(VecDuplicate(tao->solution, &nlsP->Gold)); 706 nlsP->M = NULL; 707 nlsP->bfgs_pre = NULL; 708 PetscFunctionReturn(PETSC_SUCCESS); 709 } 710 711 /*------------------------------------------------------------*/ 712 static PetscErrorCode TaoDestroy_NLS(Tao tao) 713 { 714 TAO_NLS *nlsP = (TAO_NLS *)tao->data; 715 716 PetscFunctionBegin; 717 if (tao->setupcalled) { 718 PetscCall(VecDestroy(&nlsP->D)); 719 PetscCall(VecDestroy(&nlsP->W)); 720 PetscCall(VecDestroy(&nlsP->Xold)); 721 PetscCall(VecDestroy(&nlsP->Gold)); 722 } 723 PetscCall(KSPDestroy(&tao->ksp)); 724 PetscCall(PetscFree(tao->data)); 725 PetscFunctionReturn(PETSC_SUCCESS); 726 } 727 728 /*------------------------------------------------------------*/ 729 static PetscErrorCode TaoSetFromOptions_NLS(Tao tao, PetscOptionItems *PetscOptionsObject) 730 { 731 TAO_NLS *nlsP = (TAO_NLS *)tao->data; 732 733 PetscFunctionBegin; 734 PetscOptionsHeadBegin(PetscOptionsObject, "Newton line search method for unconstrained optimization"); 735 PetscCall(PetscOptionsEList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[nlsP->init_type], &nlsP->init_type, NULL)); 736 PetscCall(PetscOptionsEList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[nlsP->update_type], &nlsP->update_type, NULL)); 737 PetscCall(PetscOptionsReal("-tao_nls_sval", "perturbation starting value", "", nlsP->sval, &nlsP->sval, NULL)); 738 PetscCall(PetscOptionsReal("-tao_nls_imin", "minimum initial perturbation", "", nlsP->imin, &nlsP->imin, NULL)); 739 PetscCall(PetscOptionsReal("-tao_nls_imax", "maximum initial perturbation", "", nlsP->imax, &nlsP->imax, NULL)); 740 PetscCall(PetscOptionsReal("-tao_nls_imfac", "initial merit factor", "", nlsP->imfac, &nlsP->imfac, NULL)); 741 PetscCall(PetscOptionsReal("-tao_nls_pmin", "minimum perturbation", "", nlsP->pmin, &nlsP->pmin, NULL)); 742 PetscCall(PetscOptionsReal("-tao_nls_pmax", "maximum perturbation", "", nlsP->pmax, &nlsP->pmax, NULL)); 743 PetscCall(PetscOptionsReal("-tao_nls_pgfac", "growth factor", "", nlsP->pgfac, &nlsP->pgfac, NULL)); 744 PetscCall(PetscOptionsReal("-tao_nls_psfac", "shrink factor", "", nlsP->psfac, &nlsP->psfac, NULL)); 745 PetscCall(PetscOptionsReal("-tao_nls_pmgfac", "merit growth factor", "", nlsP->pmgfac, &nlsP->pmgfac, NULL)); 746 PetscCall(PetscOptionsReal("-tao_nls_pmsfac", "merit shrink factor", "", nlsP->pmsfac, &nlsP->pmsfac, NULL)); 747 PetscCall(PetscOptionsReal("-tao_nls_eta1", "poor steplength; reduce radius", "", nlsP->eta1, &nlsP->eta1, NULL)); 748 PetscCall(PetscOptionsReal("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", nlsP->eta2, &nlsP->eta2, NULL)); 749 PetscCall(PetscOptionsReal("-tao_nls_eta3", "good steplength; increase radius", "", nlsP->eta3, &nlsP->eta3, NULL)); 750 PetscCall(PetscOptionsReal("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", nlsP->eta4, &nlsP->eta4, NULL)); 751 PetscCall(PetscOptionsReal("-tao_nls_alpha1", "", "", nlsP->alpha1, &nlsP->alpha1, NULL)); 752 PetscCall(PetscOptionsReal("-tao_nls_alpha2", "", "", nlsP->alpha2, &nlsP->alpha2, NULL)); 753 PetscCall(PetscOptionsReal("-tao_nls_alpha3", "", "", nlsP->alpha3, &nlsP->alpha3, NULL)); 754 PetscCall(PetscOptionsReal("-tao_nls_alpha4", "", "", nlsP->alpha4, &nlsP->alpha4, NULL)); 755 PetscCall(PetscOptionsReal("-tao_nls_alpha5", "", "", nlsP->alpha5, &nlsP->alpha5, NULL)); 756 PetscCall(PetscOptionsReal("-tao_nls_nu1", "poor steplength; reduce radius", "", nlsP->nu1, &nlsP->nu1, NULL)); 757 PetscCall(PetscOptionsReal("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", nlsP->nu2, &nlsP->nu2, NULL)); 758 PetscCall(PetscOptionsReal("-tao_nls_nu3", "good steplength; increase radius", "", nlsP->nu3, &nlsP->nu3, NULL)); 759 PetscCall(PetscOptionsReal("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", nlsP->nu4, &nlsP->nu4, NULL)); 760 PetscCall(PetscOptionsReal("-tao_nls_omega1", "", "", nlsP->omega1, &nlsP->omega1, NULL)); 761 PetscCall(PetscOptionsReal("-tao_nls_omega2", "", "", nlsP->omega2, &nlsP->omega2, NULL)); 762 PetscCall(PetscOptionsReal("-tao_nls_omega3", "", "", nlsP->omega3, &nlsP->omega3, NULL)); 763 PetscCall(PetscOptionsReal("-tao_nls_omega4", "", "", nlsP->omega4, &nlsP->omega4, NULL)); 764 PetscCall(PetscOptionsReal("-tao_nls_omega5", "", "", nlsP->omega5, &nlsP->omega5, NULL)); 765 PetscCall(PetscOptionsReal("-tao_nls_mu1_i", "", "", nlsP->mu1_i, &nlsP->mu1_i, NULL)); 766 PetscCall(PetscOptionsReal("-tao_nls_mu2_i", "", "", nlsP->mu2_i, &nlsP->mu2_i, NULL)); 767 PetscCall(PetscOptionsReal("-tao_nls_gamma1_i", "", "", nlsP->gamma1_i, &nlsP->gamma1_i, NULL)); 768 PetscCall(PetscOptionsReal("-tao_nls_gamma2_i", "", "", nlsP->gamma2_i, &nlsP->gamma2_i, NULL)); 769 PetscCall(PetscOptionsReal("-tao_nls_gamma3_i", "", "", nlsP->gamma3_i, &nlsP->gamma3_i, NULL)); 770 PetscCall(PetscOptionsReal("-tao_nls_gamma4_i", "", "", nlsP->gamma4_i, &nlsP->gamma4_i, NULL)); 771 PetscCall(PetscOptionsReal("-tao_nls_theta_i", "", "", nlsP->theta_i, &nlsP->theta_i, NULL)); 772 PetscCall(PetscOptionsReal("-tao_nls_mu1", "", "", nlsP->mu1, &nlsP->mu1, NULL)); 773 PetscCall(PetscOptionsReal("-tao_nls_mu2", "", "", nlsP->mu2, &nlsP->mu2, NULL)); 774 PetscCall(PetscOptionsReal("-tao_nls_gamma1", "", "", nlsP->gamma1, &nlsP->gamma1, NULL)); 775 PetscCall(PetscOptionsReal("-tao_nls_gamma2", "", "", nlsP->gamma2, &nlsP->gamma2, NULL)); 776 PetscCall(PetscOptionsReal("-tao_nls_gamma3", "", "", nlsP->gamma3, &nlsP->gamma3, NULL)); 777 PetscCall(PetscOptionsReal("-tao_nls_gamma4", "", "", nlsP->gamma4, &nlsP->gamma4, NULL)); 778 PetscCall(PetscOptionsReal("-tao_nls_theta", "", "", nlsP->theta, &nlsP->theta, NULL)); 779 PetscCall(PetscOptionsReal("-tao_nls_min_radius", "lower bound on initial radius", "", nlsP->min_radius, &nlsP->min_radius, NULL)); 780 PetscCall(PetscOptionsReal("-tao_nls_max_radius", "upper bound on radius", "", nlsP->max_radius, &nlsP->max_radius, NULL)); 781 PetscCall(PetscOptionsReal("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", nlsP->epsilon, &nlsP->epsilon, NULL)); 782 PetscOptionsHeadEnd(); 783 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 784 PetscCall(KSPSetFromOptions(tao->ksp)); 785 PetscFunctionReturn(PETSC_SUCCESS); 786 } 787 788 /*------------------------------------------------------------*/ 789 static PetscErrorCode TaoView_NLS(Tao tao, PetscViewer viewer) 790 { 791 TAO_NLS *nlsP = (TAO_NLS *)tao->data; 792 PetscBool isascii; 793 794 PetscFunctionBegin; 795 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 796 if (isascii) { 797 PetscCall(PetscViewerASCIIPushTab(viewer)); 798 PetscCall(PetscViewerASCIIPrintf(viewer, "Newton steps: %" PetscInt_FMT "\n", nlsP->newt)); 799 PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS steps: %" PetscInt_FMT "\n", nlsP->bfgs)); 800 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", nlsP->grad)); 801 802 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp atol: %" PetscInt_FMT "\n", nlsP->ksp_atol)); 803 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %" PetscInt_FMT "\n", nlsP->ksp_rtol)); 804 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %" PetscInt_FMT "\n", nlsP->ksp_ctol)); 805 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp negc: %" PetscInt_FMT "\n", nlsP->ksp_negc)); 806 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %" PetscInt_FMT "\n", nlsP->ksp_dtol)); 807 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp iter: %" PetscInt_FMT "\n", nlsP->ksp_iter)); 808 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp othr: %" PetscInt_FMT "\n", nlsP->ksp_othr)); 809 PetscCall(PetscViewerASCIIPopTab(viewer)); 810 } 811 PetscFunctionReturn(PETSC_SUCCESS); 812 } 813 814 /* ---------------------------------------------------------- */ 815 /*MC 816 TAONLS - Newton's method with linesearch for unconstrained minimization. 817 At each iteration, the Newton line search method solves the symmetric 818 system of equations to obtain the step direction dk: 819 Hk dk = -gk 820 a More-Thuente line search is applied on the direction dk to approximately 821 solve 822 min_t f(xk + t d_k) 823 824 Options Database Keys: 825 + -tao_nls_init_type - "constant","direction","interpolation" 826 . -tao_nls_update_type - "step","direction","interpolation" 827 . -tao_nls_sval - perturbation starting value 828 . -tao_nls_imin - minimum initial perturbation 829 . -tao_nls_imax - maximum initial perturbation 830 . -tao_nls_pmin - minimum perturbation 831 . -tao_nls_pmax - maximum perturbation 832 . -tao_nls_pgfac - growth factor 833 . -tao_nls_psfac - shrink factor 834 . -tao_nls_imfac - initial merit factor 835 . -tao_nls_pmgfac - merit growth factor 836 . -tao_nls_pmsfac - merit shrink factor 837 . -tao_nls_eta1 - poor steplength; reduce radius 838 . -tao_nls_eta2 - reasonable steplength; leave radius 839 . -tao_nls_eta3 - good steplength; increase radius 840 . -tao_nls_eta4 - excellent steplength; greatly increase radius 841 . -tao_nls_alpha1 - alpha1 reduction 842 . -tao_nls_alpha2 - alpha2 reduction 843 . -tao_nls_alpha3 - alpha3 reduction 844 . -tao_nls_alpha4 - alpha4 reduction 845 . -tao_nls_alpha - alpha5 reduction 846 . -tao_nls_mu1 - mu1 interpolation update 847 . -tao_nls_mu2 - mu2 interpolation update 848 . -tao_nls_gamma1 - gamma1 interpolation update 849 . -tao_nls_gamma2 - gamma2 interpolation update 850 . -tao_nls_gamma3 - gamma3 interpolation update 851 . -tao_nls_gamma4 - gamma4 interpolation update 852 . -tao_nls_theta - theta interpolation update 853 . -tao_nls_omega1 - omega1 step update 854 . -tao_nls_omega2 - omega2 step update 855 . -tao_nls_omega3 - omega3 step update 856 . -tao_nls_omega4 - omega4 step update 857 . -tao_nls_omega5 - omega5 step update 858 . -tao_nls_mu1_i - mu1 interpolation init factor 859 . -tao_nls_mu2_i - mu2 interpolation init factor 860 . -tao_nls_gamma1_i - gamma1 interpolation init factor 861 . -tao_nls_gamma2_i - gamma2 interpolation init factor 862 . -tao_nls_gamma3_i - gamma3 interpolation init factor 863 . -tao_nls_gamma4_i - gamma4 interpolation init factor 864 - -tao_nls_theta_i - theta interpolation init factor 865 866 Level: beginner 867 M*/ 868 869 PETSC_EXTERN PetscErrorCode TaoCreate_NLS(Tao tao) 870 { 871 TAO_NLS *nlsP; 872 const char *morethuente_type = TAOLINESEARCHMT; 873 874 PetscFunctionBegin; 875 PetscCall(PetscNew(&nlsP)); 876 877 tao->ops->setup = TaoSetUp_NLS; 878 tao->ops->solve = TaoSolve_NLS; 879 tao->ops->view = TaoView_NLS; 880 tao->ops->setfromoptions = TaoSetFromOptions_NLS; 881 tao->ops->destroy = TaoDestroy_NLS; 882 883 /* Override default settings (unless already changed) */ 884 if (!tao->max_it_changed) tao->max_it = 50; 885 if (!tao->trust0_changed) tao->trust0 = 100.0; 886 887 tao->data = (void *)nlsP; 888 889 nlsP->sval = 0.0; 890 nlsP->imin = 1.0e-4; 891 nlsP->imax = 1.0e+2; 892 nlsP->imfac = 1.0e-1; 893 894 nlsP->pmin = 1.0e-12; 895 nlsP->pmax = 1.0e+2; 896 nlsP->pgfac = 1.0e+1; 897 nlsP->psfac = 4.0e-1; 898 nlsP->pmgfac = 1.0e-1; 899 nlsP->pmsfac = 1.0e-1; 900 901 /* Default values for trust-region radius update based on steplength */ 902 nlsP->nu1 = 0.25; 903 nlsP->nu2 = 0.50; 904 nlsP->nu3 = 1.00; 905 nlsP->nu4 = 1.25; 906 907 nlsP->omega1 = 0.25; 908 nlsP->omega2 = 0.50; 909 nlsP->omega3 = 1.00; 910 nlsP->omega4 = 2.00; 911 nlsP->omega5 = 4.00; 912 913 /* Default values for trust-region radius update based on reduction */ 914 nlsP->eta1 = 1.0e-4; 915 nlsP->eta2 = 0.25; 916 nlsP->eta3 = 0.50; 917 nlsP->eta4 = 0.90; 918 919 nlsP->alpha1 = 0.25; 920 nlsP->alpha2 = 0.50; 921 nlsP->alpha3 = 1.00; 922 nlsP->alpha4 = 2.00; 923 nlsP->alpha5 = 4.00; 924 925 /* Default values for trust-region radius update based on interpolation */ 926 nlsP->mu1 = 0.10; 927 nlsP->mu2 = 0.50; 928 929 nlsP->gamma1 = 0.25; 930 nlsP->gamma2 = 0.50; 931 nlsP->gamma3 = 2.00; 932 nlsP->gamma4 = 4.00; 933 934 nlsP->theta = 0.05; 935 936 /* Default values for trust region initialization based on interpolation */ 937 nlsP->mu1_i = 0.35; 938 nlsP->mu2_i = 0.50; 939 940 nlsP->gamma1_i = 0.0625; 941 nlsP->gamma2_i = 0.5; 942 nlsP->gamma3_i = 2.0; 943 nlsP->gamma4_i = 5.0; 944 945 nlsP->theta_i = 0.25; 946 947 /* Remaining parameters */ 948 nlsP->min_radius = 1.0e-10; 949 nlsP->max_radius = 1.0e10; 950 nlsP->epsilon = 1.0e-6; 951 952 nlsP->init_type = NLS_INIT_INTERPOLATION; 953 nlsP->update_type = NLS_UPDATE_STEP; 954 955 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 956 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 957 PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type)); 958 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao)); 959 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix)); 960 961 /* Set linear solver to default for symmetric matrices */ 962 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 963 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 964 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix)); 965 PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_nls_")); 966 PetscCall(KSPSetType(tao->ksp, KSPSTCG)); 967 PetscFunctionReturn(PETSC_SUCCESS); 968 } 969