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