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 PetscCall((*tao->ops->convergencetest)(tao,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 PetscCall((*tao->ops->convergencetest)(tao,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) { 324 PetscCall(MatAssembled(tao->hessian, &hessComputed)); 325 } 326 if (hessComputed) { 327 PetscCall(MatHasOperation(tao->hessian, MATOP_GET_DIAGONAL, &diagExists)); 328 } 329 if (diagExists) { 330 /* BFGS preconditioner doesn't exist so let's invert the absolute diagonal of the Hessian instead onto the gradient */ 331 PetscCall(MatGetDiagonal(tao->hessian, bnk->Xwork)); 332 PetscCall(VecAbs(bnk->Xwork)); 333 PetscCall(VecMedian(bnk->Diag_min, bnk->Xwork, bnk->Diag_max, bnk->Xwork)); 334 PetscCall(VecReciprocal(bnk->Xwork)); 335 PetscCall(VecPointwiseMult(bnk->W, bnk->Xwork, bnk->unprojected_gradient)); 336 } else { 337 /* If the Hessian or its diagonal does not exist, we will simply use gradient step */ 338 PetscCall(VecCopy(bnk->unprojected_gradient, bnk->W)); 339 } 340 } 341 PetscCall(VecScale(bnk->W, -1.0)); 342 PetscCall(TaoEstimateActiveBounds(tao->solution, tao->XL, tao->XU, bnk->unprojected_gradient, bnk->W, bnk->Xwork, bnk->as_step, &bnk->as_tol, 343 &bnk->active_lower, &bnk->active_upper, &bnk->active_fixed, &bnk->active_idx, &bnk->inactive_idx)); 344 break; 345 346 default: 347 break; 348 } 349 bnk->resetksp = (PetscBool)(bnk->active_idx || hadactive); /* inactive Hessian size may have changed, need to reset operators */ 350 PetscFunctionReturn(0); 351 } 352 353 /*------------------------------------------------------------*/ 354 355 /* Routine for bounding the step direction */ 356 357 PetscErrorCode TaoBNKBoundStep(Tao tao, PetscInt asType, Vec step) 358 { 359 TAO_BNK *bnk = (TAO_BNK *)tao->data; 360 361 PetscFunctionBegin; 362 switch (asType) { 363 case BNK_AS_NONE: 364 PetscCall(VecISSet(step, bnk->active_idx, 0.0)); 365 break; 366 367 case BNK_AS_BERTSEKAS: 368 PetscCall(TaoBoundStep(tao->solution, tao->XL, tao->XU, bnk->active_lower, bnk->active_upper, bnk->active_fixed, 1.0, step)); 369 break; 370 371 default: 372 break; 373 } 374 PetscFunctionReturn(0); 375 } 376 377 /*------------------------------------------------------------*/ 378 379 /* Routine for taking a finite number of BNCG iterations to 380 accelerate Newton convergence. 381 382 In practice, this approach simply trades off Hessian evaluations 383 for more gradient evaluations. 384 */ 385 386 PetscErrorCode TaoBNKTakeCGSteps(Tao tao, PetscBool *terminate) 387 { 388 TAO_BNK *bnk = (TAO_BNK *)tao->data; 389 390 PetscFunctionBegin; 391 *terminate = PETSC_FALSE; 392 if (bnk->max_cg_its > 0) { 393 /* Copy the current function value (important vectors are already shared) */ 394 bnk->bncg_ctx->f = bnk->f; 395 /* Take some small finite number of BNCG iterations */ 396 PetscCall(TaoSolve(bnk->bncg)); 397 /* Add the number of gradient and function evaluations to the total */ 398 tao->nfuncs += bnk->bncg->nfuncs; 399 tao->nfuncgrads += bnk->bncg->nfuncgrads; 400 tao->ngrads += bnk->bncg->ngrads; 401 tao->nhess += bnk->bncg->nhess; 402 bnk->tot_cg_its += bnk->bncg->niter; 403 /* Extract the BNCG function value out and save it into BNK */ 404 bnk->f = bnk->bncg_ctx->f; 405 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) { 406 *terminate = PETSC_TRUE; 407 } else { 408 PetscCall(TaoBNKEstimateActiveSet(tao, bnk->as_type)); 409 } 410 } 411 PetscFunctionReturn(0); 412 } 413 414 /*------------------------------------------------------------*/ 415 416 /* Routine for computing the Newton step. */ 417 418 PetscErrorCode TaoBNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type) 419 { 420 TAO_BNK *bnk = (TAO_BNK *)tao->data; 421 PetscInt bfgsUpdates = 0; 422 PetscInt kspits; 423 PetscBool is_lmvm; 424 PetscVoidFunction kspTR; 425 426 PetscFunctionBegin; 427 /* If there are no inactive variables left, save some computation and return an adjusted zero step 428 that has (l-x) and (u-x) for lower and upper bounded variables. */ 429 if (!bnk->inactive_idx) { 430 PetscCall(VecSet(tao->stepdirection, 0.0)); 431 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 432 PetscFunctionReturn(0); 433 } 434 435 /* Shift the reduced Hessian matrix */ 436 if (shift && bnk->pert > 0) { 437 PetscCall(PetscObjectTypeCompare((PetscObject)tao->hessian, MATLMVM, &is_lmvm)); 438 if (is_lmvm) { 439 PetscCall(MatShift(tao->hessian, bnk->pert)); 440 } else { 441 PetscCall(MatShift(bnk->H_inactive, bnk->pert)); 442 if (bnk->H_inactive != bnk->Hpre_inactive) { 443 PetscCall(MatShift(bnk->Hpre_inactive, bnk->pert)); 444 } 445 } 446 } 447 448 /* Solve the Newton system of equations */ 449 tao->ksp_its = 0; 450 PetscCall(VecSet(tao->stepdirection, 0.0)); 451 if (bnk->resetksp) { 452 PetscCall(KSPReset(tao->ksp)); 453 PetscCall(KSPResetFromOptions(tao->ksp)); 454 bnk->resetksp = PETSC_FALSE; 455 } 456 PetscCall(KSPSetOperators(tao->ksp,bnk->H_inactive,bnk->Hpre_inactive)); 457 PetscCall(VecCopy(bnk->unprojected_gradient, bnk->Gwork)); 458 if (bnk->active_idx) { 459 PetscCall(VecGetSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive)); 460 PetscCall(VecGetSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive)); 461 } else { 462 bnk->G_inactive = bnk->unprojected_gradient; 463 bnk->X_inactive = tao->stepdirection; 464 } 465 PetscCall(KSPCGSetRadius(tao->ksp,tao->trust)); 466 PetscCall(KSPSolve(tao->ksp, bnk->G_inactive, bnk->X_inactive)); 467 PetscCall(KSPGetIterationNumber(tao->ksp,&kspits)); 468 tao->ksp_its += kspits; 469 tao->ksp_tot_its += kspits; 470 PetscCall(PetscObjectQueryFunction((PetscObject)tao->ksp,"KSPCGGetNormD_C",&kspTR)); 471 if (kspTR) { 472 PetscCall(KSPCGGetNormD(tao->ksp,&bnk->dnorm)); 473 474 if (0.0 == tao->trust) { 475 /* Radius was uninitialized; use the norm of the direction */ 476 if (bnk->dnorm > 0.0) { 477 tao->trust = bnk->dnorm; 478 479 /* Modify the radius if it is too large or small */ 480 tao->trust = PetscMax(tao->trust, bnk->min_radius); 481 tao->trust = PetscMin(tao->trust, bnk->max_radius); 482 } else { 483 /* The direction was bad; set radius to default value and re-solve 484 the trust-region subproblem to get a direction */ 485 tao->trust = tao->trust0; 486 487 /* Modify the radius if it is too large or small */ 488 tao->trust = PetscMax(tao->trust, bnk->min_radius); 489 tao->trust = PetscMin(tao->trust, bnk->max_radius); 490 491 PetscCall(KSPCGSetRadius(tao->ksp,tao->trust)); 492 PetscCall(KSPSolve(tao->ksp, bnk->G_inactive, bnk->X_inactive)); 493 PetscCall(KSPGetIterationNumber(tao->ksp,&kspits)); 494 tao->ksp_its += kspits; 495 tao->ksp_tot_its += kspits; 496 PetscCall(KSPCGGetNormD(tao->ksp,&bnk->dnorm)); 497 498 PetscCheck(bnk->dnorm != 0.0,PetscObjectComm((PetscObject)tao),PETSC_ERR_PLIB, "Initial direction zero"); 499 } 500 } 501 } 502 /* Restore sub vectors back */ 503 if (bnk->active_idx) { 504 PetscCall(VecRestoreSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive)); 505 PetscCall(VecRestoreSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive)); 506 } 507 /* Make sure the safeguarded fall-back step is zero for actively bounded variables */ 508 PetscCall(VecScale(tao->stepdirection, -1.0)); 509 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 510 511 /* Record convergence reasons */ 512 PetscCall(KSPGetConvergedReason(tao->ksp, ksp_reason)); 513 if (KSP_CONVERGED_ATOL == *ksp_reason) { 514 ++bnk->ksp_atol; 515 } else if (KSP_CONVERGED_RTOL == *ksp_reason) { 516 ++bnk->ksp_rtol; 517 } else if (KSP_CONVERGED_CG_CONSTRAINED == *ksp_reason) { 518 ++bnk->ksp_ctol; 519 } else if (KSP_CONVERGED_CG_NEG_CURVE == *ksp_reason) { 520 ++bnk->ksp_negc; 521 } else if (KSP_DIVERGED_DTOL == *ksp_reason) { 522 ++bnk->ksp_dtol; 523 } else if (KSP_DIVERGED_ITS == *ksp_reason) { 524 ++bnk->ksp_iter; 525 } else { 526 ++bnk->ksp_othr; 527 } 528 529 /* Make sure the BFGS preconditioner is healthy */ 530 if (bnk->M) { 531 PetscCall(MatLMVMGetUpdateCount(bnk->M, &bfgsUpdates)); 532 if ((KSP_DIVERGED_INDEFINITE_PC == *ksp_reason) && (bfgsUpdates > 0)) { 533 /* Preconditioner is numerically indefinite; reset the approximation. */ 534 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 535 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 536 } 537 } 538 *step_type = BNK_NEWTON; 539 PetscFunctionReturn(0); 540 } 541 542 /*------------------------------------------------------------*/ 543 544 /* Routine for recomputing the predicted reduction for a given step vector */ 545 546 PetscErrorCode TaoBNKRecomputePred(Tao tao, Vec S, PetscReal *prered) 547 { 548 TAO_BNK *bnk = (TAO_BNK *)tao->data; 549 550 PetscFunctionBegin; 551 /* Extract subvectors associated with the inactive set */ 552 if (bnk->active_idx) { 553 PetscCall(VecGetSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive)); 554 PetscCall(VecGetSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work)); 555 PetscCall(VecGetSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive)); 556 } else { 557 bnk->X_inactive = tao->stepdirection; 558 bnk->inactive_work = bnk->Xwork; 559 bnk->G_inactive = bnk->Gwork; 560 } 561 /* Recompute the predicted decrease based on the quadratic model */ 562 PetscCall(MatMult(bnk->H_inactive, bnk->X_inactive, bnk->inactive_work)); 563 PetscCall(VecAYPX(bnk->inactive_work, -0.5, bnk->G_inactive)); 564 PetscCall(VecDot(bnk->inactive_work, bnk->X_inactive, prered)); 565 /* Restore the sub vectors */ 566 if (bnk->active_idx) { 567 PetscCall(VecRestoreSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive)); 568 PetscCall(VecRestoreSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work)); 569 PetscCall(VecRestoreSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive)); 570 } 571 PetscFunctionReturn(0); 572 } 573 574 /*------------------------------------------------------------*/ 575 576 /* Routine for ensuring that the Newton step is a descent direction. 577 578 The step direction falls back onto BFGS, scaled gradient and gradient steps 579 in the event that the Newton step fails the test. 580 */ 581 582 PetscErrorCode TaoBNKSafeguardStep(Tao tao, KSPConvergedReason ksp_reason, PetscInt *stepType) 583 { 584 TAO_BNK *bnk = (TAO_BNK *)tao->data; 585 PetscReal gdx, e_min; 586 PetscInt bfgsUpdates; 587 588 PetscFunctionBegin; 589 switch (*stepType) { 590 case BNK_NEWTON: 591 PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx)); 592 if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) { 593 /* Newton step is not descent or direction produced Inf or NaN 594 Update the perturbation for next time */ 595 if (bnk->pert <= 0.0) { 596 PetscBool is_gltr; 597 598 /* Initialize the perturbation */ 599 bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm)); 600 PetscCall(PetscObjectTypeCompare((PetscObject)(tao->ksp),KSPGLTR,&is_gltr)); 601 if (is_gltr) { 602 PetscCall(KSPGLTRGetMinEig(tao->ksp,&e_min)); 603 bnk->pert = PetscMax(bnk->pert, -e_min); 604 } 605 } else { 606 /* Increase the perturbation */ 607 bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm)); 608 } 609 610 if (!bnk->M) { 611 /* We don't have the bfgs matrix around and updated 612 Must use gradient direction in this case */ 613 PetscCall(VecCopy(tao->gradient, tao->stepdirection)); 614 *stepType = BNK_GRADIENT; 615 } else { 616 /* Attempt to use the BFGS direction */ 617 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 618 619 /* Check for success (descent direction) 620 NOTE: Negative gdx here means not a descent direction because 621 the fall-back step is missing a negative sign. */ 622 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx)); 623 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 624 /* BFGS direction is not descent or direction produced not a number 625 We can assert bfgsUpdates > 1 in this case because 626 the first solve produces the scaled gradient direction, 627 which is guaranteed to be descent */ 628 629 /* Use steepest descent direction (scaled) */ 630 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 631 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 632 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 633 634 *stepType = BNK_SCALED_GRADIENT; 635 } else { 636 PetscCall(MatLMVMGetUpdateCount(bnk->M, &bfgsUpdates)); 637 if (1 == bfgsUpdates) { 638 /* The first BFGS direction is always the scaled gradient */ 639 *stepType = BNK_SCALED_GRADIENT; 640 } else { 641 *stepType = BNK_BFGS; 642 } 643 } 644 } 645 /* Make sure the safeguarded fall-back step is zero for actively bounded variables */ 646 PetscCall(VecScale(tao->stepdirection, -1.0)); 647 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 648 } else { 649 /* Computed Newton step is descent */ 650 switch (ksp_reason) { 651 case KSP_DIVERGED_NANORINF: 652 case KSP_DIVERGED_BREAKDOWN: 653 case KSP_DIVERGED_INDEFINITE_MAT: 654 case KSP_DIVERGED_INDEFINITE_PC: 655 case KSP_CONVERGED_CG_NEG_CURVE: 656 /* Matrix or preconditioner is indefinite; increase perturbation */ 657 if (bnk->pert <= 0.0) { 658 PetscBool is_gltr; 659 660 /* Initialize the perturbation */ 661 bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm)); 662 PetscCall(PetscObjectTypeCompare((PetscObject)(tao->ksp),KSPGLTR,&is_gltr)); 663 if (is_gltr) { 664 PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min)); 665 bnk->pert = PetscMax(bnk->pert, -e_min); 666 } 667 } else { 668 /* Increase the perturbation */ 669 bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm)); 670 } 671 break; 672 673 default: 674 /* Newton step computation is good; decrease perturbation */ 675 bnk->pert = PetscMin(bnk->psfac * bnk->pert, bnk->pmsfac * bnk->gnorm); 676 if (bnk->pert < bnk->pmin) { 677 bnk->pert = 0.0; 678 } 679 break; 680 } 681 *stepType = BNK_NEWTON; 682 } 683 break; 684 685 case BNK_BFGS: 686 /* Check for success (descent direction) */ 687 PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx)); 688 if (gdx >= 0 || PetscIsInfOrNanReal(gdx)) { 689 /* Step is not descent or solve was not successful 690 Use steepest descent direction (scaled) */ 691 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 692 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 693 PetscCall(MatSolve(bnk->M, tao->gradient, tao->stepdirection)); 694 PetscCall(VecScale(tao->stepdirection,-1.0)); 695 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 696 *stepType = BNK_SCALED_GRADIENT; 697 } else { 698 *stepType = BNK_BFGS; 699 } 700 break; 701 702 case BNK_SCALED_GRADIENT: 703 break; 704 705 default: 706 break; 707 } 708 709 PetscFunctionReturn(0); 710 } 711 712 /*------------------------------------------------------------*/ 713 714 /* Routine for performing a bound-projected More-Thuente line search. 715 716 Includes fallbacks to BFGS, scaled gradient, and unscaled gradient steps if the 717 Newton step does not produce a valid step length. 718 */ 719 720 PetscErrorCode TaoBNKPerformLineSearch(Tao tao, PetscInt *stepType, PetscReal *steplen, TaoLineSearchConvergedReason *reason) 721 { 722 TAO_BNK *bnk = (TAO_BNK *)tao->data; 723 TaoLineSearchConvergedReason ls_reason; 724 PetscReal e_min, gdx; 725 PetscInt bfgsUpdates; 726 727 PetscFunctionBegin; 728 /* Perform the linesearch */ 729 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, bnk->unprojected_gradient, tao->stepdirection, steplen, &ls_reason)); 730 PetscCall(TaoAddLineSearchCounts(tao)); 731 732 while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && *stepType != BNK_SCALED_GRADIENT && *stepType != BNK_GRADIENT) { 733 /* Linesearch failed, revert solution */ 734 bnk->f = bnk->fold; 735 PetscCall(VecCopy(bnk->Xold, tao->solution)); 736 PetscCall(VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient)); 737 738 switch(*stepType) { 739 case BNK_NEWTON: 740 /* Failed to obtain acceptable iterate with Newton step 741 Update the perturbation for next time */ 742 if (bnk->pert <= 0.0) { 743 PetscBool is_gltr; 744 745 /* Initialize the perturbation */ 746 bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm)); 747 PetscCall(PetscObjectTypeCompare((PetscObject)(tao->ksp),KSPGLTR,&is_gltr)); 748 if (is_gltr) { 749 PetscCall(KSPGLTRGetMinEig(tao->ksp,&e_min)); 750 bnk->pert = PetscMax(bnk->pert, -e_min); 751 } 752 } else { 753 /* Increase the perturbation */ 754 bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm)); 755 } 756 757 if (!bnk->M) { 758 /* We don't have the bfgs matrix around and being updated 759 Must use gradient direction in this case */ 760 PetscCall(VecCopy(bnk->unprojected_gradient, tao->stepdirection)); 761 *stepType = BNK_GRADIENT; 762 } else { 763 /* Attempt to use the BFGS direction */ 764 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 765 /* Check for success (descent direction) 766 NOTE: Negative gdx means not a descent direction because the step here is missing a negative sign. */ 767 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx)); 768 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 769 /* BFGS direction is not descent or direction produced not a number 770 We can assert bfgsUpdates > 1 in this case 771 Use steepest descent direction (scaled) */ 772 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 773 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 774 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 775 776 bfgsUpdates = 1; 777 *stepType = BNK_SCALED_GRADIENT; 778 } else { 779 PetscCall(MatLMVMGetUpdateCount(bnk->M, &bfgsUpdates)); 780 if (1 == bfgsUpdates) { 781 /* The first BFGS direction is always the scaled gradient */ 782 *stepType = BNK_SCALED_GRADIENT; 783 } else { 784 *stepType = BNK_BFGS; 785 } 786 } 787 } 788 break; 789 790 case BNK_BFGS: 791 /* Can only enter if pc_type == BNK_PC_BFGS 792 Failed to obtain acceptable iterate with BFGS step 793 Attempt to use the scaled gradient direction */ 794 PetscCall(MatLMVMReset(bnk->M, PETSC_FALSE)); 795 PetscCall(MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient)); 796 PetscCall(MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection)); 797 798 bfgsUpdates = 1; 799 *stepType = BNK_SCALED_GRADIENT; 800 break; 801 } 802 /* Make sure the safeguarded fall-back step is zero for actively bounded variables */ 803 PetscCall(VecScale(tao->stepdirection, -1.0)); 804 PetscCall(TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection)); 805 806 /* Perform one last line search with the fall-back step */ 807 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, bnk->unprojected_gradient, tao->stepdirection, steplen, &ls_reason)); 808 PetscCall(TaoAddLineSearchCounts(tao)); 809 } 810 *reason = ls_reason; 811 PetscFunctionReturn(0); 812 } 813 814 /*------------------------------------------------------------*/ 815 816 /* Routine for updating the trust radius. 817 818 Function features three different update methods: 819 1) Line-search step length based 820 2) Predicted decrease on the CG quadratic model 821 3) Interpolation 822 */ 823 824 PetscErrorCode TaoBNKUpdateTrustRadius(Tao tao, PetscReal prered, PetscReal actred, PetscInt updateType, PetscInt stepType, PetscBool *accept) 825 { 826 TAO_BNK *bnk = (TAO_BNK *)tao->data; 827 828 PetscReal step, kappa; 829 PetscReal gdx, tau_1, tau_2, tau_min, tau_max; 830 831 PetscFunctionBegin; 832 /* Update trust region radius */ 833 *accept = PETSC_FALSE; 834 switch(updateType) { 835 case BNK_UPDATE_STEP: 836 *accept = PETSC_TRUE; /* always accept here because line search succeeded */ 837 if (stepType == BNK_NEWTON) { 838 PetscCall(TaoLineSearchGetStepLength(tao->linesearch, &step)); 839 if (step < bnk->nu1) { 840 /* Very bad step taken; reduce radius */ 841 tao->trust = bnk->omega1 * PetscMin(bnk->dnorm, tao->trust); 842 } else if (step < bnk->nu2) { 843 /* Reasonably bad step taken; reduce radius */ 844 tao->trust = bnk->omega2 * PetscMin(bnk->dnorm, tao->trust); 845 } else if (step < bnk->nu3) { 846 /* Reasonable step was taken; leave radius alone */ 847 if (bnk->omega3 < 1.0) { 848 tao->trust = bnk->omega3 * PetscMin(bnk->dnorm, tao->trust); 849 } else if (bnk->omega3 > 1.0) { 850 tao->trust = PetscMax(bnk->omega3 * bnk->dnorm, tao->trust); 851 } 852 } else if (step < bnk->nu4) { 853 /* Full step taken; increase the radius */ 854 tao->trust = PetscMax(bnk->omega4 * bnk->dnorm, tao->trust); 855 } else { 856 /* More than full step taken; increase the radius */ 857 tao->trust = PetscMax(bnk->omega5 * bnk->dnorm, tao->trust); 858 } 859 } else { 860 /* Newton step was not good; reduce the radius */ 861 tao->trust = bnk->omega1 * PetscMin(bnk->dnorm, tao->trust); 862 } 863 break; 864 865 case BNK_UPDATE_REDUCTION: 866 if (stepType == BNK_NEWTON) { 867 if ((prered < 0.0) || PetscIsInfOrNanReal(prered)) { 868 /* The predicted reduction has the wrong sign. This cannot 869 happen in infinite precision arithmetic. Step should 870 be rejected! */ 871 tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm); 872 } else { 873 if (PetscIsInfOrNanReal(actred)) { 874 tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm); 875 } else { 876 if ((PetscAbsScalar(actred) <= PetscMax(1.0, PetscAbsScalar(bnk->f))*bnk->epsilon) && (PetscAbsScalar(prered) <= PetscMax(1.0, PetscAbsScalar(bnk->f))*bnk->epsilon)) { 877 kappa = 1.0; 878 } else { 879 kappa = actred / prered; 880 } 881 /* Accept or reject the step and update radius */ 882 if (kappa < bnk->eta1) { 883 /* Reject the step */ 884 tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm); 885 } else { 886 /* Accept the step */ 887 *accept = PETSC_TRUE; 888 /* Update the trust region radius only if the computed step is at the trust radius boundary */ 889 if (bnk->dnorm == tao->trust) { 890 if (kappa < bnk->eta2) { 891 /* Marginal bad step */ 892 tao->trust = bnk->alpha2 * tao->trust; 893 } else if (kappa < bnk->eta3) { 894 /* Reasonable step */ 895 tao->trust = bnk->alpha3 * tao->trust; 896 } else if (kappa < bnk->eta4) { 897 /* Good step */ 898 tao->trust = bnk->alpha4 * tao->trust; 899 } else { 900 /* Very good step */ 901 tao->trust = bnk->alpha5 * tao->trust; 902 } 903 } 904 } 905 } 906 } 907 } else { 908 /* Newton step was not good; reduce the radius */ 909 tao->trust = bnk->alpha1 * PetscMin(bnk->dnorm, tao->trust); 910 } 911 break; 912 913 default: 914 if (stepType == BNK_NEWTON) { 915 if (prered < 0.0) { 916 /* The predicted reduction has the wrong sign. This cannot */ 917 /* happen in infinite precision arithmetic. Step should */ 918 /* be rejected! */ 919 tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm); 920 } else { 921 if (PetscIsInfOrNanReal(actred)) { 922 tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm); 923 } else { 924 if ((PetscAbsScalar(actred) <= bnk->epsilon) && (PetscAbsScalar(prered) <= bnk->epsilon)) { 925 kappa = 1.0; 926 } else { 927 kappa = actred / prered; 928 } 929 930 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx)); 931 tau_1 = bnk->theta * gdx / (bnk->theta * gdx - (1.0 - bnk->theta) * prered + actred); 932 tau_2 = bnk->theta * gdx / (bnk->theta * gdx + (1.0 + bnk->theta) * prered - actred); 933 tau_min = PetscMin(tau_1, tau_2); 934 tau_max = PetscMax(tau_1, tau_2); 935 936 if (kappa >= 1.0 - bnk->mu1) { 937 /* Great agreement */ 938 *accept = PETSC_TRUE; 939 if (tau_max < 1.0) { 940 tao->trust = PetscMax(tao->trust, bnk->gamma3 * bnk->dnorm); 941 } else if (tau_max > bnk->gamma4) { 942 tao->trust = PetscMax(tao->trust, bnk->gamma4 * bnk->dnorm); 943 } else { 944 tao->trust = PetscMax(tao->trust, tau_max * bnk->dnorm); 945 } 946 } else if (kappa >= 1.0 - bnk->mu2) { 947 /* Good agreement */ 948 *accept = PETSC_TRUE; 949 if (tau_max < bnk->gamma2) { 950 tao->trust = bnk->gamma2 * PetscMin(tao->trust, bnk->dnorm); 951 } else if (tau_max > bnk->gamma3) { 952 tao->trust = PetscMax(tao->trust, bnk->gamma3 * bnk->dnorm); 953 } else if (tau_max < 1.0) { 954 tao->trust = tau_max * PetscMin(tao->trust, bnk->dnorm); 955 } else { 956 tao->trust = PetscMax(tao->trust, tau_max * bnk->dnorm); 957 } 958 } else { 959 /* Not good agreement */ 960 if (tau_min > 1.0) { 961 tao->trust = bnk->gamma2 * PetscMin(tao->trust, bnk->dnorm); 962 } else if (tau_max < bnk->gamma1) { 963 tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm); 964 } else if ((tau_min < bnk->gamma1) && (tau_max >= 1.0)) { 965 tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm); 966 } else if ((tau_1 >= bnk->gamma1) && (tau_1 < 1.0) && ((tau_2 < bnk->gamma1) || (tau_2 >= 1.0))) { 967 tao->trust = tau_1 * PetscMin(tao->trust, bnk->dnorm); 968 } else if ((tau_2 >= bnk->gamma1) && (tau_2 < 1.0) && ((tau_1 < bnk->gamma1) || (tau_2 >= 1.0))) { 969 tao->trust = tau_2 * PetscMin(tao->trust, bnk->dnorm); 970 } else { 971 tao->trust = tau_max * PetscMin(tao->trust, bnk->dnorm); 972 } 973 } 974 } 975 } 976 } else { 977 /* Newton step was not good; reduce the radius */ 978 tao->trust = bnk->gamma1 * PetscMin(bnk->dnorm, tao->trust); 979 } 980 break; 981 } 982 /* Make sure the radius does not violate min and max settings */ 983 tao->trust = PetscMin(tao->trust, bnk->max_radius); 984 tao->trust = PetscMax(tao->trust, bnk->min_radius); 985 PetscFunctionReturn(0); 986 } 987 988 /* ---------------------------------------------------------- */ 989 990 PetscErrorCode TaoBNKAddStepCounts(Tao tao, PetscInt stepType) 991 { 992 TAO_BNK *bnk = (TAO_BNK *)tao->data; 993 994 PetscFunctionBegin; 995 switch (stepType) { 996 case BNK_NEWTON: 997 ++bnk->newt; 998 break; 999 case BNK_BFGS: 1000 ++bnk->bfgs; 1001 break; 1002 case BNK_SCALED_GRADIENT: 1003 ++bnk->sgrad; 1004 break; 1005 case BNK_GRADIENT: 1006 ++bnk->grad; 1007 break; 1008 default: 1009 break; 1010 } 1011 PetscFunctionReturn(0); 1012 } 1013 1014 /* ---------------------------------------------------------- */ 1015 1016 PetscErrorCode TaoSetUp_BNK(Tao tao) 1017 { 1018 TAO_BNK *bnk = (TAO_BNK *)tao->data; 1019 PetscInt i; 1020 1021 PetscFunctionBegin; 1022 if (!tao->gradient) { 1023 PetscCall(VecDuplicate(tao->solution,&tao->gradient)); 1024 } 1025 if (!tao->stepdirection) { 1026 PetscCall(VecDuplicate(tao->solution,&tao->stepdirection)); 1027 } 1028 if (!bnk->W) { 1029 PetscCall(VecDuplicate(tao->solution,&bnk->W)); 1030 } 1031 if (!bnk->Xold) { 1032 PetscCall(VecDuplicate(tao->solution,&bnk->Xold)); 1033 } 1034 if (!bnk->Gold) { 1035 PetscCall(VecDuplicate(tao->solution,&bnk->Gold)); 1036 } 1037 if (!bnk->Xwork) { 1038 PetscCall(VecDuplicate(tao->solution,&bnk->Xwork)); 1039 } 1040 if (!bnk->Gwork) { 1041 PetscCall(VecDuplicate(tao->solution,&bnk->Gwork)); 1042 } 1043 if (!bnk->unprojected_gradient) { 1044 PetscCall(VecDuplicate(tao->solution,&bnk->unprojected_gradient)); 1045 } 1046 if (!bnk->unprojected_gradient_old) { 1047 PetscCall(VecDuplicate(tao->solution,&bnk->unprojected_gradient_old)); 1048 } 1049 if (!bnk->Diag_min) { 1050 PetscCall(VecDuplicate(tao->solution,&bnk->Diag_min)); 1051 } 1052 if (!bnk->Diag_max) { 1053 PetscCall(VecDuplicate(tao->solution,&bnk->Diag_max)); 1054 } 1055 if (bnk->max_cg_its > 0) { 1056 /* Ensure that the important common vectors are shared between BNK and embedded BNCG */ 1057 bnk->bncg_ctx = (TAO_BNCG *)bnk->bncg->data; 1058 PetscCall(PetscObjectReference((PetscObject)(bnk->unprojected_gradient_old))); 1059 PetscCall(VecDestroy(&bnk->bncg_ctx->unprojected_gradient_old)); 1060 bnk->bncg_ctx->unprojected_gradient_old = bnk->unprojected_gradient_old; 1061 PetscCall(PetscObjectReference((PetscObject)(bnk->unprojected_gradient))); 1062 PetscCall(VecDestroy(&bnk->bncg_ctx->unprojected_gradient)); 1063 bnk->bncg_ctx->unprojected_gradient = bnk->unprojected_gradient; 1064 PetscCall(PetscObjectReference((PetscObject)(bnk->Gold))); 1065 PetscCall(VecDestroy(&bnk->bncg_ctx->G_old)); 1066 bnk->bncg_ctx->G_old = bnk->Gold; 1067 PetscCall(PetscObjectReference((PetscObject)(tao->gradient))); 1068 PetscCall(VecDestroy(&bnk->bncg->gradient)); 1069 bnk->bncg->gradient = tao->gradient; 1070 PetscCall(PetscObjectReference((PetscObject)(tao->stepdirection))); 1071 PetscCall(VecDestroy(&bnk->bncg->stepdirection)); 1072 bnk->bncg->stepdirection = tao->stepdirection; 1073 PetscCall(TaoSetSolution(bnk->bncg, tao->solution)); 1074 /* Copy over some settings from BNK into BNCG */ 1075 PetscCall(TaoSetMaximumIterations(bnk->bncg, bnk->max_cg_its)); 1076 PetscCall(TaoSetTolerances(bnk->bncg, tao->gatol, tao->grtol, tao->gttol)); 1077 PetscCall(TaoSetFunctionLowerBound(bnk->bncg, tao->fmin)); 1078 PetscCall(TaoSetConvergenceTest(bnk->bncg, tao->ops->convergencetest, tao->cnvP)); 1079 PetscCall(TaoSetObjective(bnk->bncg, tao->ops->computeobjective, tao->user_objP)); 1080 PetscCall(TaoSetGradient(bnk->bncg, NULL, tao->ops->computegradient, tao->user_gradP)); 1081 PetscCall(TaoSetObjectiveAndGradient(bnk->bncg, NULL, tao->ops->computeobjectiveandgradient, tao->user_objgradP)); 1082 PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)tao, (PetscObject)(bnk->bncg))); 1083 for (i=0; i<tao->numbermonitors; ++i) { 1084 PetscCall(TaoSetMonitor(bnk->bncg, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i])); 1085 PetscCall(PetscObjectReference((PetscObject)(tao->monitorcontext[i]))); 1086 } 1087 } 1088 bnk->X_inactive = NULL; 1089 bnk->G_inactive = NULL; 1090 bnk->inactive_work = NULL; 1091 bnk->active_work = NULL; 1092 bnk->inactive_idx = NULL; 1093 bnk->active_idx = NULL; 1094 bnk->active_lower = NULL; 1095 bnk->active_upper = NULL; 1096 bnk->active_fixed = NULL; 1097 bnk->M = NULL; 1098 bnk->H_inactive = NULL; 1099 bnk->Hpre_inactive = NULL; 1100 PetscFunctionReturn(0); 1101 } 1102 1103 /*------------------------------------------------------------*/ 1104 1105 PetscErrorCode TaoDestroy_BNK(Tao tao) 1106 { 1107 TAO_BNK *bnk = (TAO_BNK *)tao->data; 1108 1109 PetscFunctionBegin; 1110 PetscCall(VecDestroy(&bnk->W)); 1111 PetscCall(VecDestroy(&bnk->Xold)); 1112 PetscCall(VecDestroy(&bnk->Gold)); 1113 PetscCall(VecDestroy(&bnk->Xwork)); 1114 PetscCall(VecDestroy(&bnk->Gwork)); 1115 PetscCall(VecDestroy(&bnk->unprojected_gradient)); 1116 PetscCall(VecDestroy(&bnk->unprojected_gradient_old)); 1117 PetscCall(VecDestroy(&bnk->Diag_min)); 1118 PetscCall(VecDestroy(&bnk->Diag_max)); 1119 PetscCall(ISDestroy(&bnk->active_lower)); 1120 PetscCall(ISDestroy(&bnk->active_upper)); 1121 PetscCall(ISDestroy(&bnk->active_fixed)); 1122 PetscCall(ISDestroy(&bnk->active_idx)); 1123 PetscCall(ISDestroy(&bnk->inactive_idx)); 1124 PetscCall(MatDestroy(&bnk->Hpre_inactive)); 1125 PetscCall(MatDestroy(&bnk->H_inactive)); 1126 PetscCall(TaoDestroy(&bnk->bncg)); 1127 PetscCall(KSPDestroy(&tao->ksp)); 1128 PetscCall(PetscFree(tao->data)); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 /*------------------------------------------------------------*/ 1133 1134 PetscErrorCode TaoSetFromOptions_BNK(PetscOptionItems *PetscOptionsObject,Tao tao) 1135 { 1136 TAO_BNK *bnk = (TAO_BNK *)tao->data; 1137 1138 PetscFunctionBegin; 1139 PetscOptionsHeadBegin(PetscOptionsObject,"Newton-Krylov method for bound constrained optimization"); 1140 PetscCall(PetscOptionsEList("-tao_bnk_init_type", "radius initialization type", "", BNK_INIT, BNK_INIT_TYPES, BNK_INIT[bnk->init_type], &bnk->init_type, NULL)); 1141 PetscCall(PetscOptionsEList("-tao_bnk_update_type", "radius update type", "", BNK_UPDATE, BNK_UPDATE_TYPES, BNK_UPDATE[bnk->update_type], &bnk->update_type, NULL)); 1142 PetscCall(PetscOptionsEList("-tao_bnk_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, NULL)); 1143 PetscCall(PetscOptionsReal("-tao_bnk_sval", "(developer) Hessian perturbation starting value", "", bnk->sval, &bnk->sval,NULL)); 1144 PetscCall(PetscOptionsReal("-tao_bnk_imin", "(developer) minimum initial Hessian perturbation", "", bnk->imin, &bnk->imin,NULL)); 1145 PetscCall(PetscOptionsReal("-tao_bnk_imax", "(developer) maximum initial Hessian perturbation", "", bnk->imax, &bnk->imax,NULL)); 1146 PetscCall(PetscOptionsReal("-tao_bnk_imfac", "(developer) initial merit factor for Hessian perturbation", "", bnk->imfac, &bnk->imfac,NULL)); 1147 PetscCall(PetscOptionsReal("-tao_bnk_pmin", "(developer) minimum Hessian perturbation", "", bnk->pmin, &bnk->pmin,NULL)); 1148 PetscCall(PetscOptionsReal("-tao_bnk_pmax", "(developer) maximum Hessian perturbation", "", bnk->pmax, &bnk->pmax,NULL)); 1149 PetscCall(PetscOptionsReal("-tao_bnk_pgfac", "(developer) Hessian perturbation growth factor", "", bnk->pgfac, &bnk->pgfac,NULL)); 1150 PetscCall(PetscOptionsReal("-tao_bnk_psfac", "(developer) Hessian perturbation shrink factor", "", bnk->psfac, &bnk->psfac,NULL)); 1151 PetscCall(PetscOptionsReal("-tao_bnk_pmgfac", "(developer) merit growth factor for Hessian perturbation", "", bnk->pmgfac, &bnk->pmgfac,NULL)); 1152 PetscCall(PetscOptionsReal("-tao_bnk_pmsfac", "(developer) merit shrink factor for Hessian perturbation", "", bnk->pmsfac, &bnk->pmsfac,NULL)); 1153 PetscCall(PetscOptionsReal("-tao_bnk_eta1", "(developer) threshold for rejecting step (-tao_bnk_update_type reduction)", "", bnk->eta1, &bnk->eta1,NULL)); 1154 PetscCall(PetscOptionsReal("-tao_bnk_eta2", "(developer) threshold for accepting marginal step (-tao_bnk_update_type reduction)", "", bnk->eta2, &bnk->eta2,NULL)); 1155 PetscCall(PetscOptionsReal("-tao_bnk_eta3", "(developer) threshold for accepting reasonable step (-tao_bnk_update_type reduction)", "", bnk->eta3, &bnk->eta3,NULL)); 1156 PetscCall(PetscOptionsReal("-tao_bnk_eta4", "(developer) threshold for accepting good step (-tao_bnk_update_type reduction)", "", bnk->eta4, &bnk->eta4,NULL)); 1157 PetscCall(PetscOptionsReal("-tao_bnk_alpha1", "(developer) radius reduction factor for rejected step (-tao_bnk_update_type reduction)", "", bnk->alpha1, &bnk->alpha1,NULL)); 1158 PetscCall(PetscOptionsReal("-tao_bnk_alpha2", "(developer) radius reduction factor for marginally accepted bad step (-tao_bnk_update_type reduction)", "", bnk->alpha2, &bnk->alpha2,NULL)); 1159 PetscCall(PetscOptionsReal("-tao_bnk_alpha3", "(developer) radius increase factor for reasonable accepted step (-tao_bnk_update_type reduction)", "", bnk->alpha3, &bnk->alpha3,NULL)); 1160 PetscCall(PetscOptionsReal("-tao_bnk_alpha4", "(developer) radius increase factor for good accepted step (-tao_bnk_update_type reduction)", "", bnk->alpha4, &bnk->alpha4,NULL)); 1161 PetscCall(PetscOptionsReal("-tao_bnk_alpha5", "(developer) radius increase factor for very good accepted step (-tao_bnk_update_type reduction)", "", bnk->alpha5, &bnk->alpha5,NULL)); 1162 PetscCall(PetscOptionsReal("-tao_bnk_nu1", "(developer) threshold for small line-search step length (-tao_bnk_update_type step)", "", bnk->nu1, &bnk->nu1,NULL)); 1163 PetscCall(PetscOptionsReal("-tao_bnk_nu2", "(developer) threshold for reasonable line-search step length (-tao_bnk_update_type step)", "", bnk->nu2, &bnk->nu2,NULL)); 1164 PetscCall(PetscOptionsReal("-tao_bnk_nu3", "(developer) threshold for large line-search step length (-tao_bnk_update_type step)", "", bnk->nu3, &bnk->nu3,NULL)); 1165 PetscCall(PetscOptionsReal("-tao_bnk_nu4", "(developer) threshold for very large line-search step length (-tao_bnk_update_type step)", "", bnk->nu4, &bnk->nu4,NULL)); 1166 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)); 1167 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)); 1168 PetscCall(PetscOptionsReal("-tao_bnk_omega3", "(developer) radius factor for decent line-search step length (-tao_bnk_update_type step)", "", bnk->omega3, &bnk->omega3,NULL)); 1169 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)); 1170 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)); 1171 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)); 1172 PetscCall(PetscOptionsReal("-tao_bnk_mu2_i", "(developer) threshold for accepting good step (-tao_bnk_init_type interpolation)", "", bnk->mu2_i, &bnk->mu2_i,NULL)); 1173 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)); 1174 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)); 1175 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)); 1176 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)); 1177 PetscCall(PetscOptionsReal("-tao_bnk_theta_i", "(developer) trust region interpolation factor (-tao_bnk_init_type interpolation)", "", bnk->theta_i, &bnk->theta_i,NULL)); 1178 PetscCall(PetscOptionsReal("-tao_bnk_mu1", "(developer) threshold for accepting very good step (-tao_bnk_update_type interpolation)", "", bnk->mu1, &bnk->mu1,NULL)); 1179 PetscCall(PetscOptionsReal("-tao_bnk_mu2", "(developer) threshold for accepting good step (-tao_bnk_update_type interpolation)", "", bnk->mu2, &bnk->mu2,NULL)); 1180 PetscCall(PetscOptionsReal("-tao_bnk_gamma1", "(developer) radius reduction factor for rejected very bad step (-tao_bnk_update_type interpolation)", "", bnk->gamma1, &bnk->gamma1,NULL)); 1181 PetscCall(PetscOptionsReal("-tao_bnk_gamma2", "(developer) radius reduction factor for rejected bad step (-tao_bnk_update_type interpolation)", "", bnk->gamma2, &bnk->gamma2,NULL)); 1182 PetscCall(PetscOptionsReal("-tao_bnk_gamma3", "(developer) radius increase factor for accepted good step (-tao_bnk_update_type interpolation)", "", bnk->gamma3, &bnk->gamma3,NULL)); 1183 PetscCall(PetscOptionsReal("-tao_bnk_gamma4", "(developer) radius increase factor for accepted very good step (-tao_bnk_update_type interpolation)", "", bnk->gamma4, &bnk->gamma4,NULL)); 1184 PetscCall(PetscOptionsReal("-tao_bnk_theta", "(developer) trust region interpolation factor (-tao_bnk_update_type interpolation)", "", bnk->theta, &bnk->theta,NULL)); 1185 PetscCall(PetscOptionsReal("-tao_bnk_min_radius", "(developer) lower bound on initial radius", "", bnk->min_radius, &bnk->min_radius,NULL)); 1186 PetscCall(PetscOptionsReal("-tao_bnk_max_radius", "(developer) upper bound on radius", "", bnk->max_radius, &bnk->max_radius,NULL)); 1187 PetscCall(PetscOptionsReal("-tao_bnk_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL)); 1188 PetscCall(PetscOptionsReal("-tao_bnk_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol,NULL)); 1189 PetscCall(PetscOptionsReal("-tao_bnk_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step,NULL)); 1190 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)); 1191 PetscOptionsHeadEnd(); 1192 1193 PetscCall(TaoSetOptionsPrefix(bnk->bncg,((PetscObject)(tao))->prefix)); 1194 PetscCall(TaoAppendOptionsPrefix(bnk->bncg,"tao_bnk_cg_")); 1195 PetscCall(TaoSetFromOptions(bnk->bncg)); 1196 1197 PetscCall(KSPSetOptionsPrefix(tao->ksp,((PetscObject)(tao))->prefix)); 1198 PetscCall(KSPAppendOptionsPrefix(tao->ksp,"tao_bnk_")); 1199 PetscCall(KSPSetFromOptions(tao->ksp)); 1200 PetscFunctionReturn(0); 1201 } 1202 1203 /*------------------------------------------------------------*/ 1204 1205 PetscErrorCode TaoView_BNK(Tao tao, PetscViewer viewer) 1206 { 1207 TAO_BNK *bnk = (TAO_BNK *)tao->data; 1208 PetscInt nrejects; 1209 PetscBool isascii; 1210 1211 PetscFunctionBegin; 1212 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 1213 if (isascii) { 1214 PetscCall(PetscViewerASCIIPushTab(viewer)); 1215 if (bnk->M) { 1216 PetscCall(MatLMVMGetRejectCount(bnk->M,&nrejects)); 1217 PetscCall(PetscViewerASCIIPrintf(viewer, "Rejected BFGS updates: %" PetscInt_FMT "\n",nrejects)); 1218 } 1219 PetscCall(PetscViewerASCIIPrintf(viewer, "CG steps: %" PetscInt_FMT "\n", bnk->tot_cg_its)); 1220 PetscCall(PetscViewerASCIIPrintf(viewer, "Newton steps: %" PetscInt_FMT "\n", bnk->newt)); 1221 if (bnk->M) { 1222 PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS steps: %" PetscInt_FMT "\n", bnk->bfgs)); 1223 } 1224 PetscCall(PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %" PetscInt_FMT "\n", bnk->sgrad)); 1225 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", bnk->grad)); 1226 PetscCall(PetscViewerASCIIPrintf(viewer, "KSP termination reasons:\n")); 1227 PetscCall(PetscViewerASCIIPrintf(viewer, " atol: %" PetscInt_FMT "\n", bnk->ksp_atol)); 1228 PetscCall(PetscViewerASCIIPrintf(viewer, " rtol: %" PetscInt_FMT "\n", bnk->ksp_rtol)); 1229 PetscCall(PetscViewerASCIIPrintf(viewer, " ctol: %" PetscInt_FMT "\n", bnk->ksp_ctol)); 1230 PetscCall(PetscViewerASCIIPrintf(viewer, " negc: %" PetscInt_FMT "\n", bnk->ksp_negc)); 1231 PetscCall(PetscViewerASCIIPrintf(viewer, " dtol: %" PetscInt_FMT "\n", bnk->ksp_dtol)); 1232 PetscCall(PetscViewerASCIIPrintf(viewer, " iter: %" PetscInt_FMT "\n", bnk->ksp_iter)); 1233 PetscCall(PetscViewerASCIIPrintf(viewer, " othr: %" PetscInt_FMT "\n", bnk->ksp_othr)); 1234 PetscCall(PetscViewerASCIIPopTab(viewer)); 1235 } 1236 PetscFunctionReturn(0); 1237 } 1238 1239 /* ---------------------------------------------------------- */ 1240 1241 /*MC 1242 TAOBNK - Shared base-type for Bounded Newton-Krylov type algorithms. 1243 At each iteration, the BNK methods solve the symmetric 1244 system of equations to obtain the step diretion dk: 1245 Hk dk = -gk 1246 for free variables only. The step can be globalized either through 1247 trust-region methods, or a line search, or a heuristic mixture of both. 1248 1249 Options Database Keys: 1250 + -tao_bnk_max_cg_its - maximum number of bounded conjugate-gradient iterations taken in each Newton loop 1251 . -tao_bnk_init_type - trust radius initialization method ("constant", "direction", "interpolation") 1252 . -tao_bnk_update_type - trust radius update method ("step", "direction", "interpolation") 1253 . -tao_bnk_as_type - active-set estimation method ("none", "bertsekas") 1254 . -tao_bnk_as_tol - (developer) initial tolerance used in estimating bounded active variables (-as_type bertsekas) 1255 . -tao_bnk_as_step - (developer) trial step length used in estimating bounded active variables (-as_type bertsekas) 1256 . -tao_bnk_sval - (developer) Hessian perturbation starting value 1257 . -tao_bnk_imin - (developer) minimum initial Hessian perturbation 1258 . -tao_bnk_imax - (developer) maximum initial Hessian perturbation 1259 . -tao_bnk_pmin - (developer) minimum Hessian perturbation 1260 . -tao_bnk_pmax - (developer) aximum Hessian perturbation 1261 . -tao_bnk_pgfac - (developer) Hessian perturbation growth factor 1262 . -tao_bnk_psfac - (developer) Hessian perturbation shrink factor 1263 . -tao_bnk_imfac - (developer) initial merit factor for Hessian perturbation 1264 . -tao_bnk_pmgfac - (developer) merit growth factor for Hessian perturbation 1265 . -tao_bnk_pmsfac - (developer) merit shrink factor for Hessian perturbation 1266 . -tao_bnk_eta1 - (developer) threshold for rejecting step (-update_type reduction) 1267 . -tao_bnk_eta2 - (developer) threshold for accepting marginal step (-update_type reduction) 1268 . -tao_bnk_eta3 - (developer) threshold for accepting reasonable step (-update_type reduction) 1269 . -tao_bnk_eta4 - (developer) threshold for accepting good step (-update_type reduction) 1270 . -tao_bnk_alpha1 - (developer) radius reduction factor for rejected step (-update_type reduction) 1271 . -tao_bnk_alpha2 - (developer) radius reduction factor for marginally accepted bad step (-update_type reduction) 1272 . -tao_bnk_alpha3 - (developer) radius increase factor for reasonable accepted step (-update_type reduction) 1273 . -tao_bnk_alpha4 - (developer) radius increase factor for good accepted step (-update_type reduction) 1274 . -tao_bnk_alpha5 - (developer) radius increase factor for very good accepted step (-update_type reduction) 1275 . -tao_bnk_epsilon - (developer) tolerance for small pred/actual ratios that trigger automatic step acceptance (-update_type reduction) 1276 . -tao_bnk_mu1 - (developer) threshold for accepting very good step (-update_type interpolation) 1277 . -tao_bnk_mu2 - (developer) threshold for accepting good step (-update_type interpolation) 1278 . -tao_bnk_gamma1 - (developer) radius reduction factor for rejected very bad step (-update_type interpolation) 1279 . -tao_bnk_gamma2 - (developer) radius reduction factor for rejected bad step (-update_type interpolation) 1280 . -tao_bnk_gamma3 - (developer) radius increase factor for accepted good step (-update_type interpolation) 1281 . -tao_bnk_gamma4 - (developer) radius increase factor for accepted very good step (-update_type interpolation) 1282 . -tao_bnk_theta - (developer) trust region interpolation factor (-update_type interpolation) 1283 . -tao_bnk_nu1 - (developer) threshold for small line-search step length (-update_type step) 1284 . -tao_bnk_nu2 - (developer) threshold for reasonable line-search step length (-update_type step) 1285 . -tao_bnk_nu3 - (developer) threshold for large line-search step length (-update_type step) 1286 . -tao_bnk_nu4 - (developer) threshold for very large line-search step length (-update_type step) 1287 . -tao_bnk_omega1 - (developer) radius reduction factor for very small line-search step length (-update_type step) 1288 . -tao_bnk_omega2 - (developer) radius reduction factor for small line-search step length (-update_type step) 1289 . -tao_bnk_omega3 - (developer) radius factor for decent line-search step length (-update_type step) 1290 . -tao_bnk_omega4 - (developer) radius increase factor for large line-search step length (-update_type step) 1291 . -tao_bnk_omega5 - (developer) radius increase factor for very large line-search step length (-update_type step) 1292 . -tao_bnk_mu1_i - (developer) threshold for accepting very good step (-init_type interpolation) 1293 . -tao_bnk_mu2_i - (developer) threshold for accepting good step (-init_type interpolation) 1294 . -tao_bnk_gamma1_i - (developer) radius reduction factor for rejected very bad step (-init_type interpolation) 1295 . -tao_bnk_gamma2_i - (developer) radius reduction factor for rejected bad step (-init_type interpolation) 1296 . -tao_bnk_gamma3_i - (developer) radius increase factor for accepted good step (-init_type interpolation) 1297 . -tao_bnk_gamma4_i - (developer) radius increase factor for accepted very good step (-init_type interpolation) 1298 - -tao_bnk_theta_i - (developer) trust region interpolation factor (-init_type interpolation) 1299 1300 Level: beginner 1301 M*/ 1302 1303 PetscErrorCode TaoCreate_BNK(Tao tao) 1304 { 1305 TAO_BNK *bnk; 1306 PC pc; 1307 1308 PetscFunctionBegin; 1309 PetscCall(PetscNewLog(tao,&bnk)); 1310 1311 tao->ops->setup = TaoSetUp_BNK; 1312 tao->ops->view = TaoView_BNK; 1313 tao->ops->setfromoptions = TaoSetFromOptions_BNK; 1314 tao->ops->destroy = TaoDestroy_BNK; 1315 1316 /* Override default settings (unless already changed) */ 1317 if (!tao->max_it_changed) tao->max_it = 50; 1318 if (!tao->trust0_changed) tao->trust0 = 100.0; 1319 1320 tao->data = (void*)bnk; 1321 1322 /* Hessian shifting parameters */ 1323 bnk->computehessian = TaoBNKComputeHessian; 1324 bnk->computestep = TaoBNKComputeStep; 1325 1326 bnk->sval = 0.0; 1327 bnk->imin = 1.0e-4; 1328 bnk->imax = 1.0e+2; 1329 bnk->imfac = 1.0e-1; 1330 1331 bnk->pmin = 1.0e-12; 1332 bnk->pmax = 1.0e+2; 1333 bnk->pgfac = 1.0e+1; 1334 bnk->psfac = 4.0e-1; 1335 bnk->pmgfac = 1.0e-1; 1336 bnk->pmsfac = 1.0e-1; 1337 1338 /* Default values for trust-region radius update based on steplength */ 1339 bnk->nu1 = 0.25; 1340 bnk->nu2 = 0.50; 1341 bnk->nu3 = 1.00; 1342 bnk->nu4 = 1.25; 1343 1344 bnk->omega1 = 0.25; 1345 bnk->omega2 = 0.50; 1346 bnk->omega3 = 1.00; 1347 bnk->omega4 = 2.00; 1348 bnk->omega5 = 4.00; 1349 1350 /* Default values for trust-region radius update based on reduction */ 1351 bnk->eta1 = 1.0e-4; 1352 bnk->eta2 = 0.25; 1353 bnk->eta3 = 0.50; 1354 bnk->eta4 = 0.90; 1355 1356 bnk->alpha1 = 0.25; 1357 bnk->alpha2 = 0.50; 1358 bnk->alpha3 = 1.00; 1359 bnk->alpha4 = 2.00; 1360 bnk->alpha5 = 4.00; 1361 1362 /* Default values for trust-region radius update based on interpolation */ 1363 bnk->mu1 = 0.10; 1364 bnk->mu2 = 0.50; 1365 1366 bnk->gamma1 = 0.25; 1367 bnk->gamma2 = 0.50; 1368 bnk->gamma3 = 2.00; 1369 bnk->gamma4 = 4.00; 1370 1371 bnk->theta = 0.05; 1372 1373 /* Default values for trust region initialization based on interpolation */ 1374 bnk->mu1_i = 0.35; 1375 bnk->mu2_i = 0.50; 1376 1377 bnk->gamma1_i = 0.0625; 1378 bnk->gamma2_i = 0.5; 1379 bnk->gamma3_i = 2.0; 1380 bnk->gamma4_i = 5.0; 1381 1382 bnk->theta_i = 0.25; 1383 1384 /* Remaining parameters */ 1385 bnk->max_cg_its = 0; 1386 bnk->min_radius = 1.0e-10; 1387 bnk->max_radius = 1.0e10; 1388 bnk->epsilon = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0); 1389 bnk->as_tol = 1.0e-3; 1390 bnk->as_step = 1.0e-3; 1391 bnk->dmin = 1.0e-6; 1392 bnk->dmax = 1.0e6; 1393 1394 bnk->M = NULL; 1395 bnk->bfgs_pre = NULL; 1396 bnk->init_type = BNK_INIT_INTERPOLATION; 1397 bnk->update_type = BNK_UPDATE_REDUCTION; 1398 bnk->as_type = BNK_AS_BERTSEKAS; 1399 1400 /* Create the embedded BNCG solver */ 1401 PetscCall(TaoCreate(PetscObjectComm((PetscObject)tao),&bnk->bncg)); 1402 PetscCall(PetscObjectIncrementTabLevel((PetscObject)bnk->bncg,(PetscObject)tao,1)); 1403 PetscCall(TaoSetType(bnk->bncg,TAOBNCG)); 1404 1405 /* Create the line search */ 1406 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch)); 1407 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch,(PetscObject)tao,1)); 1408 PetscCall(TaoLineSearchSetType(tao->linesearch,TAOLINESEARCHMT)); 1409 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch,tao)); 1410 1411 /* Set linear solver to default for symmetric matrices */ 1412 PetscCall(KSPCreate(((PetscObject)tao)->comm,&tao->ksp)); 1413 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1)); 1414 PetscCall(KSPSetType(tao->ksp,KSPSTCG)); 1415 PetscCall(KSPGetPC(tao->ksp,&pc)); 1416 PetscCall(PCSetType(pc,PCLMVM)); 1417 PetscFunctionReturn(0); 1418 } 1419