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