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