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