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