1 #include <petsctaolinesearch.h> 2 #include <../src/tao/bound/impls/bnk/bnk.h> 3 #include <petscksp.h> 4 5 static const char *BNK_INIT[64] = {"constant", "direction", "interpolation"}; 6 static const char *BNK_UPDATE[64] = {"step", "reduction", "interpolation"}; 7 static const char *BNK_AS[64] = {"none", "bertsekas"}; 8 9 /*------------------------------------------------------------*/ 10 11 /* Routine for initializing the KSP solver, the BFGS preconditioner, and the initial trust radius estimation */ 12 13 PetscErrorCode TaoBNKInitialize(Tao tao, PetscInt initType, PetscBool *needH) 14 { 15 TAO_BNK *bnk = (TAO_BNK *)tao->data; 16 PC pc; 17 PetscReal f_min, ftrial, prered, actred, kappa, sigma, resnorm; 18 PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius; 19 PetscBool is_bfgs, is_jacobi, is_symmetric, sym_set; 20 PetscInt n, N, nDiff; 21 PetscInt i_max = 5; 22 PetscInt j_max = 1; 23 PetscInt i, j; 24 PetscVoidFunction kspTR; 25 26 PetscFunctionBegin; 27 /* Project the current point onto the feasible set */ 28 PetscCall(TaoComputeVariableBounds(tao)); 29 PetscCall(TaoSetVariableBounds(bnk->bncg, tao->XL, tao->XU)); 30 if (tao->bounded) PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU)); 31 32 /* Project the initial point onto the feasible region */ 33 PetscCall(TaoBoundSolution(tao->solution, tao->XL, tao->XU, 0.0, &nDiff, tao->solution)); 34 35 /* Check convergence criteria */ 36 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &bnk->f, bnk->unprojected_gradient)); 37 PetscCall(TaoBNKEstimateActiveSet(tao, bnk->as_type)); 38 PetscCall(VecCopy(bnk->unprojected_gradient, tao->gradient)); 39 PetscCall(VecISSet(tao->gradient, bnk->active_idx, 0.0)); 40 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm)); 41 42 /* Test the initial point for convergence */ 43 PetscCall(VecFischer(tao->solution, bnk->unprojected_gradient, tao->XL, tao->XU, bnk->W)); 44 PetscCall(VecNorm(bnk->W, NORM_2, &resnorm)); 45 PetscCheck(!PetscIsInfOrNanReal(bnk->f) && !PetscIsInfOrNanReal(resnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 46 PetscCall(TaoLogConvergenceHistory(tao, bnk->f, resnorm, 0.0, tao->ksp_its)); 47 PetscCall(TaoMonitor(tao, tao->niter, bnk->f, resnorm, 0.0, 1.0)); 48 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 49 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 50 51 /* Reset KSP stopping reason counters */ 52 bnk->ksp_atol = 0; 53 bnk->ksp_rtol = 0; 54 bnk->ksp_dtol = 0; 55 bnk->ksp_ctol = 0; 56 bnk->ksp_negc = 0; 57 bnk->ksp_iter = 0; 58 bnk->ksp_othr = 0; 59 60 /* Reset accepted step type counters */ 61 bnk->tot_cg_its = 0; 62 bnk->newt = 0; 63 bnk->bfgs = 0; 64 bnk->sgrad = 0; 65 bnk->grad = 0; 66 67 /* Initialize the Hessian perturbation */ 68 bnk->pert = bnk->sval; 69 70 /* Reset initial steplength to zero (this helps BNCG reset its direction internally) */ 71 PetscCall(VecSet(tao->stepdirection, 0.0)); 72 73 /* Allocate the vectors needed for the BFGS approximation */ 74 PetscCall(KSPGetPC(tao->ksp, &pc)); 75 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs)); 76 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi)); 77 if (is_bfgs) { 78 bnk->bfgs_pre = pc; 79 PetscCall(PCLMVMGetMatLMVM(bnk->bfgs_pre, &bnk->M)); 80 PetscCall(VecGetLocalSize(tao->solution, &n)); 81 PetscCall(VecGetSize(tao->solution, &N)); 82 PetscCall(MatSetSizes(bnk->M, n, n, N, N)); 83 PetscCall(MatLMVMAllocate(bnk->M, tao->solution, bnk->unprojected_gradient)); 84 PetscCall(MatIsSymmetricKnown(bnk->M, &sym_set, &is_symmetric)); 85 PetscCheck(sym_set && is_symmetric, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric."); 86 } else if (is_jacobi) PetscCall(PCJacobiSetUseAbs(pc, PETSC_TRUE)); 87 88 /* Prepare the min/max vectors for safeguarding diagonal scales */ 89 PetscCall(VecSet(bnk->Diag_min, bnk->dmin)); 90 PetscCall(VecSet(bnk->Diag_max, bnk->dmax)); 91 92 /* Initialize trust-region radius. The initialization is only performed 93 when we are using Nash, Steihaug-Toint or the Generalized Lanczos method. */ 94 *needH = PETSC_TRUE; 95 PetscCall(PetscObjectQueryFunction((PetscObject)tao->ksp, "KSPCGSetRadius_C", &kspTR)); 96 if (kspTR) { 97 switch (initType) { 98 case BNK_INIT_CONSTANT: 99 /* Use the initial radius specified */ 100 tao->trust = tao->trust0; 101 break; 102 103 case BNK_INIT_INTERPOLATION: 104 /* Use interpolation based on the initial Hessian */ 105 max_radius = 0.0; 106 tao->trust = tao->trust0; 107 for (j = 0; j < j_max; ++j) { 108 f_min = bnk->f; 109 sigma = 0.0; 110 111 if (*needH) { 112 /* Compute the Hessian at the new step, and extract the inactive subsystem */ 113 PetscCall((*bnk->computehessian)(tao)); 114 PetscCall(TaoBNKEstimateActiveSet(tao, BNK_AS_NONE)); 115 PetscCall(MatDestroy(&bnk->H_inactive)); 116 if (bnk->active_idx) { 117 PetscCall(MatCreateSubMatrix(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, MAT_INITIAL_MATRIX, &bnk->H_inactive)); 118 } else { 119 PetscCall(PetscObjectReference((PetscObject)tao->hessian)); 120 bnk->H_inactive = tao->hessian; 121 } 122 *needH = PETSC_FALSE; 123 } 124 125 for (i = 0; i < i_max; ++i) { 126 /* Take a steepest descent step and snap it to bounds */ 127 PetscCall(VecCopy(tao->solution, bnk->Xold)); 128 PetscCall(VecAXPY(tao->solution, -tao->trust / bnk->gnorm, tao->gradient)); 129 PetscCall(TaoBoundSolution(tao->solution, tao->XL, tao->XU, 0.0, &nDiff, tao->solution)); 130 /* Compute the step we actually accepted */ 131 PetscCall(VecCopy(tao->solution, bnk->W)); 132 PetscCall(VecAXPY(bnk->W, -1.0, bnk->Xold)); 133 /* Compute the objective at the trial */ 134 PetscCall(TaoComputeObjective(tao, tao->solution, &ftrial)); 135 PetscCheck(!PetscIsInfOrNanReal(bnk->f), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 136 PetscCall(VecCopy(bnk->Xold, tao->solution)); 137 if (PetscIsInfOrNanReal(ftrial)) { 138 tau = bnk->gamma1_i; 139 } else { 140 if (ftrial < f_min) { 141 f_min = ftrial; 142 sigma = -tao->trust / bnk->gnorm; 143 } 144 145 /* Compute the predicted and actual reduction */ 146 if (bnk->active_idx) { 147 PetscCall(VecGetSubVector(bnk->W, bnk->inactive_idx, &bnk->X_inactive)); 148 PetscCall(VecGetSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work)); 149 } else { 150 bnk->X_inactive = bnk->W; 151 bnk->inactive_work = bnk->Xwork; 152 } 153 PetscCall(MatMult(bnk->H_inactive, bnk->X_inactive, bnk->inactive_work)); 154 PetscCall(VecDot(bnk->X_inactive, bnk->inactive_work, &prered)); 155 if (bnk->active_idx) { 156 PetscCall(VecRestoreSubVector(bnk->W, bnk->inactive_idx, &bnk->X_inactive)); 157 PetscCall(VecRestoreSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work)); 158 } 159 prered = tao->trust * (bnk->gnorm - 0.5 * tao->trust * prered / (bnk->gnorm * bnk->gnorm)); 160 actred = bnk->f - ftrial; 161 if ((PetscAbsScalar(actred) <= bnk->epsilon) && (PetscAbsScalar(prered) <= bnk->epsilon)) { 162 kappa = 1.0; 163 } else { 164 kappa = actred / prered; 165 } 166 167 tau_1 = bnk->theta_i * bnk->gnorm * tao->trust / (bnk->theta_i * bnk->gnorm * tao->trust + (1.0 - bnk->theta_i) * prered - actred); 168 tau_2 = bnk->theta_i * bnk->gnorm * tao->trust / (bnk->theta_i * bnk->gnorm * tao->trust - (1.0 + bnk->theta_i) * prered + actred); 169 tau_min = PetscMin(tau_1, tau_2); 170 tau_max = PetscMax(tau_1, tau_2); 171 172 if (PetscAbsScalar(kappa - (PetscReal)1.0) <= bnk->mu1_i) { 173 /* Great agreement */ 174 max_radius = PetscMax(max_radius, tao->trust); 175 176 if (tau_max < 1.0) { 177 tau = bnk->gamma3_i; 178 } else if (tau_max > bnk->gamma4_i) { 179 tau = bnk->gamma4_i; 180 } else { 181 tau = tau_max; 182 } 183 } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= bnk->mu2_i) { 184 /* Good agreement */ 185 max_radius = PetscMax(max_radius, tao->trust); 186 187 if (tau_max < bnk->gamma2_i) { 188 tau = bnk->gamma2_i; 189 } else if (tau_max > bnk->gamma3_i) { 190 tau = bnk->gamma3_i; 191 } else { 192 tau = tau_max; 193 } 194 } else { 195 /* Not good agreement */ 196 if (tau_min > 1.0) { 197 tau = bnk->gamma2_i; 198 } else if (tau_max < bnk->gamma1_i) { 199 tau = bnk->gamma1_i; 200 } else if ((tau_min < bnk->gamma1_i) && (tau_max >= 1.0)) { 201 tau = bnk->gamma1_i; 202 } else if ((tau_1 >= bnk->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < bnk->gamma1_i) || (tau_2 >= 1.0))) { 203 tau = tau_1; 204 } else if ((tau_2 >= bnk->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < bnk->gamma1_i) || (tau_2 >= 1.0))) { 205 tau = tau_2; 206 } else { 207 tau = tau_max; 208 } 209 } 210 } 211 tao->trust = tau * tao->trust; 212 } 213 214 if (f_min < bnk->f) { 215 /* We accidentally found a solution better than the initial, so accept it */ 216 bnk->f = f_min; 217 PetscCall(VecCopy(tao->solution, bnk->Xold)); 218 PetscCall(VecAXPY(tao->solution, sigma, tao->gradient)); 219 PetscCall(TaoBoundSolution(tao->solution, tao->XL, tao->XU, 0.0, &nDiff, tao->solution)); 220 PetscCall(VecCopy(tao->solution, tao->stepdirection)); 221 PetscCall(VecAXPY(tao->stepdirection, -1.0, bnk->Xold)); 222 PetscCall(TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient)); 223 PetscCall(TaoBNKEstimateActiveSet(tao, bnk->as_type)); 224 PetscCall(VecCopy(bnk->unprojected_gradient, tao->gradient)); 225 PetscCall(VecISSet(tao->gradient, bnk->active_idx, 0.0)); 226 /* Compute gradient at the new iterate and flip switch to compute the Hessian later */ 227 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm)); 228 *needH = PETSC_TRUE; 229 /* Test the new step for convergence */ 230 PetscCall(VecFischer(tao->solution, bnk->unprojected_gradient, tao->XL, tao->XU, bnk->W)); 231 PetscCall(VecNorm(bnk->W, NORM_2, &resnorm)); 232 PetscCheck(!PetscIsInfOrNanReal(resnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 233 PetscCall(TaoLogConvergenceHistory(tao, bnk->f, resnorm, 0.0, tao->ksp_its)); 234 PetscCall(TaoMonitor(tao, tao->niter, bnk->f, resnorm, 0.0, 1.0)); 235 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 236 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 237 /* active BNCG recycling early because we have a stepdirection computed */ 238 PetscCall(TaoSetRecycleHistory(bnk->bncg, PETSC_TRUE)); 239 } 240 } 241 tao->trust = PetscMax(tao->trust, max_radius); 242 243 /* Ensure that the trust radius is within the limits */ 244 tao->trust = PetscMax(tao->trust, bnk->min_radius); 245 tao->trust = PetscMin(tao->trust, bnk->max_radius); 246 break; 247 248 default: 249 /* Norm of the first direction will initialize radius */ 250 tao->trust = 0.0; 251 break; 252 } 253 } 254 PetscFunctionReturn(0); 255 } 256 257 /*------------------------------------------------------------*/ 258 259 /* Routine for computing the exact Hessian and preparing the preconditioner at the new iterate */ 260 261 PetscErrorCode TaoBNKComputeHessian(Tao tao) 262 { 263 TAO_BNK *bnk = (TAO_BNK *)tao->data; 264 265 PetscFunctionBegin; 266 /* Compute the Hessian */ 267 PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre)); 268 /* Add a correction to the BFGS preconditioner */ 269 if (bnk->M) PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 270 /* Prepare the reduced sub-matrices for the inactive set */ 271 PetscCall(MatDestroy(&bnk->Hpre_inactive)); 272 PetscCall(MatDestroy(&bnk->H_inactive)); 273 if (bnk->active_idx) { 274 PetscCall(MatCreateSubMatrix(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, MAT_INITIAL_MATRIX, &bnk->H_inactive)); 275 if (tao->hessian == tao->hessian_pre) { 276 PetscCall(PetscObjectReference((PetscObject)bnk->H_inactive)); 277 bnk->Hpre_inactive = bnk->H_inactive; 278 } else { 279 PetscCall(MatCreateSubMatrix(tao->hessian_pre, bnk->inactive_idx, bnk->inactive_idx, MAT_INITIAL_MATRIX, &bnk->Hpre_inactive)); 280 } 281 if (bnk->bfgs_pre) PetscCall(PCLMVMSetIS(bnk->bfgs_pre, bnk->inactive_idx)); 282 } else { 283 PetscCall(PetscObjectReference((PetscObject)tao->hessian)); 284 bnk->H_inactive = tao->hessian; 285 if (tao->hessian == tao->hessian_pre) { 286 PetscCall(PetscObjectReference((PetscObject)bnk->H_inactive)); 287 bnk->Hpre_inactive = bnk->H_inactive; 288 } else { 289 PetscCall(PetscObjectReference((PetscObject)tao->hessian_pre)); 290 bnk->Hpre_inactive = tao->hessian_pre; 291 } 292 if (bnk->bfgs_pre) PetscCall(PCLMVMClearIS(bnk->bfgs_pre)); 293 } 294 PetscFunctionReturn(0); 295 } 296 297 /*------------------------------------------------------------*/ 298 299 /* Routine for estimating the active set */ 300 301 PetscErrorCode TaoBNKEstimateActiveSet(Tao tao, PetscInt asType) 302 { 303 TAO_BNK *bnk = (TAO_BNK *)tao->data; 304 PetscBool hessComputed, diagExists, hadactive; 305 306 PetscFunctionBegin; 307 hadactive = bnk->active_idx ? PETSC_TRUE : PETSC_FALSE; 308 switch (asType) { 309 case BNK_AS_NONE: 310 PetscCall(ISDestroy(&bnk->inactive_idx)); 311 PetscCall(VecWhichInactive(tao->XL, tao->solution, bnk->unprojected_gradient, tao->XU, PETSC_TRUE, &bnk->inactive_idx)); 312 PetscCall(ISDestroy(&bnk->active_idx)); 313 PetscCall(ISComplementVec(bnk->inactive_idx, tao->solution, &bnk->active_idx)); 314 break; 315 316 case BNK_AS_BERTSEKAS: 317 /* Compute the trial step vector with which we will estimate the active set at the next iteration */ 318 if (bnk->M) { 319 /* If the BFGS preconditioner matrix is available, we will construct a trial step with it */ 320 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, bnk->W)); 321 } else { 322 hessComputed = diagExists = PETSC_FALSE; 323 if (tao->hessian) PetscCall(MatAssembled(tao->hessian, &hessComputed)); 324 if (hessComputed) PetscCall(MatHasOperation(tao->hessian, MATOP_GET_DIAGONAL, &diagExists)); 325 if (diagExists) { 326 /* BFGS preconditioner doesn't exist so let's invert the absolute diagonal of the Hessian instead onto the gradient */ 327 PetscCall(MatGetDiagonal(tao->hessian, bnk->Xwork)); 328 PetscCall(VecAbs(bnk->Xwork)); 329 PetscCall(VecMedian(bnk->Diag_min, bnk->Xwork, bnk->Diag_max, bnk->Xwork)); 330 PetscCall(VecReciprocal(bnk->Xwork)); 331 PetscCall(VecPointwiseMult(bnk->W, bnk->Xwork, bnk->unprojected_gradient)); 332 } else { 333 /* If the Hessian or its diagonal does not exist, we will simply use gradient step */ 334 PetscCall(VecCopy(bnk->unprojected_gradient, bnk->W)); 335 } 336 } 337 PetscCall(VecScale(bnk->W, -1.0)); 338 PetscCall(TaoEstimateActiveBounds(tao->solution, tao->XL, tao->XU, bnk->unprojected_gradient, bnk->W, bnk->Xwork, bnk->as_step, &bnk->as_tol, &bnk->active_lower, &bnk->active_upper, &bnk->active_fixed, &bnk->active_idx, &bnk->inactive_idx)); 339 break; 340 341 default: 342 break; 343 } 344 bnk->resetksp = (PetscBool)(bnk->active_idx || hadactive); /* inactive Hessian size may have changed, need to reset operators */ 345 PetscFunctionReturn(0); 346 } 347 348 /*------------------------------------------------------------*/ 349 350 /* Routine for bounding the step direction */ 351 352 PetscErrorCode TaoBNKBoundStep(Tao tao, PetscInt asType, Vec step) 353 { 354 TAO_BNK *bnk = (TAO_BNK *)tao->data; 355 356 PetscFunctionBegin; 357 switch (asType) { 358 case BNK_AS_NONE: 359 PetscCall(VecISSet(step, bnk->active_idx, 0.0)); 360 break; 361 362 case BNK_AS_BERTSEKAS: 363 PetscCall(TaoBoundStep(tao->solution, tao->XL, tao->XU, bnk->active_lower, bnk->active_upper, bnk->active_fixed, 1.0, step)); 364 break; 365 366 default: 367 break; 368 } 369 PetscFunctionReturn(0); 370 } 371 372 /*------------------------------------------------------------*/ 373 374 /* Routine for taking a finite number of BNCG iterations to 375 accelerate Newton convergence. 376 377 In practice, this approach simply trades off Hessian evaluations 378 for more gradient evaluations. 379 */ 380 381 PetscErrorCode TaoBNKTakeCGSteps(Tao tao, PetscBool *terminate) 382 { 383 TAO_BNK *bnk = (TAO_BNK *)tao->data; 384 385 PetscFunctionBegin; 386 *terminate = PETSC_FALSE; 387 if (bnk->max_cg_its > 0) { 388 /* Copy the current function value (important vectors are already shared) */ 389 bnk->bncg_ctx->f = bnk->f; 390 /* Take some small finite number of BNCG iterations */ 391 PetscCall(TaoSolve(bnk->bncg)); 392 /* Add the number of gradient and function evaluations to the total */ 393 tao->nfuncs += bnk->bncg->nfuncs; 394 tao->nfuncgrads += bnk->bncg->nfuncgrads; 395 tao->ngrads += bnk->bncg->ngrads; 396 tao->nhess += bnk->bncg->nhess; 397 bnk->tot_cg_its += bnk->bncg->niter; 398 /* Extract the BNCG function value out and save it into BNK */ 399 bnk->f = bnk->bncg_ctx->f; 400 if (bnk->bncg->reason == TAO_CONVERGED_GATOL || bnk->bncg->reason == TAO_CONVERGED_GRTOL || bnk->bncg->reason == TAO_CONVERGED_GTTOL || bnk->bncg->reason == TAO_CONVERGED_MINF) { 401 *terminate = PETSC_TRUE; 402 } else { 403 PetscCall(TaoBNKEstimateActiveSet(tao, bnk->as_type)); 404 } 405 } 406 PetscFunctionReturn(0); 407 } 408 409 /*------------------------------------------------------------*/ 410 411 /* Routine for computing the Newton step. */ 412 413 PetscErrorCode TaoBNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type) 414 { 415 TAO_BNK *bnk = (TAO_BNK *)tao->data; 416 PetscInt bfgsUpdates = 0; 417 PetscInt kspits; 418 PetscBool is_lmvm; 419 PetscVoidFunction kspTR; 420 421 PetscFunctionBegin; 422 /* If there are no inactive variables left, save some computation and return an adjusted zero step 423 that has (l-x) and (u-x) for lower and upper bounded variables. */ 424 if (!bnk->inactive_idx) { 425 PetscCall(VecSet(tao->stepdirection, 0.0)); 426 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 427 PetscFunctionReturn(0); 428 } 429 430 /* Shift the reduced Hessian matrix */ 431 if (shift && bnk->pert > 0) { 432 PetscCall(PetscObjectTypeCompare((PetscObject)tao->hessian, MATLMVM, &is_lmvm)); 433 if (is_lmvm) { 434 PetscCall(MatShift(tao->hessian, bnk->pert)); 435 } else { 436 PetscCall(MatShift(bnk->H_inactive, bnk->pert)); 437 if (bnk->H_inactive != bnk->Hpre_inactive) PetscCall(MatShift(bnk->Hpre_inactive, bnk->pert)); 438 } 439 } 440 441 /* Solve the Newton system of equations */ 442 tao->ksp_its = 0; 443 PetscCall(VecSet(tao->stepdirection, 0.0)); 444 if (bnk->resetksp) { 445 PetscCall(KSPReset(tao->ksp)); 446 PetscCall(KSPResetFromOptions(tao->ksp)); 447 bnk->resetksp = PETSC_FALSE; 448 } 449 PetscCall(KSPSetOperators(tao->ksp, bnk->H_inactive, bnk->Hpre_inactive)); 450 PetscCall(VecCopy(bnk->unprojected_gradient, bnk->Gwork)); 451 if (bnk->active_idx) { 452 PetscCall(VecGetSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive)); 453 PetscCall(VecGetSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive)); 454 } else { 455 bnk->G_inactive = bnk->unprojected_gradient; 456 bnk->X_inactive = tao->stepdirection; 457 } 458 PetscCall(KSPCGSetRadius(tao->ksp, tao->trust)); 459 PetscCall(KSPSolve(tao->ksp, bnk->G_inactive, bnk->X_inactive)); 460 PetscCall(KSPGetIterationNumber(tao->ksp, &kspits)); 461 tao->ksp_its += kspits; 462 tao->ksp_tot_its += kspits; 463 PetscCall(PetscObjectQueryFunction((PetscObject)tao->ksp, "KSPCGGetNormD_C", &kspTR)); 464 if (kspTR) { 465 PetscCall(KSPCGGetNormD(tao->ksp, &bnk->dnorm)); 466 467 if (0.0 == tao->trust) { 468 /* Radius was uninitialized; use the norm of the direction */ 469 if (bnk->dnorm > 0.0) { 470 tao->trust = bnk->dnorm; 471 472 /* Modify the radius if it is too large or small */ 473 tao->trust = PetscMax(tao->trust, bnk->min_radius); 474 tao->trust = PetscMin(tao->trust, bnk->max_radius); 475 } else { 476 /* The direction was bad; set radius to default value and re-solve 477 the trust-region subproblem to get a direction */ 478 tao->trust = tao->trust0; 479 480 /* Modify the radius if it is too large or small */ 481 tao->trust = PetscMax(tao->trust, bnk->min_radius); 482 tao->trust = PetscMin(tao->trust, bnk->max_radius); 483 484 PetscCall(KSPCGSetRadius(tao->ksp, tao->trust)); 485 PetscCall(KSPSolve(tao->ksp, bnk->G_inactive, bnk->X_inactive)); 486 PetscCall(KSPGetIterationNumber(tao->ksp, &kspits)); 487 tao->ksp_its += kspits; 488 tao->ksp_tot_its += kspits; 489 PetscCall(KSPCGGetNormD(tao->ksp, &bnk->dnorm)); 490 491 PetscCheck(bnk->dnorm != 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_PLIB, "Initial direction zero"); 492 } 493 } 494 } 495 /* Restore sub vectors back */ 496 if (bnk->active_idx) { 497 PetscCall(VecRestoreSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive)); 498 PetscCall(VecRestoreSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive)); 499 } 500 /* Make sure the safeguarded fall-back step is zero for actively bounded variables */ 501 PetscCall(VecScale(tao->stepdirection, -1.0)); 502 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 503 504 /* Record convergence reasons */ 505 PetscCall(KSPGetConvergedReason(tao->ksp, ksp_reason)); 506 if (KSP_CONVERGED_ATOL == *ksp_reason) { 507 ++bnk->ksp_atol; 508 } else if (KSP_CONVERGED_RTOL == *ksp_reason) { 509 ++bnk->ksp_rtol; 510 } else if (KSP_CONVERGED_CG_CONSTRAINED == *ksp_reason) { 511 ++bnk->ksp_ctol; 512 } else if (KSP_CONVERGED_CG_NEG_CURVE == *ksp_reason) { 513 ++bnk->ksp_negc; 514 } else if (KSP_DIVERGED_DTOL == *ksp_reason) { 515 ++bnk->ksp_dtol; 516 } else if (KSP_DIVERGED_ITS == *ksp_reason) { 517 ++bnk->ksp_iter; 518 } else { 519 ++bnk->ksp_othr; 520 } 521 522 /* Make sure the BFGS preconditioner is healthy */ 523 if (bnk->M) { 524 PetscCall(MatLMVMGetUpdateCount(bnk->M, &bfgsUpdates)); 525 if ((KSP_DIVERGED_INDEFINITE_PC == *ksp_reason) && (bfgsUpdates > 0)) { 526 /* Preconditioner is numerically indefinite; reset the approximation. */ 527 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 528 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 529 } 530 } 531 *step_type = BNK_NEWTON; 532 PetscFunctionReturn(0); 533 } 534 535 /*------------------------------------------------------------*/ 536 537 /* Routine for recomputing the predicted reduction for a given step vector */ 538 539 PetscErrorCode TaoBNKRecomputePred(Tao tao, Vec S, PetscReal *prered) 540 { 541 TAO_BNK *bnk = (TAO_BNK *)tao->data; 542 543 PetscFunctionBegin; 544 /* Extract subvectors associated with the inactive set */ 545 if (bnk->active_idx) { 546 PetscCall(VecGetSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive)); 547 PetscCall(VecGetSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work)); 548 PetscCall(VecGetSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive)); 549 } else { 550 bnk->X_inactive = tao->stepdirection; 551 bnk->inactive_work = bnk->Xwork; 552 bnk->G_inactive = bnk->Gwork; 553 } 554 /* Recompute the predicted decrease based on the quadratic model */ 555 PetscCall(MatMult(bnk->H_inactive, bnk->X_inactive, bnk->inactive_work)); 556 PetscCall(VecAYPX(bnk->inactive_work, -0.5, bnk->G_inactive)); 557 PetscCall(VecDot(bnk->inactive_work, bnk->X_inactive, prered)); 558 /* Restore the sub vectors */ 559 if (bnk->active_idx) { 560 PetscCall(VecRestoreSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive)); 561 PetscCall(VecRestoreSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work)); 562 PetscCall(VecRestoreSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive)); 563 } 564 PetscFunctionReturn(0); 565 } 566 567 /*------------------------------------------------------------*/ 568 569 /* Routine for ensuring that the Newton step is a descent direction. 570 571 The step direction falls back onto BFGS, scaled gradient and gradient steps 572 in the event that the Newton step fails the test. 573 */ 574 575 PetscErrorCode TaoBNKSafeguardStep(Tao tao, KSPConvergedReason ksp_reason, PetscInt *stepType) 576 { 577 TAO_BNK *bnk = (TAO_BNK *)tao->data; 578 PetscReal gdx, e_min; 579 PetscInt bfgsUpdates; 580 581 PetscFunctionBegin; 582 switch (*stepType) { 583 case BNK_NEWTON: 584 PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx)); 585 if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) { 586 /* Newton step is not descent or direction produced Inf or NaN 587 Update the perturbation for next time */ 588 if (bnk->pert <= 0.0) { 589 PetscBool is_gltr; 590 591 /* Initialize the perturbation */ 592 bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm)); 593 PetscCall(PetscObjectTypeCompare((PetscObject)(tao->ksp), KSPGLTR, &is_gltr)); 594 if (is_gltr) { 595 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 596 bnk->pert = PetscMax(bnk->pert, -e_min); 597 } 598 } else { 599 /* Increase the perturbation */ 600 bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm)); 601 } 602 603 if (!bnk->M) { 604 /* We don't have the bfgs matrix around and updated 605 Must use gradient direction in this case */ 606 PetscCall(VecCopy(tao->gradient, tao->stepdirection)); 607 *stepType = BNK_GRADIENT; 608 } else { 609 /* Attempt to use the BFGS direction */ 610 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 611 612 /* Check for success (descent direction) 613 NOTE: Negative gdx here means not a descent direction because 614 the fall-back step is missing a negative sign. */ 615 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx)); 616 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 617 /* BFGS direction is not descent or direction produced not a number 618 We can assert bfgsUpdates > 1 in this case because 619 the first solve produces the scaled gradient direction, 620 which is guaranteed to be descent */ 621 622 /* Use steepest descent direction (scaled) */ 623 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 624 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 625 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 626 627 *stepType = BNK_SCALED_GRADIENT; 628 } else { 629 PetscCall(MatLMVMGetUpdateCount(bnk->M, &bfgsUpdates)); 630 if (1 == bfgsUpdates) { 631 /* The first BFGS direction is always the scaled gradient */ 632 *stepType = BNK_SCALED_GRADIENT; 633 } else { 634 *stepType = BNK_BFGS; 635 } 636 } 637 } 638 /* Make sure the safeguarded fall-back step is zero for actively bounded variables */ 639 PetscCall(VecScale(tao->stepdirection, -1.0)); 640 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 641 } else { 642 /* Computed Newton step is descent */ 643 switch (ksp_reason) { 644 case KSP_DIVERGED_NANORINF: 645 case KSP_DIVERGED_BREAKDOWN: 646 case KSP_DIVERGED_INDEFINITE_MAT: 647 case KSP_DIVERGED_INDEFINITE_PC: 648 case KSP_CONVERGED_CG_NEG_CURVE: 649 /* Matrix or preconditioner is indefinite; increase perturbation */ 650 if (bnk->pert <= 0.0) { 651 PetscBool is_gltr; 652 653 /* Initialize the perturbation */ 654 bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm)); 655 PetscCall(PetscObjectTypeCompare((PetscObject)(tao->ksp), KSPGLTR, &is_gltr)); 656 if (is_gltr) { 657 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 658 bnk->pert = PetscMax(bnk->pert, -e_min); 659 } 660 } else { 661 /* Increase the perturbation */ 662 bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm)); 663 } 664 break; 665 666 default: 667 /* Newton step computation is good; decrease perturbation */ 668 bnk->pert = PetscMin(bnk->psfac * bnk->pert, bnk->pmsfac * bnk->gnorm); 669 if (bnk->pert < bnk->pmin) bnk->pert = 0.0; 670 break; 671 } 672 *stepType = BNK_NEWTON; 673 } 674 break; 675 676 case BNK_BFGS: 677 /* Check for success (descent direction) */ 678 PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx)); 679 if (gdx >= 0 || PetscIsInfOrNanReal(gdx)) { 680 /* Step is not descent or solve was not successful 681 Use steepest descent direction (scaled) */ 682 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 683 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 684 PetscCall(MatSolve(bnk->M, tao->gradient, tao->stepdirection)); 685 PetscCall(VecScale(tao->stepdirection, -1.0)); 686 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 687 *stepType = BNK_SCALED_GRADIENT; 688 } else { 689 *stepType = BNK_BFGS; 690 } 691 break; 692 693 case BNK_SCALED_GRADIENT: 694 break; 695 696 default: 697 break; 698 } 699 700 PetscFunctionReturn(0); 701 } 702 703 /*------------------------------------------------------------*/ 704 705 /* Routine for performing a bound-projected More-Thuente line search. 706 707 Includes fallbacks to BFGS, scaled gradient, and unscaled gradient steps if the 708 Newton step does not produce a valid step length. 709 */ 710 711 PetscErrorCode TaoBNKPerformLineSearch(Tao tao, PetscInt *stepType, PetscReal *steplen, TaoLineSearchConvergedReason *reason) 712 { 713 TAO_BNK *bnk = (TAO_BNK *)tao->data; 714 TaoLineSearchConvergedReason ls_reason; 715 PetscReal e_min, gdx; 716 PetscInt bfgsUpdates; 717 718 PetscFunctionBegin; 719 /* Perform the linesearch */ 720 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, bnk->unprojected_gradient, tao->stepdirection, steplen, &ls_reason)); 721 PetscCall(TaoAddLineSearchCounts(tao)); 722 723 while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && *stepType != BNK_SCALED_GRADIENT && *stepType != BNK_GRADIENT) { 724 /* Linesearch failed, revert solution */ 725 bnk->f = bnk->fold; 726 PetscCall(VecCopy(bnk->Xold, tao->solution)); 727 PetscCall(VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient)); 728 729 switch (*stepType) { 730 case BNK_NEWTON: 731 /* Failed to obtain acceptable iterate with Newton step 732 Update the perturbation for next time */ 733 if (bnk->pert <= 0.0) { 734 PetscBool is_gltr; 735 736 /* Initialize the perturbation */ 737 bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm)); 738 PetscCall(PetscObjectTypeCompare((PetscObject)(tao->ksp), KSPGLTR, &is_gltr)); 739 if (is_gltr) { 740 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 741 bnk->pert = PetscMax(bnk->pert, -e_min); 742 } 743 } else { 744 /* Increase the perturbation */ 745 bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm)); 746 } 747 748 if (!bnk->M) { 749 /* We don't have the bfgs matrix around and being updated 750 Must use gradient direction in this case */ 751 PetscCall(VecCopy(bnk->unprojected_gradient, tao->stepdirection)); 752 *stepType = BNK_GRADIENT; 753 } else { 754 /* Attempt to use the BFGS direction */ 755 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 756 /* Check for success (descent direction) 757 NOTE: Negative gdx means not a descent direction because the step here is missing a negative sign. */ 758 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx)); 759 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 760 /* BFGS direction is not descent or direction produced not a number 761 We can assert bfgsUpdates > 1 in this case 762 Use steepest descent direction (scaled) */ 763 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 764 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 765 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 766 767 bfgsUpdates = 1; 768 *stepType = BNK_SCALED_GRADIENT; 769 } else { 770 PetscCall(MatLMVMGetUpdateCount(bnk->M, &bfgsUpdates)); 771 if (1 == bfgsUpdates) { 772 /* The first BFGS direction is always the scaled gradient */ 773 *stepType = BNK_SCALED_GRADIENT; 774 } else { 775 *stepType = BNK_BFGS; 776 } 777 } 778 } 779 break; 780 781 case BNK_BFGS: 782 /* Can only enter if pc_type == BNK_PC_BFGS 783 Failed to obtain acceptable iterate with BFGS step 784 Attempt to use the scaled gradient direction */ 785 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 786 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 787 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 788 789 bfgsUpdates = 1; 790 *stepType = BNK_SCALED_GRADIENT; 791 break; 792 } 793 /* Make sure the safeguarded fall-back step is zero for actively bounded variables */ 794 PetscCall(VecScale(tao->stepdirection, -1.0)); 795 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 796 797 /* Perform one last line search with the fall-back step */ 798 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, bnk->unprojected_gradient, tao->stepdirection, steplen, &ls_reason)); 799 PetscCall(TaoAddLineSearchCounts(tao)); 800 } 801 *reason = ls_reason; 802 PetscFunctionReturn(0); 803 } 804 805 /*------------------------------------------------------------*/ 806 807 /* Routine for updating the trust radius. 808 809 Function features three different update methods: 810 1) Line-search step length based 811 2) Predicted decrease on the CG quadratic model 812 3) Interpolation 813 */ 814 815 PetscErrorCode TaoBNKUpdateTrustRadius(Tao tao, PetscReal prered, PetscReal actred, PetscInt updateType, PetscInt stepType, PetscBool *accept) 816 { 817 TAO_BNK *bnk = (TAO_BNK *)tao->data; 818 819 PetscReal step, kappa; 820 PetscReal gdx, tau_1, tau_2, tau_min, tau_max; 821 822 PetscFunctionBegin; 823 /* Update trust region radius */ 824 *accept = PETSC_FALSE; 825 switch (updateType) { 826 case BNK_UPDATE_STEP: 827 *accept = PETSC_TRUE; /* always accept here because line search succeeded */ 828 if (stepType == BNK_NEWTON) { 829 PetscCall(TaoLineSearchGetStepLength(tao->linesearch, &step)); 830 if (step < bnk->nu1) { 831 /* Very bad step taken; reduce radius */ 832 tao->trust = bnk->omega1 * PetscMin(bnk->dnorm, tao->trust); 833 } else if (step < bnk->nu2) { 834 /* Reasonably bad step taken; reduce radius */ 835 tao->trust = bnk->omega2 * PetscMin(bnk->dnorm, tao->trust); 836 } else if (step < bnk->nu3) { 837 /* Reasonable step was taken; leave radius alone */ 838 if (bnk->omega3 < 1.0) { 839 tao->trust = bnk->omega3 * PetscMin(bnk->dnorm, tao->trust); 840 } else if (bnk->omega3 > 1.0) { 841 tao->trust = PetscMax(bnk->omega3 * bnk->dnorm, tao->trust); 842 } 843 } else if (step < bnk->nu4) { 844 /* Full step taken; increase the radius */ 845 tao->trust = PetscMax(bnk->omega4 * bnk->dnorm, tao->trust); 846 } else { 847 /* More than full step taken; increase the radius */ 848 tao->trust = PetscMax(bnk->omega5 * bnk->dnorm, tao->trust); 849 } 850 } else { 851 /* Newton step was not good; reduce the radius */ 852 tao->trust = bnk->omega1 * PetscMin(bnk->dnorm, tao->trust); 853 } 854 break; 855 856 case BNK_UPDATE_REDUCTION: 857 if (stepType == BNK_NEWTON) { 858 if ((prered < 0.0) || PetscIsInfOrNanReal(prered)) { 859 /* The predicted reduction has the wrong sign. This cannot 860 happen in infinite precision arithmetic. Step should 861 be rejected! */ 862 tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm); 863 } else { 864 if (PetscIsInfOrNanReal(actred)) { 865 tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm); 866 } else { 867 if ((PetscAbsScalar(actred) <= PetscMax(1.0, PetscAbsScalar(bnk->f)) * bnk->epsilon) && (PetscAbsScalar(prered) <= PetscMax(1.0, PetscAbsScalar(bnk->f)) * bnk->epsilon)) { 868 kappa = 1.0; 869 } else { 870 kappa = actred / prered; 871 } 872 /* Accept or reject the step and update radius */ 873 if (kappa < bnk->eta1) { 874 /* Reject the step */ 875 tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm); 876 } else { 877 /* Accept the step */ 878 *accept = PETSC_TRUE; 879 /* Update the trust region radius only if the computed step is at the trust radius boundary */ 880 if (bnk->dnorm == tao->trust) { 881 if (kappa < bnk->eta2) { 882 /* Marginal bad step */ 883 tao->trust = bnk->alpha2 * tao->trust; 884 } else if (kappa < bnk->eta3) { 885 /* Reasonable step */ 886 tao->trust = bnk->alpha3 * tao->trust; 887 } else if (kappa < bnk->eta4) { 888 /* Good step */ 889 tao->trust = bnk->alpha4 * tao->trust; 890 } else { 891 /* Very good step */ 892 tao->trust = bnk->alpha5 * tao->trust; 893 } 894 } 895 } 896 } 897 } 898 } else { 899 /* Newton step was not good; reduce the radius */ 900 tao->trust = bnk->alpha1 * PetscMin(bnk->dnorm, tao->trust); 901 } 902 break; 903 904 default: 905 if (stepType == BNK_NEWTON) { 906 if (prered < 0.0) { 907 /* The predicted reduction has the wrong sign. This cannot */ 908 /* happen in infinite precision arithmetic. Step should */ 909 /* be rejected! */ 910 tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm); 911 } else { 912 if (PetscIsInfOrNanReal(actred)) { 913 tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm); 914 } else { 915 if ((PetscAbsScalar(actred) <= bnk->epsilon) && (PetscAbsScalar(prered) <= bnk->epsilon)) { 916 kappa = 1.0; 917 } else { 918 kappa = actred / prered; 919 } 920 921 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx)); 922 tau_1 = bnk->theta * gdx / (bnk->theta * gdx - (1.0 - bnk->theta) * prered + actred); 923 tau_2 = bnk->theta * gdx / (bnk->theta * gdx + (1.0 + bnk->theta) * prered - actred); 924 tau_min = PetscMin(tau_1, tau_2); 925 tau_max = PetscMax(tau_1, tau_2); 926 927 if (kappa >= 1.0 - bnk->mu1) { 928 /* Great agreement */ 929 *accept = PETSC_TRUE; 930 if (tau_max < 1.0) { 931 tao->trust = PetscMax(tao->trust, bnk->gamma3 * bnk->dnorm); 932 } else if (tau_max > bnk->gamma4) { 933 tao->trust = PetscMax(tao->trust, bnk->gamma4 * bnk->dnorm); 934 } else { 935 tao->trust = PetscMax(tao->trust, tau_max * bnk->dnorm); 936 } 937 } else if (kappa >= 1.0 - bnk->mu2) { 938 /* Good agreement */ 939 *accept = PETSC_TRUE; 940 if (tau_max < bnk->gamma2) { 941 tao->trust = bnk->gamma2 * PetscMin(tao->trust, bnk->dnorm); 942 } else if (tau_max > bnk->gamma3) { 943 tao->trust = PetscMax(tao->trust, bnk->gamma3 * bnk->dnorm); 944 } else if (tau_max < 1.0) { 945 tao->trust = tau_max * PetscMin(tao->trust, bnk->dnorm); 946 } else { 947 tao->trust = PetscMax(tao->trust, tau_max * bnk->dnorm); 948 } 949 } else { 950 /* Not good agreement */ 951 if (tau_min > 1.0) { 952 tao->trust = bnk->gamma2 * PetscMin(tao->trust, bnk->dnorm); 953 } else if (tau_max < bnk->gamma1) { 954 tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm); 955 } else if ((tau_min < bnk->gamma1) && (tau_max >= 1.0)) { 956 tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm); 957 } else if ((tau_1 >= bnk->gamma1) && (tau_1 < 1.0) && ((tau_2 < bnk->gamma1) || (tau_2 >= 1.0))) { 958 tao->trust = tau_1 * PetscMin(tao->trust, bnk->dnorm); 959 } else if ((tau_2 >= bnk->gamma1) && (tau_2 < 1.0) && ((tau_1 < bnk->gamma1) || (tau_2 >= 1.0))) { 960 tao->trust = tau_2 * PetscMin(tao->trust, bnk->dnorm); 961 } else { 962 tao->trust = tau_max * PetscMin(tao->trust, bnk->dnorm); 963 } 964 } 965 } 966 } 967 } else { 968 /* Newton step was not good; reduce the radius */ 969 tao->trust = bnk->gamma1 * PetscMin(bnk->dnorm, tao->trust); 970 } 971 break; 972 } 973 /* Make sure the radius does not violate min and max settings */ 974 tao->trust = PetscMin(tao->trust, bnk->max_radius); 975 tao->trust = PetscMax(tao->trust, bnk->min_radius); 976 PetscFunctionReturn(0); 977 } 978 979 /* ---------------------------------------------------------- */ 980 981 PetscErrorCode TaoBNKAddStepCounts(Tao tao, PetscInt stepType) 982 { 983 TAO_BNK *bnk = (TAO_BNK *)tao->data; 984 985 PetscFunctionBegin; 986 switch (stepType) { 987 case BNK_NEWTON: 988 ++bnk->newt; 989 break; 990 case BNK_BFGS: 991 ++bnk->bfgs; 992 break; 993 case BNK_SCALED_GRADIENT: 994 ++bnk->sgrad; 995 break; 996 case BNK_GRADIENT: 997 ++bnk->grad; 998 break; 999 default: 1000 break; 1001 } 1002 PetscFunctionReturn(0); 1003 } 1004 1005 /* ---------------------------------------------------------- */ 1006 1007 PetscErrorCode TaoSetUp_BNK(Tao tao) 1008 { 1009 TAO_BNK *bnk = (TAO_BNK *)tao->data; 1010 PetscInt i; 1011 1012 PetscFunctionBegin; 1013 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 1014 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 1015 if (!bnk->W) PetscCall(VecDuplicate(tao->solution, &bnk->W)); 1016 if (!bnk->Xold) PetscCall(VecDuplicate(tao->solution, &bnk->Xold)); 1017 if (!bnk->Gold) PetscCall(VecDuplicate(tao->solution, &bnk->Gold)); 1018 if (!bnk->Xwork) PetscCall(VecDuplicate(tao->solution, &bnk->Xwork)); 1019 if (!bnk->Gwork) PetscCall(VecDuplicate(tao->solution, &bnk->Gwork)); 1020 if (!bnk->unprojected_gradient) PetscCall(VecDuplicate(tao->solution, &bnk->unprojected_gradient)); 1021 if (!bnk->unprojected_gradient_old) PetscCall(VecDuplicate(tao->solution, &bnk->unprojected_gradient_old)); 1022 if (!bnk->Diag_min) PetscCall(VecDuplicate(tao->solution, &bnk->Diag_min)); 1023 if (!bnk->Diag_max) PetscCall(VecDuplicate(tao->solution, &bnk->Diag_max)); 1024 if (bnk->max_cg_its > 0) { 1025 /* Ensure that the important common vectors are shared between BNK and embedded BNCG */ 1026 bnk->bncg_ctx = (TAO_BNCG *)bnk->bncg->data; 1027 PetscCall(PetscObjectReference((PetscObject)(bnk->unprojected_gradient_old))); 1028 PetscCall(VecDestroy(&bnk->bncg_ctx->unprojected_gradient_old)); 1029 bnk->bncg_ctx->unprojected_gradient_old = bnk->unprojected_gradient_old; 1030 PetscCall(PetscObjectReference((PetscObject)(bnk->unprojected_gradient))); 1031 PetscCall(VecDestroy(&bnk->bncg_ctx->unprojected_gradient)); 1032 bnk->bncg_ctx->unprojected_gradient = bnk->unprojected_gradient; 1033 PetscCall(PetscObjectReference((PetscObject)(bnk->Gold))); 1034 PetscCall(VecDestroy(&bnk->bncg_ctx->G_old)); 1035 bnk->bncg_ctx->G_old = bnk->Gold; 1036 PetscCall(PetscObjectReference((PetscObject)(tao->gradient))); 1037 PetscCall(VecDestroy(&bnk->bncg->gradient)); 1038 bnk->bncg->gradient = tao->gradient; 1039 PetscCall(PetscObjectReference((PetscObject)(tao->stepdirection))); 1040 PetscCall(VecDestroy(&bnk->bncg->stepdirection)); 1041 bnk->bncg->stepdirection = tao->stepdirection; 1042 PetscCall(TaoSetSolution(bnk->bncg, tao->solution)); 1043 /* Copy over some settings from BNK into BNCG */ 1044 PetscCall(TaoSetMaximumIterations(bnk->bncg, bnk->max_cg_its)); 1045 PetscCall(TaoSetTolerances(bnk->bncg, tao->gatol, tao->grtol, tao->gttol)); 1046 PetscCall(TaoSetFunctionLowerBound(bnk->bncg, tao->fmin)); 1047 PetscCall(TaoSetConvergenceTest(bnk->bncg, tao->ops->convergencetest, tao->cnvP)); 1048 PetscCall(TaoSetObjective(bnk->bncg, tao->ops->computeobjective, tao->user_objP)); 1049 PetscCall(TaoSetGradient(bnk->bncg, NULL, tao->ops->computegradient, tao->user_gradP)); 1050 PetscCall(TaoSetObjectiveAndGradient(bnk->bncg, NULL, tao->ops->computeobjectiveandgradient, tao->user_objgradP)); 1051 PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)tao, (PetscObject)(bnk->bncg))); 1052 for (i = 0; i < tao->numbermonitors; ++i) { 1053 PetscCall(TaoSetMonitor(bnk->bncg, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i])); 1054 PetscCall(PetscObjectReference((PetscObject)(tao->monitorcontext[i]))); 1055 } 1056 } 1057 bnk->X_inactive = NULL; 1058 bnk->G_inactive = NULL; 1059 bnk->inactive_work = NULL; 1060 bnk->active_work = NULL; 1061 bnk->inactive_idx = NULL; 1062 bnk->active_idx = NULL; 1063 bnk->active_lower = NULL; 1064 bnk->active_upper = NULL; 1065 bnk->active_fixed = NULL; 1066 bnk->M = NULL; 1067 bnk->H_inactive = NULL; 1068 bnk->Hpre_inactive = NULL; 1069 PetscFunctionReturn(0); 1070 } 1071 1072 /*------------------------------------------------------------*/ 1073 1074 PetscErrorCode TaoDestroy_BNK(Tao tao) 1075 { 1076 TAO_BNK *bnk = (TAO_BNK *)tao->data; 1077 1078 PetscFunctionBegin; 1079 PetscCall(VecDestroy(&bnk->W)); 1080 PetscCall(VecDestroy(&bnk->Xold)); 1081 PetscCall(VecDestroy(&bnk->Gold)); 1082 PetscCall(VecDestroy(&bnk->Xwork)); 1083 PetscCall(VecDestroy(&bnk->Gwork)); 1084 PetscCall(VecDestroy(&bnk->unprojected_gradient)); 1085 PetscCall(VecDestroy(&bnk->unprojected_gradient_old)); 1086 PetscCall(VecDestroy(&bnk->Diag_min)); 1087 PetscCall(VecDestroy(&bnk->Diag_max)); 1088 PetscCall(ISDestroy(&bnk->active_lower)); 1089 PetscCall(ISDestroy(&bnk->active_upper)); 1090 PetscCall(ISDestroy(&bnk->active_fixed)); 1091 PetscCall(ISDestroy(&bnk->active_idx)); 1092 PetscCall(ISDestroy(&bnk->inactive_idx)); 1093 PetscCall(MatDestroy(&bnk->Hpre_inactive)); 1094 PetscCall(MatDestroy(&bnk->H_inactive)); 1095 PetscCall(TaoDestroy(&bnk->bncg)); 1096 PetscCall(KSPDestroy(&tao->ksp)); 1097 PetscCall(PetscFree(tao->data)); 1098 PetscFunctionReturn(0); 1099 } 1100 1101 /*------------------------------------------------------------*/ 1102 1103 PetscErrorCode TaoSetFromOptions_BNK(Tao tao, PetscOptionItems *PetscOptionsObject) 1104 { 1105 TAO_BNK *bnk = (TAO_BNK *)tao->data; 1106 1107 PetscFunctionBegin; 1108 PetscOptionsHeadBegin(PetscOptionsObject, "Newton-Krylov method for bound constrained optimization"); 1109 PetscCall(PetscOptionsEList("-tao_bnk_init_type", "radius initialization type", "", BNK_INIT, BNK_INIT_TYPES, BNK_INIT[bnk->init_type], &bnk->init_type, NULL)); 1110 PetscCall(PetscOptionsEList("-tao_bnk_update_type", "radius update type", "", BNK_UPDATE, BNK_UPDATE_TYPES, BNK_UPDATE[bnk->update_type], &bnk->update_type, NULL)); 1111 PetscCall(PetscOptionsEList("-tao_bnk_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, NULL)); 1112 PetscCall(PetscOptionsReal("-tao_bnk_sval", "(developer) Hessian perturbation starting value", "", bnk->sval, &bnk->sval, NULL)); 1113 PetscCall(PetscOptionsReal("-tao_bnk_imin", "(developer) minimum initial Hessian perturbation", "", bnk->imin, &bnk->imin, NULL)); 1114 PetscCall(PetscOptionsReal("-tao_bnk_imax", "(developer) maximum initial Hessian perturbation", "", bnk->imax, &bnk->imax, NULL)); 1115 PetscCall(PetscOptionsReal("-tao_bnk_imfac", "(developer) initial merit factor for Hessian perturbation", "", bnk->imfac, &bnk->imfac, NULL)); 1116 PetscCall(PetscOptionsReal("-tao_bnk_pmin", "(developer) minimum Hessian perturbation", "", bnk->pmin, &bnk->pmin, NULL)); 1117 PetscCall(PetscOptionsReal("-tao_bnk_pmax", "(developer) maximum Hessian perturbation", "", bnk->pmax, &bnk->pmax, NULL)); 1118 PetscCall(PetscOptionsReal("-tao_bnk_pgfac", "(developer) Hessian perturbation growth factor", "", bnk->pgfac, &bnk->pgfac, NULL)); 1119 PetscCall(PetscOptionsReal("-tao_bnk_psfac", "(developer) Hessian perturbation shrink factor", "", bnk->psfac, &bnk->psfac, NULL)); 1120 PetscCall(PetscOptionsReal("-tao_bnk_pmgfac", "(developer) merit growth factor for Hessian perturbation", "", bnk->pmgfac, &bnk->pmgfac, NULL)); 1121 PetscCall(PetscOptionsReal("-tao_bnk_pmsfac", "(developer) merit shrink factor for Hessian perturbation", "", bnk->pmsfac, &bnk->pmsfac, NULL)); 1122 PetscCall(PetscOptionsReal("-tao_bnk_eta1", "(developer) threshold for rejecting step (-tao_bnk_update_type reduction)", "", bnk->eta1, &bnk->eta1, NULL)); 1123 PetscCall(PetscOptionsReal("-tao_bnk_eta2", "(developer) threshold for accepting marginal step (-tao_bnk_update_type reduction)", "", bnk->eta2, &bnk->eta2, NULL)); 1124 PetscCall(PetscOptionsReal("-tao_bnk_eta3", "(developer) threshold for accepting reasonable step (-tao_bnk_update_type reduction)", "", bnk->eta3, &bnk->eta3, NULL)); 1125 PetscCall(PetscOptionsReal("-tao_bnk_eta4", "(developer) threshold for accepting good step (-tao_bnk_update_type reduction)", "", bnk->eta4, &bnk->eta4, NULL)); 1126 PetscCall(PetscOptionsReal("-tao_bnk_alpha1", "(developer) radius reduction factor for rejected step (-tao_bnk_update_type reduction)", "", bnk->alpha1, &bnk->alpha1, NULL)); 1127 PetscCall(PetscOptionsReal("-tao_bnk_alpha2", "(developer) radius reduction factor for marginally accepted bad step (-tao_bnk_update_type reduction)", "", bnk->alpha2, &bnk->alpha2, NULL)); 1128 PetscCall(PetscOptionsReal("-tao_bnk_alpha3", "(developer) radius increase factor for reasonable accepted step (-tao_bnk_update_type reduction)", "", bnk->alpha3, &bnk->alpha3, NULL)); 1129 PetscCall(PetscOptionsReal("-tao_bnk_alpha4", "(developer) radius increase factor for good accepted step (-tao_bnk_update_type reduction)", "", bnk->alpha4, &bnk->alpha4, NULL)); 1130 PetscCall(PetscOptionsReal("-tao_bnk_alpha5", "(developer) radius increase factor for very good accepted step (-tao_bnk_update_type reduction)", "", bnk->alpha5, &bnk->alpha5, NULL)); 1131 PetscCall(PetscOptionsReal("-tao_bnk_nu1", "(developer) threshold for small line-search step length (-tao_bnk_update_type step)", "", bnk->nu1, &bnk->nu1, NULL)); 1132 PetscCall(PetscOptionsReal("-tao_bnk_nu2", "(developer) threshold for reasonable line-search step length (-tao_bnk_update_type step)", "", bnk->nu2, &bnk->nu2, NULL)); 1133 PetscCall(PetscOptionsReal("-tao_bnk_nu3", "(developer) threshold for large line-search step length (-tao_bnk_update_type step)", "", bnk->nu3, &bnk->nu3, NULL)); 1134 PetscCall(PetscOptionsReal("-tao_bnk_nu4", "(developer) threshold for very large line-search step length (-tao_bnk_update_type step)", "", bnk->nu4, &bnk->nu4, NULL)); 1135 PetscCall(PetscOptionsReal("-tao_bnk_omega1", "(developer) radius reduction factor for very small line-search step length (-tao_bnk_update_type step)", "", bnk->omega1, &bnk->omega1, NULL)); 1136 PetscCall(PetscOptionsReal("-tao_bnk_omega2", "(developer) radius reduction factor for small line-search step length (-tao_bnk_update_type step)", "", bnk->omega2, &bnk->omega2, NULL)); 1137 PetscCall(PetscOptionsReal("-tao_bnk_omega3", "(developer) radius factor for decent line-search step length (-tao_bnk_update_type step)", "", bnk->omega3, &bnk->omega3, NULL)); 1138 PetscCall(PetscOptionsReal("-tao_bnk_omega4", "(developer) radius increase factor for large line-search step length (-tao_bnk_update_type step)", "", bnk->omega4, &bnk->omega4, NULL)); 1139 PetscCall(PetscOptionsReal("-tao_bnk_omega5", "(developer) radius increase factor for very large line-search step length (-tao_bnk_update_type step)", "", bnk->omega5, &bnk->omega5, NULL)); 1140 PetscCall(PetscOptionsReal("-tao_bnk_mu1_i", "(developer) threshold for accepting very good step (-tao_bnk_init_type interpolation)", "", bnk->mu1_i, &bnk->mu1_i, NULL)); 1141 PetscCall(PetscOptionsReal("-tao_bnk_mu2_i", "(developer) threshold for accepting good step (-tao_bnk_init_type interpolation)", "", bnk->mu2_i, &bnk->mu2_i, NULL)); 1142 PetscCall(PetscOptionsReal("-tao_bnk_gamma1_i", "(developer) radius reduction factor for rejected very bad step (-tao_bnk_init_type interpolation)", "", bnk->gamma1_i, &bnk->gamma1_i, NULL)); 1143 PetscCall(PetscOptionsReal("-tao_bnk_gamma2_i", "(developer) radius reduction factor for rejected bad step (-tao_bnk_init_type interpolation)", "", bnk->gamma2_i, &bnk->gamma2_i, NULL)); 1144 PetscCall(PetscOptionsReal("-tao_bnk_gamma3_i", "(developer) radius increase factor for accepted good step (-tao_bnk_init_type interpolation)", "", bnk->gamma3_i, &bnk->gamma3_i, NULL)); 1145 PetscCall(PetscOptionsReal("-tao_bnk_gamma4_i", "(developer) radius increase factor for accepted very good step (-tao_bnk_init_type interpolation)", "", bnk->gamma4_i, &bnk->gamma4_i, NULL)); 1146 PetscCall(PetscOptionsReal("-tao_bnk_theta_i", "(developer) trust region interpolation factor (-tao_bnk_init_type interpolation)", "", bnk->theta_i, &bnk->theta_i, NULL)); 1147 PetscCall(PetscOptionsReal("-tao_bnk_mu1", "(developer) threshold for accepting very good step (-tao_bnk_update_type interpolation)", "", bnk->mu1, &bnk->mu1, NULL)); 1148 PetscCall(PetscOptionsReal("-tao_bnk_mu2", "(developer) threshold for accepting good step (-tao_bnk_update_type interpolation)", "", bnk->mu2, &bnk->mu2, NULL)); 1149 PetscCall(PetscOptionsReal("-tao_bnk_gamma1", "(developer) radius reduction factor for rejected very bad step (-tao_bnk_update_type interpolation)", "", bnk->gamma1, &bnk->gamma1, NULL)); 1150 PetscCall(PetscOptionsReal("-tao_bnk_gamma2", "(developer) radius reduction factor for rejected bad step (-tao_bnk_update_type interpolation)", "", bnk->gamma2, &bnk->gamma2, NULL)); 1151 PetscCall(PetscOptionsReal("-tao_bnk_gamma3", "(developer) radius increase factor for accepted good step (-tao_bnk_update_type interpolation)", "", bnk->gamma3, &bnk->gamma3, NULL)); 1152 PetscCall(PetscOptionsReal("-tao_bnk_gamma4", "(developer) radius increase factor for accepted very good step (-tao_bnk_update_type interpolation)", "", bnk->gamma4, &bnk->gamma4, NULL)); 1153 PetscCall(PetscOptionsReal("-tao_bnk_theta", "(developer) trust region interpolation factor (-tao_bnk_update_type interpolation)", "", bnk->theta, &bnk->theta, NULL)); 1154 PetscCall(PetscOptionsReal("-tao_bnk_min_radius", "(developer) lower bound on initial radius", "", bnk->min_radius, &bnk->min_radius, NULL)); 1155 PetscCall(PetscOptionsReal("-tao_bnk_max_radius", "(developer) upper bound on radius", "", bnk->max_radius, &bnk->max_radius, NULL)); 1156 PetscCall(PetscOptionsReal("-tao_bnk_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon, NULL)); 1157 PetscCall(PetscOptionsReal("-tao_bnk_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol, NULL)); 1158 PetscCall(PetscOptionsReal("-tao_bnk_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step, NULL)); 1159 PetscCall(PetscOptionsInt("-tao_bnk_max_cg_its", "number of BNCG iterations to take for each Newton step", "", bnk->max_cg_its, &bnk->max_cg_its, NULL)); 1160 PetscOptionsHeadEnd(); 1161 1162 PetscCall(TaoSetOptionsPrefix(bnk->bncg, ((PetscObject)(tao))->prefix)); 1163 PetscCall(TaoAppendOptionsPrefix(bnk->bncg, "tao_bnk_cg_")); 1164 PetscCall(TaoSetFromOptions(bnk->bncg)); 1165 1166 PetscCall(KSPSetOptionsPrefix(tao->ksp, ((PetscObject)(tao))->prefix)); 1167 PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_bnk_")); 1168 PetscCall(KSPSetFromOptions(tao->ksp)); 1169 PetscFunctionReturn(0); 1170 } 1171 1172 /*------------------------------------------------------------*/ 1173 1174 PetscErrorCode TaoView_BNK(Tao tao, PetscViewer viewer) 1175 { 1176 TAO_BNK *bnk = (TAO_BNK *)tao->data; 1177 PetscInt nrejects; 1178 PetscBool isascii; 1179 1180 PetscFunctionBegin; 1181 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 1182 if (isascii) { 1183 PetscCall(PetscViewerASCIIPushTab(viewer)); 1184 PetscCall(TaoView(bnk->bncg, viewer)); 1185 if (bnk->M) { 1186 PetscCall(MatLMVMGetRejectCount(bnk->M, &nrejects)); 1187 PetscCall(PetscViewerASCIIPrintf(viewer, "Rejected BFGS updates: %" PetscInt_FMT "\n", nrejects)); 1188 } 1189 PetscCall(PetscViewerASCIIPrintf(viewer, "CG steps: %" PetscInt_FMT "\n", bnk->tot_cg_its)); 1190 PetscCall(PetscViewerASCIIPrintf(viewer, "Newton steps: %" PetscInt_FMT "\n", bnk->newt)); 1191 if (bnk->M) PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS steps: %" PetscInt_FMT "\n", bnk->bfgs)); 1192 PetscCall(PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %" PetscInt_FMT "\n", bnk->sgrad)); 1193 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", bnk->grad)); 1194 PetscCall(PetscViewerASCIIPrintf(viewer, "KSP termination reasons:\n")); 1195 PetscCall(PetscViewerASCIIPrintf(viewer, " atol: %" PetscInt_FMT "\n", bnk->ksp_atol)); 1196 PetscCall(PetscViewerASCIIPrintf(viewer, " rtol: %" PetscInt_FMT "\n", bnk->ksp_rtol)); 1197 PetscCall(PetscViewerASCIIPrintf(viewer, " ctol: %" PetscInt_FMT "\n", bnk->ksp_ctol)); 1198 PetscCall(PetscViewerASCIIPrintf(viewer, " negc: %" PetscInt_FMT "\n", bnk->ksp_negc)); 1199 PetscCall(PetscViewerASCIIPrintf(viewer, " dtol: %" PetscInt_FMT "\n", bnk->ksp_dtol)); 1200 PetscCall(PetscViewerASCIIPrintf(viewer, " iter: %" PetscInt_FMT "\n", bnk->ksp_iter)); 1201 PetscCall(PetscViewerASCIIPrintf(viewer, " othr: %" PetscInt_FMT "\n", bnk->ksp_othr)); 1202 PetscCall(PetscViewerASCIIPopTab(viewer)); 1203 } 1204 PetscFunctionReturn(0); 1205 } 1206 1207 /* ---------------------------------------------------------- */ 1208 1209 /*MC 1210 TAOBNK - Shared base-type for Bounded Newton-Krylov type algorithms. 1211 At each iteration, the BNK methods solve the symmetric 1212 system of equations to obtain the step diretion dk: 1213 Hk dk = -gk 1214 for free variables only. The step can be globalized either through 1215 trust-region methods, or a line search, or a heuristic mixture of both. 1216 1217 Options Database Keys: 1218 + -tao_bnk_max_cg_its - maximum number of bounded conjugate-gradient iterations taken in each Newton loop 1219 . -tao_bnk_init_type - trust radius initialization method ("constant", "direction", "interpolation") 1220 . -tao_bnk_update_type - trust radius update method ("step", "direction", "interpolation") 1221 . -tao_bnk_as_type - active-set estimation method ("none", "bertsekas") 1222 . -tao_bnk_as_tol - (developer) initial tolerance used in estimating bounded active variables (-as_type bertsekas) 1223 . -tao_bnk_as_step - (developer) trial step length used in estimating bounded active variables (-as_type bertsekas) 1224 . -tao_bnk_sval - (developer) Hessian perturbation starting value 1225 . -tao_bnk_imin - (developer) minimum initial Hessian perturbation 1226 . -tao_bnk_imax - (developer) maximum initial Hessian perturbation 1227 . -tao_bnk_pmin - (developer) minimum Hessian perturbation 1228 . -tao_bnk_pmax - (developer) aximum Hessian perturbation 1229 . -tao_bnk_pgfac - (developer) Hessian perturbation growth factor 1230 . -tao_bnk_psfac - (developer) Hessian perturbation shrink factor 1231 . -tao_bnk_imfac - (developer) initial merit factor for Hessian perturbation 1232 . -tao_bnk_pmgfac - (developer) merit growth factor for Hessian perturbation 1233 . -tao_bnk_pmsfac - (developer) merit shrink factor for Hessian perturbation 1234 . -tao_bnk_eta1 - (developer) threshold for rejecting step (-update_type reduction) 1235 . -tao_bnk_eta2 - (developer) threshold for accepting marginal step (-update_type reduction) 1236 . -tao_bnk_eta3 - (developer) threshold for accepting reasonable step (-update_type reduction) 1237 . -tao_bnk_eta4 - (developer) threshold for accepting good step (-update_type reduction) 1238 . -tao_bnk_alpha1 - (developer) radius reduction factor for rejected step (-update_type reduction) 1239 . -tao_bnk_alpha2 - (developer) radius reduction factor for marginally accepted bad step (-update_type reduction) 1240 . -tao_bnk_alpha3 - (developer) radius increase factor for reasonable accepted step (-update_type reduction) 1241 . -tao_bnk_alpha4 - (developer) radius increase factor for good accepted step (-update_type reduction) 1242 . -tao_bnk_alpha5 - (developer) radius increase factor for very good accepted step (-update_type reduction) 1243 . -tao_bnk_epsilon - (developer) tolerance for small pred/actual ratios that trigger automatic step acceptance (-update_type reduction) 1244 . -tao_bnk_mu1 - (developer) threshold for accepting very good step (-update_type interpolation) 1245 . -tao_bnk_mu2 - (developer) threshold for accepting good step (-update_type interpolation) 1246 . -tao_bnk_gamma1 - (developer) radius reduction factor for rejected very bad step (-update_type interpolation) 1247 . -tao_bnk_gamma2 - (developer) radius reduction factor for rejected bad step (-update_type interpolation) 1248 . -tao_bnk_gamma3 - (developer) radius increase factor for accepted good step (-update_type interpolation) 1249 . -tao_bnk_gamma4 - (developer) radius increase factor for accepted very good step (-update_type interpolation) 1250 . -tao_bnk_theta - (developer) trust region interpolation factor (-update_type interpolation) 1251 . -tao_bnk_nu1 - (developer) threshold for small line-search step length (-update_type step) 1252 . -tao_bnk_nu2 - (developer) threshold for reasonable line-search step length (-update_type step) 1253 . -tao_bnk_nu3 - (developer) threshold for large line-search step length (-update_type step) 1254 . -tao_bnk_nu4 - (developer) threshold for very large line-search step length (-update_type step) 1255 . -tao_bnk_omega1 - (developer) radius reduction factor for very small line-search step length (-update_type step) 1256 . -tao_bnk_omega2 - (developer) radius reduction factor for small line-search step length (-update_type step) 1257 . -tao_bnk_omega3 - (developer) radius factor for decent line-search step length (-update_type step) 1258 . -tao_bnk_omega4 - (developer) radius increase factor for large line-search step length (-update_type step) 1259 . -tao_bnk_omega5 - (developer) radius increase factor for very large line-search step length (-update_type step) 1260 . -tao_bnk_mu1_i - (developer) threshold for accepting very good step (-init_type interpolation) 1261 . -tao_bnk_mu2_i - (developer) threshold for accepting good step (-init_type interpolation) 1262 . -tao_bnk_gamma1_i - (developer) radius reduction factor for rejected very bad step (-init_type interpolation) 1263 . -tao_bnk_gamma2_i - (developer) radius reduction factor for rejected bad step (-init_type interpolation) 1264 . -tao_bnk_gamma3_i - (developer) radius increase factor for accepted good step (-init_type interpolation) 1265 . -tao_bnk_gamma4_i - (developer) radius increase factor for accepted very good step (-init_type interpolation) 1266 - -tao_bnk_theta_i - (developer) trust region interpolation factor (-init_type interpolation) 1267 1268 Level: beginner 1269 M*/ 1270 1271 PetscErrorCode TaoCreate_BNK(Tao tao) 1272 { 1273 TAO_BNK *bnk; 1274 PC pc; 1275 1276 PetscFunctionBegin; 1277 PetscCall(PetscNew(&bnk)); 1278 1279 tao->ops->setup = TaoSetUp_BNK; 1280 tao->ops->view = TaoView_BNK; 1281 tao->ops->setfromoptions = TaoSetFromOptions_BNK; 1282 tao->ops->destroy = TaoDestroy_BNK; 1283 1284 /* Override default settings (unless already changed) */ 1285 if (!tao->max_it_changed) tao->max_it = 50; 1286 if (!tao->trust0_changed) tao->trust0 = 100.0; 1287 1288 tao->data = (void *)bnk; 1289 1290 /* Hessian shifting parameters */ 1291 bnk->computehessian = TaoBNKComputeHessian; 1292 bnk->computestep = TaoBNKComputeStep; 1293 1294 bnk->sval = 0.0; 1295 bnk->imin = 1.0e-4; 1296 bnk->imax = 1.0e+2; 1297 bnk->imfac = 1.0e-1; 1298 1299 bnk->pmin = 1.0e-12; 1300 bnk->pmax = 1.0e+2; 1301 bnk->pgfac = 1.0e+1; 1302 bnk->psfac = 4.0e-1; 1303 bnk->pmgfac = 1.0e-1; 1304 bnk->pmsfac = 1.0e-1; 1305 1306 /* Default values for trust-region radius update based on steplength */ 1307 bnk->nu1 = 0.25; 1308 bnk->nu2 = 0.50; 1309 bnk->nu3 = 1.00; 1310 bnk->nu4 = 1.25; 1311 1312 bnk->omega1 = 0.25; 1313 bnk->omega2 = 0.50; 1314 bnk->omega3 = 1.00; 1315 bnk->omega4 = 2.00; 1316 bnk->omega5 = 4.00; 1317 1318 /* Default values for trust-region radius update based on reduction */ 1319 bnk->eta1 = 1.0e-4; 1320 bnk->eta2 = 0.25; 1321 bnk->eta3 = 0.50; 1322 bnk->eta4 = 0.90; 1323 1324 bnk->alpha1 = 0.25; 1325 bnk->alpha2 = 0.50; 1326 bnk->alpha3 = 1.00; 1327 bnk->alpha4 = 2.00; 1328 bnk->alpha5 = 4.00; 1329 1330 /* Default values for trust-region radius update based on interpolation */ 1331 bnk->mu1 = 0.10; 1332 bnk->mu2 = 0.50; 1333 1334 bnk->gamma1 = 0.25; 1335 bnk->gamma2 = 0.50; 1336 bnk->gamma3 = 2.00; 1337 bnk->gamma4 = 4.00; 1338 1339 bnk->theta = 0.05; 1340 1341 /* Default values for trust region initialization based on interpolation */ 1342 bnk->mu1_i = 0.35; 1343 bnk->mu2_i = 0.50; 1344 1345 bnk->gamma1_i = 0.0625; 1346 bnk->gamma2_i = 0.5; 1347 bnk->gamma3_i = 2.0; 1348 bnk->gamma4_i = 5.0; 1349 1350 bnk->theta_i = 0.25; 1351 1352 /* Remaining parameters */ 1353 bnk->max_cg_its = 0; 1354 bnk->min_radius = 1.0e-10; 1355 bnk->max_radius = 1.0e10; 1356 bnk->epsilon = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0); 1357 bnk->as_tol = 1.0e-3; 1358 bnk->as_step = 1.0e-3; 1359 bnk->dmin = 1.0e-6; 1360 bnk->dmax = 1.0e6; 1361 1362 bnk->M = NULL; 1363 bnk->bfgs_pre = NULL; 1364 bnk->init_type = BNK_INIT_INTERPOLATION; 1365 bnk->update_type = BNK_UPDATE_REDUCTION; 1366 bnk->as_type = BNK_AS_BERTSEKAS; 1367 1368 /* Create the embedded BNCG solver */ 1369 PetscCall(TaoCreate(PetscObjectComm((PetscObject)tao), &bnk->bncg)); 1370 PetscCall(PetscObjectIncrementTabLevel((PetscObject)bnk->bncg, (PetscObject)tao, 1)); 1371 PetscCall(TaoSetType(bnk->bncg, TAOBNCG)); 1372 1373 /* Create the line search */ 1374 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 1375 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 1376 PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT)); 1377 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao)); 1378 1379 /* Set linear solver to default for symmetric matrices */ 1380 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 1381 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 1382 PetscCall(KSPSetType(tao->ksp, KSPSTCG)); 1383 PetscCall(KSPGetPC(tao->ksp, &pc)); 1384 PetscCall(PCSetType(pc, PCLMVM)); 1385 PetscFunctionReturn(0); 1386 } 1387