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 CHKERRQ(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 CHKERRQ(KSPGetType(tao->ksp,&ksp_type)); 84 CHKERRQ(PetscStrcmp(ksp_type,KSPNASH,&is_nash)); 85 CHKERRQ(PetscStrcmp(ksp_type,KSPSTCG,&is_stcg)); 86 CHKERRQ(PetscStrcmp(ksp_type,KSPGLTR,&is_gltr)); 87 88 CHKERRQ(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 CHKERRQ(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient)); 99 CHKERRQ(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 CHKERRQ(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 104 CHKERRQ(TaoMonitor(tao,tao->niter,f,gnorm,0.0,step)); 105 CHKERRQ((*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 CHKERRQ(KSPGetPC(tao->ksp, &pc)); 110 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs)); 111 CHKERRQ(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi)); 112 if (is_bfgs) { 113 nlsP->bfgs_pre = pc; 114 CHKERRQ(PCLMVMGetMatLMVM(nlsP->bfgs_pre, &nlsP->M)); 115 CHKERRQ(VecGetLocalSize(tao->solution, &n)); 116 CHKERRQ(VecGetSize(tao->solution, &N)); 117 CHKERRQ(MatSetSizes(nlsP->M, n, n, N, N)); 118 CHKERRQ(MatLMVMAllocate(nlsP->M, tao->solution, tao->gradient)); 119 CHKERRQ(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 CHKERRQ(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 CHKERRQ(TaoComputeHessian(tao, tao->solution,tao->hessian,tao->hessian_pre)); 143 needH = 0; 144 } 145 146 for (i = 0; i < i_max; ++i) { 147 CHKERRQ(VecCopy(tao->solution,nlsP->W)); 148 CHKERRQ(VecAXPY(nlsP->W,-tao->trust/gnorm,tao->gradient)); 149 CHKERRQ(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 CHKERRQ(MatMult(tao->hessian, tao->gradient, nlsP->D)); 159 CHKERRQ(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 CHKERRQ(VecAXPY(tao->solution,sigma,tao->gradient)); 223 CHKERRQ(TaoComputeGradient(tao,tao->solution,tao->gradient)); 224 225 CHKERRQ(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 CHKERRQ(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 230 CHKERRQ(TaoMonitor(tao,tao->niter,f,gnorm,0.0,step)); 231 CHKERRQ((*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 CHKERRQ((*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 CHKERRQ(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 266 } 267 268 /* Shift the Hessian matrix */ 269 if (pert > 0) { 270 CHKERRQ(MatShift(tao->hessian, pert)); 271 if (tao->hessian != tao->hessian_pre) { 272 CHKERRQ(MatShift(tao->hessian_pre, pert)); 273 } 274 } 275 276 if (nlsP->bfgs_pre) { 277 CHKERRQ(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 278 ++bfgsUpdates; 279 } 280 281 /* Solve the Newton system of equations */ 282 CHKERRQ(KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre)); 283 if (is_nash || is_stcg || is_gltr) { 284 CHKERRQ(KSPCGSetRadius(tao->ksp,nlsP->max_radius)); 285 CHKERRQ(KSPSolve(tao->ksp, tao->gradient, nlsP->D)); 286 CHKERRQ(KSPGetIterationNumber(tao->ksp,&kspits)); 287 tao->ksp_its += kspits; 288 tao->ksp_tot_its += kspits; 289 CHKERRQ(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 CHKERRQ(KSPCGSetRadius(tao->ksp,nlsP->max_radius)); 309 CHKERRQ(KSPSolve(tao->ksp, tao->gradient, nlsP->D)); 310 CHKERRQ(KSPGetIterationNumber(tao->ksp,&kspits)); 311 tao->ksp_its += kspits; 312 tao->ksp_tot_its += kspits; 313 CHKERRQ(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 CHKERRQ(KSPSolve(tao->ksp, tao->gradient, nlsP->D)); 320 CHKERRQ(KSPGetIterationNumber(tao->ksp, &kspits)); 321 tao->ksp_its += kspits; 322 tao->ksp_tot_its += kspits; 323 } 324 CHKERRQ(VecScale(nlsP->D, -1.0)); 325 CHKERRQ(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 CHKERRQ(MatLMVMReset(nlsP->M, PETSC_FALSE)); 330 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecCopy(tao->gradient, nlsP->D)); 371 CHKERRQ(VecScale(nlsP->D, -1.0)); 372 ++nlsP->grad; 373 stepType = NLS_GRADIENT; 374 } else { 375 /* Attempt to use the BFGS direction */ 376 CHKERRQ(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 377 CHKERRQ(VecScale(nlsP->D, -1.0)); 378 379 /* Check for success (descent direction) */ 380 CHKERRQ(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 CHKERRQ(MatLMVMReset(nlsP->M, PETSC_FALSE)); 389 CHKERRQ(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 390 CHKERRQ(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 391 CHKERRQ(VecScale(nlsP->D, -1.0)); 392 393 bfgsUpdates = 1; 394 ++nlsP->grad; 395 stepType = NLS_GRADIENT; 396 } else { 397 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecCopy(tao->solution, nlsP->Xold)); 446 CHKERRQ(VecCopy(tao->gradient, nlsP->Gold)); 447 448 CHKERRQ(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason)); 449 CHKERRQ(TaoAddLineSearchCounts(tao)); 450 451 while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NLS_GRADIENT) { /* Linesearch failed */ 452 f = fold; 453 CHKERRQ(VecCopy(nlsP->Xold, tao->solution)); 454 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecCopy(tao->gradient, nlsP->D)); 476 ++nlsP->grad; 477 stepType = NLS_GRADIENT; 478 } else { 479 /* Attempt to use the BFGS direction */ 480 CHKERRQ(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 481 /* Check for success (descent direction) */ 482 CHKERRQ(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 CHKERRQ(MatLMVMReset(nlsP->M, PETSC_FALSE)); 488 CHKERRQ(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 489 CHKERRQ(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 CHKERRQ(MatLMVMReset(nlsP->M, PETSC_FALSE)); 512 CHKERRQ(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient)); 513 CHKERRQ(MatSolve(nlsP->M, tao->gradient, nlsP->D)); 514 515 bfgsUpdates = 1; 516 ++nlsP->grad; 517 stepType = NLS_GRADIENT; 518 break; 519 } 520 CHKERRQ(VecScale(nlsP->D, -1.0)); 521 522 CHKERRQ(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason)); 523 CHKERRQ(TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full)); 524 CHKERRQ(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 CHKERRQ(VecCopy(nlsP->Xold, tao->solution)); 531 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 701 CHKERRQ(TaoMonitor(tao,tao->niter,f,gnorm,0.0,step)); 702 CHKERRQ((*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) CHKERRQ(VecDuplicate(tao->solution,&tao->gradient)); 714 if (!tao->stepdirection) CHKERRQ(VecDuplicate(tao->solution,&tao->stepdirection)); 715 if (!nlsP->W) CHKERRQ(VecDuplicate(tao->solution,&nlsP->W)); 716 if (!nlsP->D) CHKERRQ(VecDuplicate(tao->solution,&nlsP->D)); 717 if (!nlsP->Xold) CHKERRQ(VecDuplicate(tao->solution,&nlsP->Xold)); 718 if (!nlsP->Gold) CHKERRQ(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 CHKERRQ(VecDestroy(&nlsP->D)); 732 CHKERRQ(VecDestroy(&nlsP->W)); 733 CHKERRQ(VecDestroy(&nlsP->Xold)); 734 CHKERRQ(VecDestroy(&nlsP->Gold)); 735 } 736 CHKERRQ(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 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"Newton line search method for unconstrained optimization")); 747 CHKERRQ(PetscOptionsEList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[nlsP->init_type], &nlsP->init_type, NULL)); 748 CHKERRQ(PetscOptionsEList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[nlsP->update_type], &nlsP->update_type, NULL)); 749 CHKERRQ(PetscOptionsReal("-tao_nls_sval", "perturbation starting value", "", nlsP->sval, &nlsP->sval,NULL)); 750 CHKERRQ(PetscOptionsReal("-tao_nls_imin", "minimum initial perturbation", "", nlsP->imin, &nlsP->imin,NULL)); 751 CHKERRQ(PetscOptionsReal("-tao_nls_imax", "maximum initial perturbation", "", nlsP->imax, &nlsP->imax,NULL)); 752 CHKERRQ(PetscOptionsReal("-tao_nls_imfac", "initial merit factor", "", nlsP->imfac, &nlsP->imfac,NULL)); 753 CHKERRQ(PetscOptionsReal("-tao_nls_pmin", "minimum perturbation", "", nlsP->pmin, &nlsP->pmin,NULL)); 754 CHKERRQ(PetscOptionsReal("-tao_nls_pmax", "maximum perturbation", "", nlsP->pmax, &nlsP->pmax,NULL)); 755 CHKERRQ(PetscOptionsReal("-tao_nls_pgfac", "growth factor", "", nlsP->pgfac, &nlsP->pgfac,NULL)); 756 CHKERRQ(PetscOptionsReal("-tao_nls_psfac", "shrink factor", "", nlsP->psfac, &nlsP->psfac,NULL)); 757 CHKERRQ(PetscOptionsReal("-tao_nls_pmgfac", "merit growth factor", "", nlsP->pmgfac, &nlsP->pmgfac,NULL)); 758 CHKERRQ(PetscOptionsReal("-tao_nls_pmsfac", "merit shrink factor", "", nlsP->pmsfac, &nlsP->pmsfac,NULL)); 759 CHKERRQ(PetscOptionsReal("-tao_nls_eta1", "poor steplength; reduce radius", "", nlsP->eta1, &nlsP->eta1,NULL)); 760 CHKERRQ(PetscOptionsReal("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", nlsP->eta2, &nlsP->eta2,NULL)); 761 CHKERRQ(PetscOptionsReal("-tao_nls_eta3", "good steplength; increase radius", "", nlsP->eta3, &nlsP->eta3,NULL)); 762 CHKERRQ(PetscOptionsReal("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", nlsP->eta4, &nlsP->eta4,NULL)); 763 CHKERRQ(PetscOptionsReal("-tao_nls_alpha1", "", "", nlsP->alpha1, &nlsP->alpha1,NULL)); 764 CHKERRQ(PetscOptionsReal("-tao_nls_alpha2", "", "", nlsP->alpha2, &nlsP->alpha2,NULL)); 765 CHKERRQ(PetscOptionsReal("-tao_nls_alpha3", "", "", nlsP->alpha3, &nlsP->alpha3,NULL)); 766 CHKERRQ(PetscOptionsReal("-tao_nls_alpha4", "", "", nlsP->alpha4, &nlsP->alpha4,NULL)); 767 CHKERRQ(PetscOptionsReal("-tao_nls_alpha5", "", "", nlsP->alpha5, &nlsP->alpha5,NULL)); 768 CHKERRQ(PetscOptionsReal("-tao_nls_nu1", "poor steplength; reduce radius", "", nlsP->nu1, &nlsP->nu1,NULL)); 769 CHKERRQ(PetscOptionsReal("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", nlsP->nu2, &nlsP->nu2,NULL)); 770 CHKERRQ(PetscOptionsReal("-tao_nls_nu3", "good steplength; increase radius", "", nlsP->nu3, &nlsP->nu3,NULL)); 771 CHKERRQ(PetscOptionsReal("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", nlsP->nu4, &nlsP->nu4,NULL)); 772 CHKERRQ(PetscOptionsReal("-tao_nls_omega1", "", "", nlsP->omega1, &nlsP->omega1,NULL)); 773 CHKERRQ(PetscOptionsReal("-tao_nls_omega2", "", "", nlsP->omega2, &nlsP->omega2,NULL)); 774 CHKERRQ(PetscOptionsReal("-tao_nls_omega3", "", "", nlsP->omega3, &nlsP->omega3,NULL)); 775 CHKERRQ(PetscOptionsReal("-tao_nls_omega4", "", "", nlsP->omega4, &nlsP->omega4,NULL)); 776 CHKERRQ(PetscOptionsReal("-tao_nls_omega5", "", "", nlsP->omega5, &nlsP->omega5,NULL)); 777 CHKERRQ(PetscOptionsReal("-tao_nls_mu1_i", "", "", nlsP->mu1_i, &nlsP->mu1_i,NULL)); 778 CHKERRQ(PetscOptionsReal("-tao_nls_mu2_i", "", "", nlsP->mu2_i, &nlsP->mu2_i,NULL)); 779 CHKERRQ(PetscOptionsReal("-tao_nls_gamma1_i", "", "", nlsP->gamma1_i, &nlsP->gamma1_i,NULL)); 780 CHKERRQ(PetscOptionsReal("-tao_nls_gamma2_i", "", "", nlsP->gamma2_i, &nlsP->gamma2_i,NULL)); 781 CHKERRQ(PetscOptionsReal("-tao_nls_gamma3_i", "", "", nlsP->gamma3_i, &nlsP->gamma3_i,NULL)); 782 CHKERRQ(PetscOptionsReal("-tao_nls_gamma4_i", "", "", nlsP->gamma4_i, &nlsP->gamma4_i,NULL)); 783 CHKERRQ(PetscOptionsReal("-tao_nls_theta_i", "", "", nlsP->theta_i, &nlsP->theta_i,NULL)); 784 CHKERRQ(PetscOptionsReal("-tao_nls_mu1", "", "", nlsP->mu1, &nlsP->mu1,NULL)); 785 CHKERRQ(PetscOptionsReal("-tao_nls_mu2", "", "", nlsP->mu2, &nlsP->mu2,NULL)); 786 CHKERRQ(PetscOptionsReal("-tao_nls_gamma1", "", "", nlsP->gamma1, &nlsP->gamma1,NULL)); 787 CHKERRQ(PetscOptionsReal("-tao_nls_gamma2", "", "", nlsP->gamma2, &nlsP->gamma2,NULL)); 788 CHKERRQ(PetscOptionsReal("-tao_nls_gamma3", "", "", nlsP->gamma3, &nlsP->gamma3,NULL)); 789 CHKERRQ(PetscOptionsReal("-tao_nls_gamma4", "", "", nlsP->gamma4, &nlsP->gamma4,NULL)); 790 CHKERRQ(PetscOptionsReal("-tao_nls_theta", "", "", nlsP->theta, &nlsP->theta,NULL)); 791 CHKERRQ(PetscOptionsReal("-tao_nls_min_radius", "lower bound on initial radius", "", nlsP->min_radius, &nlsP->min_radius,NULL)); 792 CHKERRQ(PetscOptionsReal("-tao_nls_max_radius", "upper bound on radius", "", nlsP->max_radius, &nlsP->max_radius,NULL)); 793 CHKERRQ(PetscOptionsReal("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", nlsP->epsilon, &nlsP->epsilon,NULL)); 794 CHKERRQ(PetscOptionsTail()); 795 CHKERRQ(TaoLineSearchSetFromOptions(tao->linesearch)); 796 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 808 if (isascii) { 809 CHKERRQ(PetscViewerASCIIPushTab(viewer)); 810 CHKERRQ(PetscViewerASCIIPrintf(viewer, "Newton steps: %D\n", nlsP->newt)); 811 CHKERRQ(PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", nlsP->bfgs)); 812 CHKERRQ(PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", nlsP->grad)); 813 814 CHKERRQ(PetscViewerASCIIPrintf(viewer, "nls ksp atol: %D\n", nlsP->ksp_atol)); 815 CHKERRQ(PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %D\n", nlsP->ksp_rtol)); 816 CHKERRQ(PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %D\n", nlsP->ksp_ctol)); 817 CHKERRQ(PetscViewerASCIIPrintf(viewer, "nls ksp negc: %D\n", nlsP->ksp_negc)); 818 CHKERRQ(PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %D\n", nlsP->ksp_dtol)); 819 CHKERRQ(PetscViewerASCIIPrintf(viewer, "nls ksp iter: %D\n", nlsP->ksp_iter)); 820 CHKERRQ(PetscViewerASCIIPrintf(viewer, "nls ksp othr: %D\n", nlsP->ksp_othr)); 821 CHKERRQ(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 CHKERRQ(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 CHKERRQ(TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch)); 968 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch,(PetscObject)tao,1)); 969 CHKERRQ(TaoLineSearchSetType(tao->linesearch,morethuente_type)); 970 CHKERRQ(TaoLineSearchUseTaoRoutines(tao->linesearch,tao)); 971 CHKERRQ(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 972 973 /* Set linear solver to default for symmetric matrices */ 974 CHKERRQ(KSPCreate(((PetscObject)tao)->comm,&tao->ksp)); 975 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1)); 976 CHKERRQ(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 977 CHKERRQ(KSPAppendOptionsPrefix(tao->ksp,"tao_nls_")); 978 CHKERRQ(KSPSetType(tao->ksp,KSPSTCG)); 979 PetscFunctionReturn(0); 980 } 981