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