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