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