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