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