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