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