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