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