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