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