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