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