1 #include <../src/tao/matrix/lmvmmat.h> 2 #include <../src/tao/unconstrained/impls/ntl/ntl.h> 3 4 #include <petscksp.h> 5 #include <petscpc.h> 6 #include <petsc-private/kspimpl.h> 7 #include <petsc-private/pcimpl.h> 8 9 #define NTL_KSP_NASH 0 10 #define NTL_KSP_STCG 1 11 #define NTL_KSP_GLTR 2 12 #define NTL_KSP_TYPES 3 13 14 #define NTL_PC_NONE 0 15 #define NTL_PC_AHESS 1 16 #define NTL_PC_BFGS 2 17 #define NTL_PC_PETSC 3 18 #define NTL_PC_TYPES 4 19 20 #define BFGS_SCALE_AHESS 0 21 #define BFGS_SCALE_BFGS 1 22 #define BFGS_SCALE_TYPES 2 23 24 #define NTL_INIT_CONSTANT 0 25 #define NTL_INIT_DIRECTION 1 26 #define NTL_INIT_INTERPOLATION 2 27 #define NTL_INIT_TYPES 3 28 29 #define NTL_UPDATE_REDUCTION 0 30 #define NTL_UPDATE_INTERPOLATION 1 31 #define NTL_UPDATE_TYPES 2 32 33 static const char *NTL_KSP[64] = {"nash", "stcg", "gltr"}; 34 35 static const char *NTL_PC[64] = {"none", "ahess", "bfgs", "petsc"}; 36 37 static const char *BFGS_SCALE[64] = {"ahess", "bfgs"}; 38 39 static const char *NTL_INIT[64] = {"constant", "direction", "interpolation"}; 40 41 static const char *NTL_UPDATE[64] = {"reduction", "interpolation"}; 42 43 /* Routine for BFGS preconditioner */ 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "MatLMVMSolveShell" 47 static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x) 48 { 49 PetscErrorCode ierr; 50 Mat M; 51 52 PetscFunctionBegin; 53 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 54 PetscValidHeaderSpecific(b,VEC_CLASSID,2); 55 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 56 ierr = PCShellGetContext(pc,(void**)&M);CHKERRQ(ierr); 57 ierr = MatLMVMSolve(M, b, x);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 /* Implements Newton's Method with a trust-region, line-search approach for 62 solving unconstrained minimization problems. A More'-Thuente line search 63 is used to guarantee that the bfgs preconditioner remains positive 64 definite. */ 65 66 #define NTL_NEWTON 0 67 #define NTL_BFGS 1 68 #define NTL_SCALED_GRADIENT 2 69 #define NTL_GRADIENT 3 70 71 #undef __FUNCT__ 72 #define __FUNCT__ "TaoSolve_NTL" 73 static PetscErrorCode TaoSolve_NTL(Tao tao) 74 { 75 TAO_NTL *tl = (TAO_NTL *)tao->data; 76 PC pc; 77 KSPConvergedReason ksp_reason; 78 TaoConvergedReason reason; 79 TaoLineSearchConvergedReason ls_reason; 80 81 PetscReal fmin, ftrial, prered, actred, kappa, sigma; 82 PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius; 83 PetscReal f, fold, gdx, gnorm; 84 PetscReal step = 1.0; 85 86 PetscReal delta; 87 PetscReal norm_d = 0.0; 88 MatStructure matflag; 89 PetscErrorCode ierr; 90 PetscInt stepType; 91 PetscInt iter = 0,its; 92 93 PetscInt bfgsUpdates = 0; 94 PetscInt needH; 95 96 PetscInt i_max = 5; 97 PetscInt j_max = 1; 98 PetscInt i, j, n, N; 99 100 PetscInt tr_reject; 101 102 PetscFunctionBegin; 103 if (tao->XL || tao->XU || tao->ops->computebounds) { 104 ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n");CHKERRQ(ierr); 105 } 106 107 /* Initialize trust-region radius */ 108 tao->trust = tao->trust0; 109 110 /* Modify the radius if it is too large or small */ 111 tao->trust = PetscMax(tao->trust, tl->min_radius); 112 tao->trust = PetscMin(tao->trust, tl->max_radius); 113 114 if (NTL_PC_BFGS == tl->pc_type && !tl->M) { 115 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 116 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 117 ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tl->M);CHKERRQ(ierr); 118 ierr = MatLMVMAllocateVectors(tl->M,tao->solution);CHKERRQ(ierr); 119 } 120 121 /* Check convergence criteria */ 122 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr); 123 ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr); 124 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 125 needH = 1; 126 127 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr); 128 if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 129 130 /* Create vectors for the limited memory preconditioner */ 131 if ((NTL_PC_BFGS == tl->pc_type) && (BFGS_SCALE_BFGS != tl->bfgs_scale_type)) { 132 if (!tl->Diag) { 133 ierr = VecDuplicate(tao->solution, &tl->Diag);CHKERRQ(ierr); 134 } 135 } 136 137 /* Modify the linear solver to a conjugate gradient method */ 138 switch(tl->ksp_type) { 139 case NTL_KSP_NASH: 140 ierr = KSPSetType(tao->ksp, KSPNASH);CHKERRQ(ierr); 141 if (tao->ksp->ops->setfromoptions) { 142 (*tao->ksp->ops->setfromoptions)(tao->ksp); 143 } 144 break; 145 146 case NTL_KSP_STCG: 147 ierr = KSPSetType(tao->ksp, KSPSTCG);CHKERRQ(ierr); 148 if (tao->ksp->ops->setfromoptions) { 149 (*tao->ksp->ops->setfromoptions)(tao->ksp); 150 } 151 break; 152 153 default: 154 ierr = KSPSetType(tao->ksp, KSPGLTR);CHKERRQ(ierr); 155 if (tao->ksp->ops->setfromoptions) { 156 (*tao->ksp->ops->setfromoptions)(tao->ksp); 157 } 158 break; 159 } 160 161 /* Modify the preconditioner to use the bfgs approximation */ 162 ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr); 163 switch(tl->pc_type) { 164 case NTL_PC_NONE: 165 ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr); 166 if (pc->ops->setfromoptions) { 167 (*pc->ops->setfromoptions)(pc); 168 } 169 break; 170 171 case NTL_PC_AHESS: 172 ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr); 173 if (pc->ops->setfromoptions) { 174 (*pc->ops->setfromoptions)(pc); 175 } 176 ierr = PCJacobiSetUseAbs(pc);CHKERRQ(ierr); 177 break; 178 179 case NTL_PC_BFGS: 180 ierr = PCSetType(pc, PCSHELL);CHKERRQ(ierr); 181 if (pc->ops->setfromoptions) { 182 (*pc->ops->setfromoptions)(pc); 183 } 184 ierr = PCShellSetName(pc, "bfgs");CHKERRQ(ierr); 185 ierr = PCShellSetContext(pc, tl->M);CHKERRQ(ierr); 186 ierr = PCShellSetApply(pc, MatLMVMSolveShell);CHKERRQ(ierr); 187 break; 188 189 default: 190 /* Use the pc method set by pc_type */ 191 break; 192 } 193 194 /* Initialize trust-region radius. The initialization is only performed 195 when we are using Steihaug-Toint or the Generalized Lanczos method. */ 196 switch(tl->init_type) { 197 case NTL_INIT_CONSTANT: 198 /* Use the initial radius specified */ 199 break; 200 201 case NTL_INIT_INTERPOLATION: 202 /* Use the initial radius specified */ 203 max_radius = 0.0; 204 205 for (j = 0; j < j_max; ++j) { 206 fmin = f; 207 sigma = 0.0; 208 209 if (needH) { 210 ierr = TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag);CHKERRQ(ierr); 211 needH = 0; 212 } 213 214 for (i = 0; i < i_max; ++i) { 215 ierr = VecCopy(tao->solution, tl->W);CHKERRQ(ierr); 216 ierr = VecAXPY(tl->W, -tao->trust/gnorm, tao->gradient);CHKERRQ(ierr); 217 218 ierr = TaoComputeObjective(tao, tl->W, &ftrial);CHKERRQ(ierr); 219 if (PetscIsInfOrNanReal(ftrial)) { 220 tau = tl->gamma1_i; 221 } else { 222 if (ftrial < fmin) { 223 fmin = ftrial; 224 sigma = -tao->trust / gnorm; 225 } 226 227 ierr = MatMult(tao->hessian, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 228 ierr = VecDot(tao->gradient, tao->stepdirection, &prered);CHKERRQ(ierr); 229 230 prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm)); 231 actred = f - ftrial; 232 if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) { 233 kappa = 1.0; 234 } else { 235 kappa = actred / prered; 236 } 237 238 tau_1 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred); 239 tau_2 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred); 240 tau_min = PetscMin(tau_1, tau_2); 241 tau_max = PetscMax(tau_1, tau_2); 242 243 if (PetscAbsScalar(kappa - 1.0) <= tl->mu1_i) { 244 /* Great agreement */ 245 max_radius = PetscMax(max_radius, tao->trust); 246 247 if (tau_max < 1.0) { 248 tau = tl->gamma3_i; 249 } else if (tau_max > tl->gamma4_i) { 250 tau = tl->gamma4_i; 251 } else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) { 252 tau = tau_1; 253 } else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) { 254 tau = tau_2; 255 } else { 256 tau = tau_max; 257 } 258 } else if (PetscAbsScalar(kappa - 1.0) <= tl->mu2_i) { 259 /* Good agreement */ 260 max_radius = PetscMax(max_radius, tao->trust); 261 262 if (tau_max < tl->gamma2_i) { 263 tau = tl->gamma2_i; 264 } else if (tau_max > tl->gamma3_i) { 265 tau = tl->gamma3_i; 266 } else { 267 tau = tau_max; 268 } 269 } else { 270 /* Not good agreement */ 271 if (tau_min > 1.0) { 272 tau = tl->gamma2_i; 273 } else if (tau_max < tl->gamma1_i) { 274 tau = tl->gamma1_i; 275 } else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) { 276 tau = tl->gamma1_i; 277 } else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) { 278 tau = tau_1; 279 } else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) { 280 tau = tau_2; 281 } else { 282 tau = tau_max; 283 } 284 } 285 } 286 tao->trust = tau * tao->trust; 287 } 288 289 if (fmin < f) { 290 f = fmin; 291 ierr = VecAXPY(tao->solution, sigma, tao->gradient);CHKERRQ(ierr); 292 ierr = TaoComputeGradient(tao, tao->solution, tao->gradient);CHKERRQ(ierr); 293 294 ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr); 295 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 296 needH = 1; 297 298 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr); 299 if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 300 } 301 } 302 tao->trust = PetscMax(tao->trust, max_radius); 303 304 /* Modify the radius if it is too large or small */ 305 tao->trust = PetscMax(tao->trust, tl->min_radius); 306 tao->trust = PetscMin(tao->trust, tl->max_radius); 307 break; 308 309 default: 310 /* Norm of the first direction will initialize radius */ 311 tao->trust = 0.0; 312 break; 313 } 314 315 /* Set initial scaling for the BFGS preconditioner 316 This step is done after computing the initial trust-region radius 317 since the function value may have decreased */ 318 if (NTL_PC_BFGS == tl->pc_type) { 319 if (f != 0.0) { 320 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 321 } else { 322 delta = 2.0 / (gnorm*gnorm); 323 } 324 ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr); 325 } 326 327 /* Set counter for gradient/reset steps */ 328 tl->ntrust = 0; 329 tl->newt = 0; 330 tl->bfgs = 0; 331 tl->sgrad = 0; 332 tl->grad = 0; 333 334 /* Have not converged; continue with Newton method */ 335 while (reason == TAO_CONTINUE_ITERATING) { 336 ++iter; 337 338 /* Compute the Hessian */ 339 if (needH) { 340 ierr = TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag);CHKERRQ(ierr); 341 needH = 0; 342 } 343 344 if (NTL_PC_BFGS == tl->pc_type) { 345 if (BFGS_SCALE_AHESS == tl->bfgs_scale_type) { 346 /* Obtain diagonal for the bfgs preconditioner */ 347 ierr = MatGetDiagonal(tao->hessian, tl->Diag);CHKERRQ(ierr); 348 ierr = VecAbs(tl->Diag);CHKERRQ(ierr); 349 ierr = VecReciprocal(tl->Diag);CHKERRQ(ierr); 350 ierr = MatLMVMSetScale(tl->M, tl->Diag);CHKERRQ(ierr); 351 } 352 353 /* Update the limited memory preconditioner */ 354 ierr = MatLMVMUpdate(tl->M,tao->solution, tao->gradient);CHKERRQ(ierr); 355 ++bfgsUpdates; 356 } 357 ierr = KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre, matflag);CHKERRQ(ierr); 358 /* Solve the Newton system of equations */ 359 if (NTL_KSP_NASH == tl->ksp_type) { 360 ierr = KSPNASHSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr); 361 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 362 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 363 tao->ksp_its+=its; 364 ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 365 } else if (NTL_KSP_STCG == tl->ksp_type) { 366 ierr = KSPSTCGSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr); 367 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 368 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 369 tao->ksp_its+=its; 370 ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 371 } else { /* NTL_KSP_GLTR */ 372 ierr = KSPGLTRSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr); 373 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 374 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 375 tao->ksp_its+=its; 376 ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 377 } 378 379 if (0.0 == tao->trust) { 380 /* Radius was uninitialized; use the norm of the direction */ 381 if (norm_d > 0.0) { 382 tao->trust = norm_d; 383 384 /* Modify the radius if it is too large or small */ 385 tao->trust = PetscMax(tao->trust, tl->min_radius); 386 tao->trust = PetscMin(tao->trust, tl->max_radius); 387 } else { 388 /* The direction was bad; set radius to default value and re-solve 389 the trust-region subproblem to get a direction */ 390 tao->trust = tao->trust0; 391 392 /* Modify the radius if it is too large or small */ 393 tao->trust = PetscMax(tao->trust, tl->min_radius); 394 tao->trust = PetscMin(tao->trust, tl->max_radius); 395 396 if (NTL_KSP_NASH == tl->ksp_type) { 397 ierr = KSPNASHSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr); 398 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 399 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 400 tao->ksp_its+=its; 401 ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 402 } else if (NTL_KSP_STCG == tl->ksp_type) { 403 ierr = KSPSTCGSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr); 404 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 405 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 406 tao->ksp_its+=its; 407 ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 408 } else { /* NTL_KSP_GLTR */ 409 ierr = KSPGLTRSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr); 410 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 411 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 412 tao->ksp_its+=its; 413 ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 414 } 415 416 417 if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero"); 418 } 419 } 420 421 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 422 ierr = KSPGetConvergedReason(tao->ksp, &ksp_reason);CHKERRQ(ierr); 423 if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (NTL_PC_BFGS == tl->pc_type) && (bfgsUpdates > 1)) { 424 /* Preconditioner is numerically indefinite; reset the 425 approximate if using BFGS preconditioning. */ 426 427 if (f != 0.0) { 428 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 429 } else { 430 delta = 2.0 / (gnorm*gnorm); 431 } 432 ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr); 433 ierr = MatLMVMReset(tl->M);CHKERRQ(ierr); 434 ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr); 435 bfgsUpdates = 1; 436 } 437 438 /* Check trust-region reduction conditions */ 439 tr_reject = 0; 440 if (NTL_UPDATE_REDUCTION == tl->update_type) { 441 /* Get predicted reduction */ 442 if (NTL_KSP_NASH == tl->ksp_type) { 443 ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 444 } else if (NTL_KSP_STCG == tl->ksp_type) { 445 ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 446 } else { /* gltr */ 447 ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 448 } 449 450 if (prered >= 0.0) { 451 /* The predicted reduction has the wrong sign. This cannot 452 happen in infinite precision arithmetic. Step should 453 be rejected! */ 454 tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d); 455 tr_reject = 1; 456 } else { 457 /* Compute trial step and function value */ 458 ierr = VecCopy(tao->solution, tl->W);CHKERRQ(ierr); 459 ierr = VecAXPY(tl->W, 1.0, tao->stepdirection);CHKERRQ(ierr); 460 ierr = TaoComputeObjective(tao, tl->W, &ftrial);CHKERRQ(ierr); 461 462 if (PetscIsInfOrNanReal(ftrial)) { 463 tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d); 464 tr_reject = 1; 465 } else { 466 /* Compute and actual reduction */ 467 actred = f - ftrial; 468 prered = -prered; 469 if ((PetscAbsScalar(actred) <= tl->epsilon) && 470 (PetscAbsScalar(prered) <= tl->epsilon)) { 471 kappa = 1.0; 472 } else { 473 kappa = actred / prered; 474 } 475 476 /* Accept of reject the step and update radius */ 477 if (kappa < tl->eta1) { 478 /* Reject the step */ 479 tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d); 480 tr_reject = 1; 481 } else { 482 /* Accept the step */ 483 if (kappa < tl->eta2) { 484 /* Marginal bad step */ 485 tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d); 486 } else if (kappa < tl->eta3) { 487 /* Reasonable step */ 488 tao->trust = tl->alpha3 * tao->trust; 489 } else if (kappa < tl->eta4) { 490 /* Good step */ 491 tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust); 492 } else { 493 /* Very good step */ 494 tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust); 495 } 496 } 497 } 498 } 499 } else { 500 /* Get predicted reduction */ 501 if (NTL_KSP_NASH == tl->ksp_type) { 502 ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 503 } else if (NTL_KSP_STCG == tl->ksp_type) { 504 ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 505 } else { /* gltr */ 506 ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 507 } 508 509 if (prered >= 0.0) { 510 /* The predicted reduction has the wrong sign. This cannot 511 happen in infinite precision arithmetic. Step should 512 be rejected! */ 513 tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d); 514 tr_reject = 1; 515 } else { 516 ierr = VecCopy(tao->solution, tl->W);CHKERRQ(ierr); 517 ierr = VecAXPY(tl->W, 1.0, tao->stepdirection);CHKERRQ(ierr); 518 ierr = TaoComputeObjective(tao, tl->W, &ftrial);CHKERRQ(ierr); 519 if (PetscIsInfOrNanReal(ftrial)) { 520 tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d); 521 tr_reject = 1; 522 } else { 523 ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr); 524 525 actred = f - ftrial; 526 prered = -prered; 527 if ((PetscAbsScalar(actred) <= tl->epsilon) && 528 (PetscAbsScalar(prered) <= tl->epsilon)) { 529 kappa = 1.0; 530 } else { 531 kappa = actred / prered; 532 } 533 534 tau_1 = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred); 535 tau_2 = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred); 536 tau_min = PetscMin(tau_1, tau_2); 537 tau_max = PetscMax(tau_1, tau_2); 538 539 if (kappa >= 1.0 - tl->mu1) { 540 /* Great agreement; accept step and update radius */ 541 if (tau_max < 1.0) { 542 tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d); 543 } else if (tau_max > tl->gamma4) { 544 tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d); 545 } else { 546 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 547 } 548 } else if (kappa >= 1.0 - tl->mu2) { 549 /* Good agreement */ 550 551 if (tau_max < tl->gamma2) { 552 tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d); 553 } else if (tau_max > tl->gamma3) { 554 tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d); 555 } else if (tau_max < 1.0) { 556 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 557 } else { 558 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 559 } 560 } else { 561 /* Not good agreement */ 562 if (tau_min > 1.0) { 563 tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d); 564 } else if (tau_max < tl->gamma1) { 565 tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d); 566 } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) { 567 tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d); 568 } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) { 569 tao->trust = tau_1 * PetscMin(tao->trust, norm_d); 570 } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) { 571 tao->trust = tau_2 * PetscMin(tao->trust, norm_d); 572 } else { 573 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 574 } 575 tr_reject = 1; 576 } 577 } 578 } 579 } 580 581 if (tr_reject) { 582 /* The trust-region constraints rejected the step. Apply a linesearch. 583 Check for descent direction. */ 584 ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr); 585 if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) { 586 /* Newton step is not descent or direction produced Inf or NaN */ 587 588 if (NTL_PC_BFGS != tl->pc_type) { 589 /* We don't have the bfgs matrix around and updated 590 Must use gradient direction in this case */ 591 ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr); 592 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 593 ++tl->grad; 594 stepType = NTL_GRADIENT; 595 } else { 596 /* Attempt to use the BFGS direction */ 597 ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 598 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 599 600 /* Check for success (descent direction) */ 601 ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr); 602 if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) { 603 /* BFGS direction is not descent or direction produced not a number 604 We can assert bfgsUpdates > 1 in this case because 605 the first solve produces the scaled gradient direction, 606 which is guaranteed to be descent */ 607 608 /* Use steepest descent direction (scaled) */ 609 if (f != 0.0) { 610 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 611 } else { 612 delta = 2.0 / (gnorm*gnorm); 613 } 614 ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr); 615 ierr = MatLMVMReset(tl->M);CHKERRQ(ierr); 616 ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr); 617 ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 618 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 619 620 bfgsUpdates = 1; 621 ++tl->sgrad; 622 stepType = NTL_SCALED_GRADIENT; 623 } else { 624 if (1 == bfgsUpdates) { 625 /* The first BFGS direction is always the scaled gradient */ 626 ++tl->sgrad; 627 stepType = NTL_SCALED_GRADIENT; 628 } else { 629 ++tl->bfgs; 630 stepType = NTL_BFGS; 631 } 632 } 633 } 634 } else { 635 /* Computed Newton step is descent */ 636 ++tl->newt; 637 stepType = NTL_NEWTON; 638 } 639 640 /* Perform the linesearch */ 641 fold = f; 642 ierr = VecCopy(tao->solution, tl->Xold);CHKERRQ(ierr); 643 ierr = VecCopy(tao->gradient, tl->Gold);CHKERRQ(ierr); 644 645 step = 1.0; 646 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);CHKERRQ(ierr); 647 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 648 649 while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) { /* Linesearch failed */ 650 /* Linesearch failed */ 651 f = fold; 652 ierr = VecCopy(tl->Xold, tao->solution);CHKERRQ(ierr); 653 ierr = VecCopy(tl->Gold, tao->gradient);CHKERRQ(ierr); 654 655 switch(stepType) { 656 case NTL_NEWTON: 657 /* Failed to obtain acceptable iterate with Newton step */ 658 659 if (NTL_PC_BFGS != tl->pc_type) { 660 /* We don't have the bfgs matrix around and being updated 661 Must use gradient direction in this case */ 662 ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr); 663 ++tl->grad; 664 stepType = NTL_GRADIENT; 665 } else { 666 /* Attempt to use the BFGS direction */ 667 ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 668 669 670 /* Check for success (descent direction) */ 671 ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr); 672 if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) { 673 /* BFGS direction is not descent or direction produced 674 not a number. We can assert bfgsUpdates > 1 in this case 675 Use steepest descent direction (scaled) */ 676 677 if (f != 0.0) { 678 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 679 } else { 680 delta = 2.0 / (gnorm*gnorm); 681 } 682 ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr); 683 ierr = MatLMVMReset(tl->M);CHKERRQ(ierr); 684 ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr); 685 ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 686 687 bfgsUpdates = 1; 688 ++tl->sgrad; 689 stepType = NTL_SCALED_GRADIENT; 690 } else { 691 if (1 == bfgsUpdates) { 692 /* The first BFGS direction is always the scaled gradient */ 693 ++tl->sgrad; 694 stepType = NTL_SCALED_GRADIENT; 695 } else { 696 ++tl->bfgs; 697 stepType = NTL_BFGS; 698 } 699 } 700 } 701 break; 702 703 case NTL_BFGS: 704 /* Can only enter if pc_type == NTL_PC_BFGS 705 Failed to obtain acceptable iterate with BFGS step 706 Attempt to use the scaled gradient direction */ 707 708 if (f != 0.0) { 709 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 710 } else { 711 delta = 2.0 / (gnorm*gnorm); 712 } 713 ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr); 714 ierr = MatLMVMReset(tl->M);CHKERRQ(ierr); 715 ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr); 716 ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 717 718 bfgsUpdates = 1; 719 ++tl->sgrad; 720 stepType = NTL_SCALED_GRADIENT; 721 break; 722 723 case NTL_SCALED_GRADIENT: 724 /* Can only enter if pc_type == NTL_PC_BFGS 725 The scaled gradient step did not produce a new iterate; 726 attemp to use the gradient direction. 727 Need to make sure we are not using a different diagonal scaling */ 728 ierr = MatLMVMSetScale(tl->M, tl->Diag);CHKERRQ(ierr); 729 ierr = MatLMVMSetDelta(tl->M, 1.0);CHKERRQ(ierr); 730 ierr = MatLMVMReset(tl->M);CHKERRQ(ierr); 731 ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr); 732 ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 733 734 bfgsUpdates = 1; 735 ++tl->grad; 736 stepType = NTL_GRADIENT; 737 break; 738 } 739 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 740 741 /* This may be incorrect; linesearch has values for stepmax and stepmin 742 that should be reset. */ 743 step = 1.0; 744 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);CHKERRQ(ierr); 745 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 746 } 747 748 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 749 /* Failed to find an improving point */ 750 f = fold; 751 ierr = VecCopy(tl->Xold, tao->solution);CHKERRQ(ierr); 752 ierr = VecCopy(tl->Gold, tao->gradient);CHKERRQ(ierr); 753 tao->trust = 0.0; 754 step = 0.0; 755 reason = TAO_DIVERGED_LS_FAILURE; 756 tao->reason = TAO_DIVERGED_LS_FAILURE; 757 break; 758 } else if (stepType == NTL_NEWTON) { 759 if (step < tl->nu1) { 760 /* Very bad step taken; reduce radius */ 761 tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust); 762 } else if (step < tl->nu2) { 763 /* Reasonably bad step taken; reduce radius */ 764 tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust); 765 } else if (step < tl->nu3) { 766 /* Reasonable step was taken; leave radius alone */ 767 if (tl->omega3 < 1.0) { 768 tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust); 769 } else if (tl->omega3 > 1.0) { 770 tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust); 771 } 772 } else if (step < tl->nu4) { 773 /* Full step taken; increase the radius */ 774 tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust); 775 } else { 776 /* More than full step taken; increase the radius */ 777 tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust); 778 } 779 } else { 780 /* Newton step was not good; reduce the radius */ 781 tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust); 782 } 783 } else { 784 /* Trust-region step is accepted */ 785 ierr = VecCopy(tl->W, tao->solution);CHKERRQ(ierr); 786 f = ftrial; 787 ierr = TaoComputeGradient(tao, tao->solution, tao->gradient);CHKERRQ(ierr); 788 ++tl->ntrust; 789 } 790 791 /* The radius may have been increased; modify if it is too large */ 792 tao->trust = PetscMin(tao->trust, tl->max_radius); 793 794 /* Check for converged */ 795 ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr); 796 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 797 needH = 1; 798 799 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, tao->trust, &reason);CHKERRQ(ierr); 800 } 801 PetscFunctionReturn(0); 802 } 803 804 /* ---------------------------------------------------------- */ 805 #undef __FUNCT__ 806 #define __FUNCT__ "TaoSetUp_NTL" 807 static PetscErrorCode TaoSetUp_NTL(Tao tao) 808 { 809 TAO_NTL *tl = (TAO_NTL *)tao->data; 810 PetscErrorCode ierr; 811 812 PetscFunctionBegin; 813 if (!tao->gradient) {ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr); } 814 if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);} 815 if (!tl->W) { ierr = VecDuplicate(tao->solution, &tl->W);CHKERRQ(ierr);} 816 if (!tl->Xold) { ierr = VecDuplicate(tao->solution, &tl->Xold);CHKERRQ(ierr);} 817 if (!tl->Gold) { ierr = VecDuplicate(tao->solution, &tl->Gold);CHKERRQ(ierr);} 818 tl->Diag = 0; 819 tl->M = 0; 820 PetscFunctionReturn(0); 821 } 822 823 /*------------------------------------------------------------*/ 824 #undef __FUNCT__ 825 #define __FUNCT__ "TaoDestroy_NTL" 826 static PetscErrorCode TaoDestroy_NTL(Tao tao) 827 { 828 TAO_NTL *tl = (TAO_NTL *)tao->data; 829 PetscErrorCode ierr; 830 831 PetscFunctionBegin; 832 if (tao->setupcalled) { 833 ierr = VecDestroy(&tl->W);CHKERRQ(ierr); 834 ierr = VecDestroy(&tl->Xold);CHKERRQ(ierr); 835 ierr = VecDestroy(&tl->Gold);CHKERRQ(ierr); 836 } 837 ierr = VecDestroy(&tl->Diag);CHKERRQ(ierr); 838 ierr = MatDestroy(&tl->M);CHKERRQ(ierr); 839 ierr = PetscFree(tao->data);CHKERRQ(ierr); 840 PetscFunctionReturn(0); 841 } 842 843 /*------------------------------------------------------------*/ 844 #undef __FUNCT__ 845 #define __FUNCT__ "TaoSetFromOptions_NTL" 846 static PetscErrorCode TaoSetFromOptions_NTL(Tao tao) 847 { 848 TAO_NTL *tl = (TAO_NTL *)tao->data; 849 PetscErrorCode ierr; 850 851 PetscFunctionBegin; 852 ierr = PetscOptionsHead("Newton line search method for unconstrained optimization");CHKERRQ(ierr); 853 ierr = PetscOptionsEList("-tao_ntl_ksp_type", "ksp type", "", NTL_KSP, NTL_KSP_TYPES, NTL_KSP[tl->ksp_type], &tl->ksp_type, 0);CHKERRQ(ierr); 854 ierr = PetscOptionsEList("-tao_ntl_pc_type", "pc type", "", NTL_PC, NTL_PC_TYPES, NTL_PC[tl->pc_type], &tl->pc_type, 0);CHKERRQ(ierr); 855 ierr = PetscOptionsEList("-tao_ntl_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tl->bfgs_scale_type], &tl->bfgs_scale_type, 0);CHKERRQ(ierr); 856 ierr = PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type, 0);CHKERRQ(ierr); 857 ierr = PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type, 0);CHKERRQ(ierr); 858 ierr = PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1, 0);CHKERRQ(ierr); 859 ierr = PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2, 0);CHKERRQ(ierr); 860 ierr = PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3, 0);CHKERRQ(ierr); 861 ierr = PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4, 0);CHKERRQ(ierr); 862 ierr = PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1, 0);CHKERRQ(ierr); 863 ierr = PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2, 0);CHKERRQ(ierr); 864 ierr = PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3, 0);CHKERRQ(ierr); 865 ierr = PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4, 0);CHKERRQ(ierr); 866 ierr = PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5, 0);CHKERRQ(ierr); 867 ierr = PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1, 0);CHKERRQ(ierr); 868 ierr = PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2, 0);CHKERRQ(ierr); 869 ierr = PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3, 0);CHKERRQ(ierr); 870 ierr = PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4, 0);CHKERRQ(ierr); 871 ierr = PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1, 0);CHKERRQ(ierr); 872 ierr = PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2, 0);CHKERRQ(ierr); 873 ierr = PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3, 0);CHKERRQ(ierr); 874 ierr = PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4, 0);CHKERRQ(ierr); 875 ierr = PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5, 0);CHKERRQ(ierr); 876 ierr = PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i, 0);CHKERRQ(ierr); 877 ierr = PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i, 0);CHKERRQ(ierr); 878 ierr = PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i, 0);CHKERRQ(ierr); 879 ierr = PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i, 0);CHKERRQ(ierr); 880 ierr = PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i, 0);CHKERRQ(ierr); 881 ierr = PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i, 0);CHKERRQ(ierr); 882 ierr = PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i, 0);CHKERRQ(ierr); 883 ierr = PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1, 0);CHKERRQ(ierr); 884 ierr = PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2, 0);CHKERRQ(ierr); 885 ierr = PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1, 0);CHKERRQ(ierr); 886 ierr = PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2, 0);CHKERRQ(ierr); 887 ierr = PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3, 0);CHKERRQ(ierr); 888 ierr = PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4, 0);CHKERRQ(ierr); 889 ierr = PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta, 0);CHKERRQ(ierr); 890 ierr = PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius, 0);CHKERRQ(ierr); 891 ierr = PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius, 0);CHKERRQ(ierr); 892 ierr = PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon, 0);CHKERRQ(ierr); 893 ierr = PetscOptionsTail();CHKERRQ(ierr); 894 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 895 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 896 PetscFunctionReturn(0); 897 } 898 899 /*------------------------------------------------------------*/ 900 #undef __FUNCT__ 901 #define __FUNCT__ "TaoView_NTL" 902 static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer) 903 { 904 TAO_NTL *tl = (TAO_NTL *)tao->data; 905 PetscInt nrejects; 906 PetscBool isascii; 907 PetscErrorCode ierr; 908 909 PetscFunctionBegin; 910 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 911 if (isascii) { 912 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 913 if (NTL_PC_BFGS == tl->pc_type && tl->M) { 914 ierr = MatLMVMGetRejects(tl->M, &nrejects);CHKERRQ(ierr); 915 ierr = PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n", nrejects);CHKERRQ(ierr); 916 } 917 ierr = PetscViewerASCIIPrintf(viewer, "Trust-region steps: %D\n", tl->ntrust);CHKERRQ(ierr); 918 ierr = PetscViewerASCIIPrintf(viewer, "Newton search steps: %D\n", tl->newt);CHKERRQ(ierr); 919 ierr = PetscViewerASCIIPrintf(viewer, "BFGS search steps: %D\n", tl->bfgs);CHKERRQ(ierr); 920 ierr = PetscViewerASCIIPrintf(viewer, "Scaled gradient search steps: %D\n", tl->sgrad);CHKERRQ(ierr); 921 ierr = PetscViewerASCIIPrintf(viewer, "Gradient search steps: %D\n", tl->grad);CHKERRQ(ierr); 922 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 923 } 924 PetscFunctionReturn(0); 925 } 926 927 /* ---------------------------------------------------------- */ 928 EXTERN_C_BEGIN 929 #undef __FUNCT__ 930 #define __FUNCT__ "TaoCreate_NTL" 931 PetscErrorCode TaoCreate_NTL(Tao tao) 932 { 933 TAO_NTL *tl; 934 PetscErrorCode ierr; 935 const char *morethuente_type = TAOLINESEARCH_MT; 936 937 PetscFunctionBegin; 938 ierr = PetscNewLog(tao,&tl);CHKERRQ(ierr); 939 tao->ops->setup = TaoSetUp_NTL; 940 tao->ops->solve = TaoSolve_NTL; 941 tao->ops->view = TaoView_NTL; 942 tao->ops->setfromoptions = TaoSetFromOptions_NTL; 943 tao->ops->destroy = TaoDestroy_NTL; 944 945 tao->max_it = 50; 946 #if defined(PETSC_USE_REAL_SINGLE) 947 tao->fatol = 1e-5; 948 tao->frtol = 1e-5; 949 #else 950 tao->fatol = 1e-10; 951 tao->frtol = 1e-10; 952 #endif 953 tao->data = (void*)tl; 954 955 tao->trust0 = 100.0; 956 957 958 /* Default values for trust-region radius update based on steplength */ 959 tl->nu1 = 0.25; 960 tl->nu2 = 0.50; 961 tl->nu3 = 1.00; 962 tl->nu4 = 1.25; 963 964 tl->omega1 = 0.25; 965 tl->omega2 = 0.50; 966 tl->omega3 = 1.00; 967 tl->omega4 = 2.00; 968 tl->omega5 = 4.00; 969 970 /* Default values for trust-region radius update based on reduction */ 971 tl->eta1 = 1.0e-4; 972 tl->eta2 = 0.25; 973 tl->eta3 = 0.50; 974 tl->eta4 = 0.90; 975 976 tl->alpha1 = 0.25; 977 tl->alpha2 = 0.50; 978 tl->alpha3 = 1.00; 979 tl->alpha4 = 2.00; 980 tl->alpha5 = 4.00; 981 982 /* Default values for trust-region radius update based on interpolation */ 983 tl->mu1 = 0.10; 984 tl->mu2 = 0.50; 985 986 tl->gamma1 = 0.25; 987 tl->gamma2 = 0.50; 988 tl->gamma3 = 2.00; 989 tl->gamma4 = 4.00; 990 991 tl->theta = 0.05; 992 993 /* Default values for trust region initialization based on interpolation */ 994 tl->mu1_i = 0.35; 995 tl->mu2_i = 0.50; 996 997 tl->gamma1_i = 0.0625; 998 tl->gamma2_i = 0.5; 999 tl->gamma3_i = 2.0; 1000 tl->gamma4_i = 5.0; 1001 1002 tl->theta_i = 0.25; 1003 1004 /* Remaining parameters */ 1005 tl->min_radius = 1.0e-10; 1006 tl->max_radius = 1.0e10; 1007 tl->epsilon = 1.0e-6; 1008 1009 tl->ksp_type = NTL_KSP_STCG; 1010 tl->pc_type = NTL_PC_BFGS; 1011 tl->bfgs_scale_type = BFGS_SCALE_AHESS; 1012 tl->init_type = NTL_INIT_INTERPOLATION; 1013 tl->update_type = NTL_UPDATE_REDUCTION; 1014 1015 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr); 1016 ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type);CHKERRQ(ierr); 1017 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch, tao);CHKERRQ(ierr); 1018 ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr); 1019 PetscFunctionReturn(0); 1020 } 1021 EXTERN_C_END 1022 1023 1024 1025