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