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