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 infinity 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 infinity 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 if (tao->ops->update) { 254 PetscUseTypeMethod(tao, update, tao->niter, tao->user_update); 255 PetscCall(TaoComputeObjective(tao, tao->solution, &f)); 256 } 257 ++tao->niter; 258 tao->ksp_its = 0; 259 260 /* Compute the Hessian */ 261 if (needH) PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre)); 262 263 /* Shift the Hessian matrix */ 264 if (pert > 0) { 265 PetscCall(MatShift(tao->hessian, pert)); 266 if (tao->hessian != tao->hessian_pre) PetscCall(MatShift(tao->hessian_pre, pert)); 267 } 268 269 if (nlsP->bfgs_pre) { 270 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 271 ++bfgsUpdates; 272 } 273 274 /* Solve the Newton system of equations */ 275 PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre)); 276 if (is_nash || is_stcg || is_gltr) { 277 PetscCall(KSPCGSetRadius(tao->ksp, nlsP->max_radius)); 278 PetscCall(KSPSolve(tao->ksp, tao->gradient, nlsP->D)); 279 PetscCall(KSPGetIterationNumber(tao->ksp, &kspits)); 280 tao->ksp_its += kspits; 281 tao->ksp_tot_its += kspits; 282 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 283 284 if (0.0 == tao->trust) { 285 /* Radius was uninitialized; use the norm of the direction */ 286 if (norm_d > 0.0) { 287 tao->trust = norm_d; 288 289 /* Modify the radius if it is too large or small */ 290 tao->trust = PetscMax(tao->trust, nlsP->min_radius); 291 tao->trust = PetscMin(tao->trust, nlsP->max_radius); 292 } else { 293 /* The direction was bad; set radius to default value and re-solve 294 the trust-region subproblem to get a direction */ 295 tao->trust = tao->trust0; 296 297 /* Modify the radius if it is too large or small */ 298 tao->trust = PetscMax(tao->trust, nlsP->min_radius); 299 tao->trust = PetscMin(tao->trust, nlsP->max_radius); 300 301 PetscCall(KSPCGSetRadius(tao->ksp, nlsP->max_radius)); 302 PetscCall(KSPSolve(tao->ksp, tao->gradient, nlsP->D)); 303 PetscCall(KSPGetIterationNumber(tao->ksp, &kspits)); 304 tao->ksp_its += kspits; 305 tao->ksp_tot_its += kspits; 306 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 307 308 PetscCheck(norm_d != 0.0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Initial direction zero"); 309 } 310 } 311 } else { 312 PetscCall(KSPSolve(tao->ksp, tao->gradient, nlsP->D)); 313 PetscCall(KSPGetIterationNumber(tao->ksp, &kspits)); 314 tao->ksp_its += kspits; 315 tao->ksp_tot_its += kspits; 316 } 317 PetscCall(VecScale(nlsP->D, -1.0)); 318 PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason)); 319 if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (nlsP->bfgs_pre)) { 320 /* Preconditioner is numerically indefinite; reset the 321 approximate if using BFGS preconditioning. */ 322 PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE)); 323 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 324 bfgsUpdates = 1; 325 } 326 327 if (KSP_CONVERGED_ATOL == ksp_reason) { 328 ++nlsP->ksp_atol; 329 } else if (KSP_CONVERGED_RTOL == ksp_reason) { 330 ++nlsP->ksp_rtol; 331 } else if (KSP_CONVERGED_STEP_LENGTH == ksp_reason) { 332 ++nlsP->ksp_ctol; 333 } else if (KSP_CONVERGED_NEG_CURVE == ksp_reason) { 334 ++nlsP->ksp_negc; 335 } else if (KSP_DIVERGED_DTOL == ksp_reason) { 336 ++nlsP->ksp_dtol; 337 } else if (KSP_DIVERGED_ITS == ksp_reason) { 338 ++nlsP->ksp_iter; 339 } else { 340 ++nlsP->ksp_othr; 341 } 342 343 /* Check for success (descent direction) */ 344 PetscCall(VecDot(nlsP->D, tao->gradient, &gdx)); 345 if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) { 346 /* Newton step is not descent or direction produced infinity or NaN 347 Update the perturbation for next time */ 348 if (pert <= 0.0) { 349 /* Initialize the perturbation */ 350 pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm)); 351 if (is_gltr) { 352 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 353 pert = PetscMax(pert, -e_min); 354 } 355 } else { 356 /* Increase the perturbation */ 357 pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm)); 358 } 359 360 if (!nlsP->bfgs_pre) { 361 /* We don't have the bfgs matrix around and updated 362 Must use gradient direction in this case */ 363 PetscCall(VecCopy(tao->gradient, nlsP->D)); 364 PetscCall(VecScale(nlsP->D, -1.0)); 365 ++nlsP->grad; 366 stepType = NLS_GRADIENT; 367 } else { 368 /* Attempt to use the BFGS direction */ 369 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 370 PetscCall(VecScale(nlsP->D, -1.0)); 371 372 /* Check for success (descent direction) */ 373 PetscCall(VecDot(tao->gradient, nlsP->D, &gdx)); 374 if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) { 375 /* BFGS direction is not descent or direction produced not a number 376 We can assert bfgsUpdates > 1 in this case because 377 the first solve produces the scaled gradient direction, 378 which is guaranteed to be descent */ 379 380 /* Use steepest descent direction (scaled) */ 381 PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE)); 382 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 383 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 384 PetscCall(VecScale(nlsP->D, -1.0)); 385 386 bfgsUpdates = 1; 387 ++nlsP->grad; 388 stepType = NLS_GRADIENT; 389 } else { 390 PetscCall(MatLMVMGetUpdateCount(nlsP->M, &bfgsUpdates)); 391 if (1 == bfgsUpdates) { 392 /* The first BFGS direction is always the scaled gradient */ 393 ++nlsP->grad; 394 stepType = NLS_GRADIENT; 395 } else { 396 ++nlsP->bfgs; 397 stepType = NLS_BFGS; 398 } 399 } 400 } 401 } else { 402 /* Computed Newton step is descent */ 403 switch (ksp_reason) { 404 case KSP_DIVERGED_NANORINF: 405 case KSP_DIVERGED_BREAKDOWN: 406 case KSP_DIVERGED_INDEFINITE_MAT: 407 case KSP_DIVERGED_INDEFINITE_PC: 408 case KSP_CONVERGED_NEG_CURVE: 409 /* Matrix or preconditioner is indefinite; increase perturbation */ 410 if (pert <= 0.0) { 411 /* Initialize the perturbation */ 412 pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm)); 413 if (is_gltr) { 414 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 415 pert = PetscMax(pert, -e_min); 416 } 417 } else { 418 /* Increase the perturbation */ 419 pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm)); 420 } 421 break; 422 423 default: 424 /* Newton step computation is good; decrease perturbation */ 425 pert = PetscMin(nlsP->psfac * pert, nlsP->pmsfac * gnorm); 426 if (pert < nlsP->pmin) pert = 0.0; 427 break; 428 } 429 430 ++nlsP->newt; 431 stepType = NLS_NEWTON; 432 } 433 434 /* Perform the linesearch */ 435 fold = f; 436 PetscCall(VecCopy(tao->solution, nlsP->Xold)); 437 PetscCall(VecCopy(tao->gradient, nlsP->Gold)); 438 439 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason)); 440 PetscCall(TaoAddLineSearchCounts(tao)); 441 442 while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NLS_GRADIENT) { /* Linesearch failed */ 443 f = fold; 444 PetscCall(VecCopy(nlsP->Xold, tao->solution)); 445 PetscCall(VecCopy(nlsP->Gold, tao->gradient)); 446 447 switch (stepType) { 448 case NLS_NEWTON: 449 /* Failed to obtain acceptable iterate with Newton 1step 450 Update the perturbation for next time */ 451 if (pert <= 0.0) { 452 /* Initialize the perturbation */ 453 pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm)); 454 if (is_gltr) { 455 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 456 pert = PetscMax(pert, -e_min); 457 } 458 } else { 459 /* Increase the perturbation */ 460 pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm)); 461 } 462 463 if (!nlsP->bfgs_pre) { 464 /* We don't have the bfgs matrix around and being updated 465 Must use gradient direction in this case */ 466 PetscCall(VecCopy(tao->gradient, nlsP->D)); 467 ++nlsP->grad; 468 stepType = NLS_GRADIENT; 469 } else { 470 /* Attempt to use the BFGS direction */ 471 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 472 /* Check for success (descent direction) */ 473 PetscCall(VecDot(tao->solution, nlsP->D, &gdx)); 474 if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) { 475 /* BFGS direction is not descent or direction produced not a number 476 We can assert bfgsUpdates > 1 in this case 477 Use steepest descent direction (scaled) */ 478 PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE)); 479 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 480 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 481 482 bfgsUpdates = 1; 483 ++nlsP->grad; 484 stepType = NLS_GRADIENT; 485 } else { 486 if (1 == bfgsUpdates) { 487 /* The first BFGS direction is always the scaled gradient */ 488 ++nlsP->grad; 489 stepType = NLS_GRADIENT; 490 } else { 491 ++nlsP->bfgs; 492 stepType = NLS_BFGS; 493 } 494 } 495 } 496 break; 497 498 case NLS_BFGS: 499 /* Can only enter if pc_type == NLS_PC_BFGS 500 Failed to obtain acceptable iterate with BFGS step 501 Attempt to use the scaled gradient direction */ 502 PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE)); 503 PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 504 PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 505 506 bfgsUpdates = 1; 507 ++nlsP->grad; 508 stepType = NLS_GRADIENT; 509 break; 510 } 511 PetscCall(VecScale(nlsP->D, -1.0)); 512 513 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason)); 514 PetscCall(TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full)); 515 PetscCall(TaoAddLineSearchCounts(tao)); 516 } 517 518 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 519 /* Failed to find an improving point */ 520 f = fold; 521 PetscCall(VecCopy(nlsP->Xold, tao->solution)); 522 PetscCall(VecCopy(nlsP->Gold, tao->gradient)); 523 step = 0.0; 524 tao->reason = TAO_DIVERGED_LS_FAILURE; 525 break; 526 } 527 528 /* Update trust region radius */ 529 if (is_nash || is_stcg || is_gltr) { 530 switch (nlsP->update_type) { 531 case NLS_UPDATE_STEP: 532 if (stepType == NLS_NEWTON) { 533 if (step < nlsP->nu1) { 534 /* Very bad step taken; reduce radius */ 535 tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust); 536 } else if (step < nlsP->nu2) { 537 /* Reasonably bad step taken; reduce radius */ 538 tao->trust = nlsP->omega2 * PetscMin(norm_d, tao->trust); 539 } else if (step < nlsP->nu3) { 540 /* Reasonable step was taken; leave radius alone */ 541 if (nlsP->omega3 < 1.0) { 542 tao->trust = nlsP->omega3 * PetscMin(norm_d, tao->trust); 543 } else if (nlsP->omega3 > 1.0) { 544 tao->trust = PetscMax(nlsP->omega3 * norm_d, tao->trust); 545 } 546 } else if (step < nlsP->nu4) { 547 /* Full step taken; increase the radius */ 548 tao->trust = PetscMax(nlsP->omega4 * norm_d, tao->trust); 549 } else { 550 /* More than full step taken; increase the radius */ 551 tao->trust = PetscMax(nlsP->omega5 * norm_d, tao->trust); 552 } 553 } else { 554 /* Newton step was not good; reduce the radius */ 555 tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust); 556 } 557 break; 558 559 case NLS_UPDATE_REDUCTION: 560 if (stepType == NLS_NEWTON) { 561 /* Get predicted reduction */ 562 563 PetscCall(KSPCGGetObjFcn(tao->ksp, &prered)); 564 if (prered >= 0.0) { 565 /* The predicted reduction has the wrong sign. This cannot */ 566 /* happen in infinite precision arithmetic. Step should */ 567 /* be rejected! */ 568 tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d); 569 } else { 570 if (PetscIsInfOrNanReal(f_full)) { 571 tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d); 572 } else { 573 /* Compute and actual reduction */ 574 actred = fold - f_full; 575 prered = -prered; 576 if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) { 577 kappa = 1.0; 578 } else { 579 kappa = actred / prered; 580 } 581 582 /* Accept of reject the step and update radius */ 583 if (kappa < nlsP->eta1) { 584 /* Very bad step */ 585 tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d); 586 } else if (kappa < nlsP->eta2) { 587 /* Marginal bad step */ 588 tao->trust = nlsP->alpha2 * PetscMin(tao->trust, norm_d); 589 } else if (kappa < nlsP->eta3) { 590 /* Reasonable step */ 591 if (nlsP->alpha3 < 1.0) { 592 tao->trust = nlsP->alpha3 * PetscMin(norm_d, tao->trust); 593 } else if (nlsP->alpha3 > 1.0) { 594 tao->trust = PetscMax(nlsP->alpha3 * norm_d, tao->trust); 595 } 596 } else if (kappa < nlsP->eta4) { 597 /* Good step */ 598 tao->trust = PetscMax(nlsP->alpha4 * norm_d, tao->trust); 599 } else { 600 /* Very good step */ 601 tao->trust = PetscMax(nlsP->alpha5 * norm_d, tao->trust); 602 } 603 } 604 } 605 } else { 606 /* Newton step was not good; reduce the radius */ 607 tao->trust = nlsP->alpha1 * PetscMin(norm_d, tao->trust); 608 } 609 break; 610 611 default: 612 if (stepType == NLS_NEWTON) { 613 PetscCall(KSPCGGetObjFcn(tao->ksp, &prered)); 614 if (prered >= 0.0) { 615 /* The predicted reduction has the wrong sign. This cannot */ 616 /* happen in infinite precision arithmetic. Step should */ 617 /* be rejected! */ 618 tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d); 619 } else { 620 if (PetscIsInfOrNanReal(f_full)) { 621 tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d); 622 } else { 623 actred = fold - f_full; 624 prered = -prered; 625 if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) { 626 kappa = 1.0; 627 } else { 628 kappa = actred / prered; 629 } 630 631 tau_1 = nlsP->theta * gdx / (nlsP->theta * gdx - (1.0 - nlsP->theta) * prered + actred); 632 tau_2 = nlsP->theta * gdx / (nlsP->theta * gdx + (1.0 + nlsP->theta) * prered - actred); 633 tau_min = PetscMin(tau_1, tau_2); 634 tau_max = PetscMax(tau_1, tau_2); 635 636 if (kappa >= 1.0 - nlsP->mu1) { 637 /* Great agreement */ 638 if (tau_max < 1.0) { 639 tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d); 640 } else if (tau_max > nlsP->gamma4) { 641 tao->trust = PetscMax(tao->trust, nlsP->gamma4 * norm_d); 642 } else { 643 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 644 } 645 } else if (kappa >= 1.0 - nlsP->mu2) { 646 /* Good agreement */ 647 648 if (tau_max < nlsP->gamma2) { 649 tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d); 650 } else if (tau_max > nlsP->gamma3) { 651 tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d); 652 } else if (tau_max < 1.0) { 653 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 654 } else { 655 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 656 } 657 } else { 658 /* Not good agreement */ 659 if (tau_min > 1.0) { 660 tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d); 661 } else if (tau_max < nlsP->gamma1) { 662 tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d); 663 } else if ((tau_min < nlsP->gamma1) && (tau_max >= 1.0)) { 664 tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d); 665 } else if ((tau_1 >= nlsP->gamma1) && (tau_1 < 1.0) && ((tau_2 < nlsP->gamma1) || (tau_2 >= 1.0))) { 666 tao->trust = tau_1 * PetscMin(tao->trust, norm_d); 667 } else if ((tau_2 >= nlsP->gamma1) && (tau_2 < 1.0) && ((tau_1 < nlsP->gamma1) || (tau_2 >= 1.0))) { 668 tao->trust = tau_2 * PetscMin(tao->trust, norm_d); 669 } else { 670 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 671 } 672 } 673 } 674 } 675 } else { 676 /* Newton step was not good; reduce the radius */ 677 tao->trust = nlsP->gamma1 * PetscMin(norm_d, tao->trust); 678 } 679 break; 680 } 681 682 /* The radius may have been increased; modify if it is too large */ 683 tao->trust = PetscMin(tao->trust, nlsP->max_radius); 684 } 685 686 /* Check for termination */ 687 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm)); 688 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number"); 689 needH = 1; 690 PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its)); 691 PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step)); 692 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 693 } 694 PetscFunctionReturn(PETSC_SUCCESS); 695 } 696 697 /* ---------------------------------------------------------- */ 698 static PetscErrorCode TaoSetUp_NLS(Tao tao) 699 { 700 TAO_NLS *nlsP = (TAO_NLS *)tao->data; 701 702 PetscFunctionBegin; 703 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 704 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 705 if (!nlsP->W) PetscCall(VecDuplicate(tao->solution, &nlsP->W)); 706 if (!nlsP->D) PetscCall(VecDuplicate(tao->solution, &nlsP->D)); 707 if (!nlsP->Xold) PetscCall(VecDuplicate(tao->solution, &nlsP->Xold)); 708 if (!nlsP->Gold) PetscCall(VecDuplicate(tao->solution, &nlsP->Gold)); 709 nlsP->M = NULL; 710 nlsP->bfgs_pre = NULL; 711 PetscFunctionReturn(PETSC_SUCCESS); 712 } 713 714 /*------------------------------------------------------------*/ 715 static PetscErrorCode TaoDestroy_NLS(Tao tao) 716 { 717 TAO_NLS *nlsP = (TAO_NLS *)tao->data; 718 719 PetscFunctionBegin; 720 if (tao->setupcalled) { 721 PetscCall(VecDestroy(&nlsP->D)); 722 PetscCall(VecDestroy(&nlsP->W)); 723 PetscCall(VecDestroy(&nlsP->Xold)); 724 PetscCall(VecDestroy(&nlsP->Gold)); 725 } 726 PetscCall(KSPDestroy(&tao->ksp)); 727 PetscCall(PetscFree(tao->data)); 728 PetscFunctionReturn(PETSC_SUCCESS); 729 } 730 731 /*------------------------------------------------------------*/ 732 static PetscErrorCode TaoSetFromOptions_NLS(Tao tao, PetscOptionItems PetscOptionsObject) 733 { 734 TAO_NLS *nlsP = (TAO_NLS *)tao->data; 735 736 PetscFunctionBegin; 737 PetscOptionsHeadBegin(PetscOptionsObject, "Newton line search method for unconstrained optimization"); 738 PetscCall(PetscOptionsEList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[nlsP->init_type], &nlsP->init_type, NULL)); 739 PetscCall(PetscOptionsEList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[nlsP->update_type], &nlsP->update_type, NULL)); 740 PetscCall(PetscOptionsReal("-tao_nls_sval", "perturbation starting value", "", nlsP->sval, &nlsP->sval, NULL)); 741 PetscCall(PetscOptionsReal("-tao_nls_imin", "minimum initial perturbation", "", nlsP->imin, &nlsP->imin, NULL)); 742 PetscCall(PetscOptionsReal("-tao_nls_imax", "maximum initial perturbation", "", nlsP->imax, &nlsP->imax, NULL)); 743 PetscCall(PetscOptionsReal("-tao_nls_imfac", "initial merit factor", "", nlsP->imfac, &nlsP->imfac, NULL)); 744 PetscCall(PetscOptionsReal("-tao_nls_pmin", "minimum perturbation", "", nlsP->pmin, &nlsP->pmin, NULL)); 745 PetscCall(PetscOptionsReal("-tao_nls_pmax", "maximum perturbation", "", nlsP->pmax, &nlsP->pmax, NULL)); 746 PetscCall(PetscOptionsReal("-tao_nls_pgfac", "growth factor", "", nlsP->pgfac, &nlsP->pgfac, NULL)); 747 PetscCall(PetscOptionsReal("-tao_nls_psfac", "shrink factor", "", nlsP->psfac, &nlsP->psfac, NULL)); 748 PetscCall(PetscOptionsReal("-tao_nls_pmgfac", "merit growth factor", "", nlsP->pmgfac, &nlsP->pmgfac, NULL)); 749 PetscCall(PetscOptionsReal("-tao_nls_pmsfac", "merit shrink factor", "", nlsP->pmsfac, &nlsP->pmsfac, NULL)); 750 PetscCall(PetscOptionsReal("-tao_nls_eta1", "poor steplength; reduce radius", "", nlsP->eta1, &nlsP->eta1, NULL)); 751 PetscCall(PetscOptionsReal("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", nlsP->eta2, &nlsP->eta2, NULL)); 752 PetscCall(PetscOptionsReal("-tao_nls_eta3", "good steplength; increase radius", "", nlsP->eta3, &nlsP->eta3, NULL)); 753 PetscCall(PetscOptionsReal("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", nlsP->eta4, &nlsP->eta4, NULL)); 754 PetscCall(PetscOptionsReal("-tao_nls_alpha1", "", "", nlsP->alpha1, &nlsP->alpha1, NULL)); 755 PetscCall(PetscOptionsReal("-tao_nls_alpha2", "", "", nlsP->alpha2, &nlsP->alpha2, NULL)); 756 PetscCall(PetscOptionsReal("-tao_nls_alpha3", "", "", nlsP->alpha3, &nlsP->alpha3, NULL)); 757 PetscCall(PetscOptionsReal("-tao_nls_alpha4", "", "", nlsP->alpha4, &nlsP->alpha4, NULL)); 758 PetscCall(PetscOptionsReal("-tao_nls_alpha5", "", "", nlsP->alpha5, &nlsP->alpha5, NULL)); 759 PetscCall(PetscOptionsReal("-tao_nls_nu1", "poor steplength; reduce radius", "", nlsP->nu1, &nlsP->nu1, NULL)); 760 PetscCall(PetscOptionsReal("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", nlsP->nu2, &nlsP->nu2, NULL)); 761 PetscCall(PetscOptionsReal("-tao_nls_nu3", "good steplength; increase radius", "", nlsP->nu3, &nlsP->nu3, NULL)); 762 PetscCall(PetscOptionsReal("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", nlsP->nu4, &nlsP->nu4, NULL)); 763 PetscCall(PetscOptionsReal("-tao_nls_omega1", "", "", nlsP->omega1, &nlsP->omega1, NULL)); 764 PetscCall(PetscOptionsReal("-tao_nls_omega2", "", "", nlsP->omega2, &nlsP->omega2, NULL)); 765 PetscCall(PetscOptionsReal("-tao_nls_omega3", "", "", nlsP->omega3, &nlsP->omega3, NULL)); 766 PetscCall(PetscOptionsReal("-tao_nls_omega4", "", "", nlsP->omega4, &nlsP->omega4, NULL)); 767 PetscCall(PetscOptionsReal("-tao_nls_omega5", "", "", nlsP->omega5, &nlsP->omega5, NULL)); 768 PetscCall(PetscOptionsReal("-tao_nls_mu1_i", "", "", nlsP->mu1_i, &nlsP->mu1_i, NULL)); 769 PetscCall(PetscOptionsReal("-tao_nls_mu2_i", "", "", nlsP->mu2_i, &nlsP->mu2_i, NULL)); 770 PetscCall(PetscOptionsReal("-tao_nls_gamma1_i", "", "", nlsP->gamma1_i, &nlsP->gamma1_i, NULL)); 771 PetscCall(PetscOptionsReal("-tao_nls_gamma2_i", "", "", nlsP->gamma2_i, &nlsP->gamma2_i, NULL)); 772 PetscCall(PetscOptionsReal("-tao_nls_gamma3_i", "", "", nlsP->gamma3_i, &nlsP->gamma3_i, NULL)); 773 PetscCall(PetscOptionsReal("-tao_nls_gamma4_i", "", "", nlsP->gamma4_i, &nlsP->gamma4_i, NULL)); 774 PetscCall(PetscOptionsReal("-tao_nls_theta_i", "", "", nlsP->theta_i, &nlsP->theta_i, NULL)); 775 PetscCall(PetscOptionsReal("-tao_nls_mu1", "", "", nlsP->mu1, &nlsP->mu1, NULL)); 776 PetscCall(PetscOptionsReal("-tao_nls_mu2", "", "", nlsP->mu2, &nlsP->mu2, NULL)); 777 PetscCall(PetscOptionsReal("-tao_nls_gamma1", "", "", nlsP->gamma1, &nlsP->gamma1, NULL)); 778 PetscCall(PetscOptionsReal("-tao_nls_gamma2", "", "", nlsP->gamma2, &nlsP->gamma2, NULL)); 779 PetscCall(PetscOptionsReal("-tao_nls_gamma3", "", "", nlsP->gamma3, &nlsP->gamma3, NULL)); 780 PetscCall(PetscOptionsReal("-tao_nls_gamma4", "", "", nlsP->gamma4, &nlsP->gamma4, NULL)); 781 PetscCall(PetscOptionsReal("-tao_nls_theta", "", "", nlsP->theta, &nlsP->theta, NULL)); 782 PetscCall(PetscOptionsReal("-tao_nls_min_radius", "lower bound on initial radius", "", nlsP->min_radius, &nlsP->min_radius, NULL)); 783 PetscCall(PetscOptionsReal("-tao_nls_max_radius", "upper bound on radius", "", nlsP->max_radius, &nlsP->max_radius, NULL)); 784 PetscCall(PetscOptionsReal("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", nlsP->epsilon, &nlsP->epsilon, NULL)); 785 PetscOptionsHeadEnd(); 786 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 787 PetscCall(KSPSetFromOptions(tao->ksp)); 788 PetscFunctionReturn(PETSC_SUCCESS); 789 } 790 791 /*------------------------------------------------------------*/ 792 static PetscErrorCode TaoView_NLS(Tao tao, PetscViewer viewer) 793 { 794 TAO_NLS *nlsP = (TAO_NLS *)tao->data; 795 PetscBool isascii; 796 797 PetscFunctionBegin; 798 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 799 if (isascii) { 800 PetscCall(PetscViewerASCIIPushTab(viewer)); 801 PetscCall(PetscViewerASCIIPrintf(viewer, "Newton steps: %" PetscInt_FMT "\n", nlsP->newt)); 802 PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS steps: %" PetscInt_FMT "\n", nlsP->bfgs)); 803 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", nlsP->grad)); 804 805 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp atol: %" PetscInt_FMT "\n", nlsP->ksp_atol)); 806 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %" PetscInt_FMT "\n", nlsP->ksp_rtol)); 807 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %" PetscInt_FMT "\n", nlsP->ksp_ctol)); 808 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp negc: %" PetscInt_FMT "\n", nlsP->ksp_negc)); 809 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %" PetscInt_FMT "\n", nlsP->ksp_dtol)); 810 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp iter: %" PetscInt_FMT "\n", nlsP->ksp_iter)); 811 PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp othr: %" PetscInt_FMT "\n", nlsP->ksp_othr)); 812 PetscCall(PetscViewerASCIIPopTab(viewer)); 813 } 814 PetscFunctionReturn(PETSC_SUCCESS); 815 } 816 817 /* ---------------------------------------------------------- */ 818 /*MC 819 TAONLS - Newton's method with linesearch for unconstrained minimization. 820 At each iteration, the Newton line search method solves the symmetric 821 system of equations to obtain the step direction dk: 822 Hk dk = -gk 823 a More-Thuente line search is applied on the direction dk to approximately 824 solve 825 min_t f(xk + t d_k) 826 827 Options Database Keys: 828 + -tao_nls_init_type - "constant","direction","interpolation" 829 . -tao_nls_update_type - "step","direction","interpolation" 830 . -tao_nls_sval - perturbation starting value 831 . -tao_nls_imin - minimum initial perturbation 832 . -tao_nls_imax - maximum initial perturbation 833 . -tao_nls_pmin - minimum perturbation 834 . -tao_nls_pmax - maximum perturbation 835 . -tao_nls_pgfac - growth factor 836 . -tao_nls_psfac - shrink factor 837 . -tao_nls_imfac - initial merit factor 838 . -tao_nls_pmgfac - merit growth factor 839 . -tao_nls_pmsfac - merit shrink factor 840 . -tao_nls_eta1 - poor steplength; reduce radius 841 . -tao_nls_eta2 - reasonable steplength; leave radius 842 . -tao_nls_eta3 - good steplength; increase radius 843 . -tao_nls_eta4 - excellent steplength; greatly increase radius 844 . -tao_nls_alpha1 - alpha1 reduction 845 . -tao_nls_alpha2 - alpha2 reduction 846 . -tao_nls_alpha3 - alpha3 reduction 847 . -tao_nls_alpha4 - alpha4 reduction 848 . -tao_nls_alpha - alpha5 reduction 849 . -tao_nls_mu1 - mu1 interpolation update 850 . -tao_nls_mu2 - mu2 interpolation update 851 . -tao_nls_gamma1 - gamma1 interpolation update 852 . -tao_nls_gamma2 - gamma2 interpolation update 853 . -tao_nls_gamma3 - gamma3 interpolation update 854 . -tao_nls_gamma4 - gamma4 interpolation update 855 . -tao_nls_theta - theta interpolation update 856 . -tao_nls_omega1 - omega1 step update 857 . -tao_nls_omega2 - omega2 step update 858 . -tao_nls_omega3 - omega3 step update 859 . -tao_nls_omega4 - omega4 step update 860 . -tao_nls_omega5 - omega5 step update 861 . -tao_nls_mu1_i - mu1 interpolation init factor 862 . -tao_nls_mu2_i - mu2 interpolation init factor 863 . -tao_nls_gamma1_i - gamma1 interpolation init factor 864 . -tao_nls_gamma2_i - gamma2 interpolation init factor 865 . -tao_nls_gamma3_i - gamma3 interpolation init factor 866 . -tao_nls_gamma4_i - gamma4 interpolation init factor 867 - -tao_nls_theta_i - theta interpolation init factor 868 869 Level: beginner 870 M*/ 871 872 PETSC_EXTERN PetscErrorCode TaoCreate_NLS(Tao tao) 873 { 874 TAO_NLS *nlsP; 875 const char *morethuente_type = TAOLINESEARCHMT; 876 877 PetscFunctionBegin; 878 PetscCall(PetscNew(&nlsP)); 879 880 tao->ops->setup = TaoSetUp_NLS; 881 tao->ops->solve = TaoSolve_NLS; 882 tao->ops->view = TaoView_NLS; 883 tao->ops->setfromoptions = TaoSetFromOptions_NLS; 884 tao->ops->destroy = TaoDestroy_NLS; 885 886 /* Override default settings (unless already changed) */ 887 PetscCall(TaoParametersInitialize(tao)); 888 PetscObjectParameterSetDefault(tao, max_it, 50); 889 PetscObjectParameterSetDefault(tao, trust0, 100.0); 890 891 tao->data = (void *)nlsP; 892 893 nlsP->sval = 0.0; 894 nlsP->imin = 1.0e-4; 895 nlsP->imax = 1.0e+2; 896 nlsP->imfac = 1.0e-1; 897 898 nlsP->pmin = 1.0e-12; 899 nlsP->pmax = 1.0e+2; 900 nlsP->pgfac = 1.0e+1; 901 nlsP->psfac = 4.0e-1; 902 nlsP->pmgfac = 1.0e-1; 903 nlsP->pmsfac = 1.0e-1; 904 905 /* Default values for trust-region radius update based on steplength */ 906 nlsP->nu1 = 0.25; 907 nlsP->nu2 = 0.50; 908 nlsP->nu3 = 1.00; 909 nlsP->nu4 = 1.25; 910 911 nlsP->omega1 = 0.25; 912 nlsP->omega2 = 0.50; 913 nlsP->omega3 = 1.00; 914 nlsP->omega4 = 2.00; 915 nlsP->omega5 = 4.00; 916 917 /* Default values for trust-region radius update based on reduction */ 918 nlsP->eta1 = 1.0e-4; 919 nlsP->eta2 = 0.25; 920 nlsP->eta3 = 0.50; 921 nlsP->eta4 = 0.90; 922 923 nlsP->alpha1 = 0.25; 924 nlsP->alpha2 = 0.50; 925 nlsP->alpha3 = 1.00; 926 nlsP->alpha4 = 2.00; 927 nlsP->alpha5 = 4.00; 928 929 /* Default values for trust-region radius update based on interpolation */ 930 nlsP->mu1 = 0.10; 931 nlsP->mu2 = 0.50; 932 933 nlsP->gamma1 = 0.25; 934 nlsP->gamma2 = 0.50; 935 nlsP->gamma3 = 2.00; 936 nlsP->gamma4 = 4.00; 937 938 nlsP->theta = 0.05; 939 940 /* Default values for trust region initialization based on interpolation */ 941 nlsP->mu1_i = 0.35; 942 nlsP->mu2_i = 0.50; 943 944 nlsP->gamma1_i = 0.0625; 945 nlsP->gamma2_i = 0.5; 946 nlsP->gamma3_i = 2.0; 947 nlsP->gamma4_i = 5.0; 948 949 nlsP->theta_i = 0.25; 950 951 /* Remaining parameters */ 952 nlsP->min_radius = 1.0e-10; 953 nlsP->max_radius = 1.0e10; 954 nlsP->epsilon = 1.0e-6; 955 956 nlsP->init_type = NLS_INIT_INTERPOLATION; 957 nlsP->update_type = NLS_UPDATE_STEP; 958 959 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 960 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 961 PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type)); 962 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao)); 963 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix)); 964 965 /* Set linear solver to default for symmetric matrices */ 966 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 967 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 968 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix)); 969 PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_nls_")); 970 PetscCall(KSPSetType(tao->ksp, KSPSTCG)); 971 PetscFunctionReturn(PETSC_SUCCESS); 972 } 973