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