xref: /petsc/src/tao/bound/impls/bnk/bnk.c (revision 8d5ead369b7916272d211d106c850283b970bae2) 
1 #include <petsctaolinesearch.h>
2 #include <../src/tao/bound/impls/bnk/bnk.h>
3 
4 #include <petscksp.h>
5 
6 /* Routine for BFGS preconditioner */
7 
8 PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x)
9 {
10   PetscErrorCode ierr;
11   Mat            M;
12 
13   PetscFunctionBegin;
14   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
15   PetscValidHeaderSpecific(b,VEC_CLASSID,2);
16   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
17   ierr = PCShellGetContext(pc,(void**)&M);CHKERRQ(ierr);
18   ierr = MatLMVMSolveInactive(M, b, x);CHKERRQ(ierr);
19   PetscFunctionReturn(0);
20 }
21 
22 /*------------------------------------------------------------*/
23 
24 /* Routine for initializing the KSP solver, the BFGS preconditioner, and the initial trust radius estimation */
25 
26 PetscErrorCode TaoBNKInitialize(Tao tao)
27 {
28   PetscErrorCode               ierr;
29   TAO_BNK                      *bnk = (TAO_BNK *)tao->data;
30   KSPType                      ksp_type;
31   PC                           pc;
32 
33   PetscReal                    fmin, ftrial, prered, actred, kappa, sigma, dnorm;
34   PetscReal                    tau, tau_1, tau_2, tau_max, tau_min, max_radius;
35   PetscReal                    delta, step = 1.0;
36 
37   PetscInt                     n,N,needH = 1;
38 
39   PetscInt                     i_max = 5;
40   PetscInt                     j_max = 1;
41   PetscInt                     i, j;
42 
43   PetscFunctionBegin;
44   /* Number of times ksp stopped because of these reasons */
45   bnk->ksp_atol = 0;
46   bnk->ksp_rtol = 0;
47   bnk->ksp_dtol = 0;
48   bnk->ksp_ctol = 0;
49   bnk->ksp_negc = 0;
50   bnk->ksp_iter = 0;
51   bnk->ksp_othr = 0;
52 
53   /* Initialize the Hessian perturbation */
54   bnk->pert = bnk->sval;
55 
56   /* Initialize trust-region radius when using nash, stcg, or gltr
57      Command automatically ignored for other methods
58      Will be reset during the first iteration
59   */
60   ierr = KSPGetType(tao->ksp,&ksp_type);CHKERRQ(ierr);
61   ierr = PetscStrcmp(ksp_type,KSPCGNASH,&bnk->is_nash);CHKERRQ(ierr);
62   ierr = PetscStrcmp(ksp_type,KSPCGSTCG,&bnk->is_stcg);CHKERRQ(ierr);
63   ierr = PetscStrcmp(ksp_type,KSPCGGLTR,&bnk->is_gltr);CHKERRQ(ierr);
64 
65   ierr = KSPCGSetRadius(tao->ksp,bnk->max_radius);CHKERRQ(ierr);
66 
67   if (bnk->is_nash || bnk->is_stcg || bnk->is_gltr) {
68     if (tao->trust0 < 0.0) SETERRQ(PETSC_COMM_SELF,1,"Initial radius negative");
69     tao->trust = tao->trust0;
70     tao->trust = PetscMax(tao->trust, bnk->min_radius);
71     tao->trust = PetscMin(tao->trust, bnk->max_radius);
72   }
73 
74   /* Get vectors we will need */
75   if (BNK_PC_BFGS == bnk->pc_type && !bnk->M) {
76     ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr);
77     ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);
78     ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&bnk->M);CHKERRQ(ierr);
79     ierr = MatLMVMAllocateVectors(bnk->M,tao->solution);CHKERRQ(ierr);
80   }
81 
82   /* create vectors for the limited memory preconditioner */
83   if ((BNK_PC_BFGS == bnk->pc_type) && (BFGS_SCALE_BFGS != bnk->bfgs_scale_type)) {
84     if (!bnk->Diag) {
85       ierr = VecDuplicate(tao->solution,&bnk->Diag);CHKERRQ(ierr);
86     }
87   }
88 
89   /* Modify the preconditioner to use the bfgs approximation */
90   ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr);
91   switch(bnk->pc_type) {
92   case BNK_PC_NONE:
93     ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr);
94     ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
95     break;
96 
97   case BNK_PC_AHESS:
98     ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr);
99     ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
100     ierr = PCJacobiSetUseAbs(pc,PETSC_TRUE);CHKERRQ(ierr);
101     break;
102 
103   case BNK_PC_BFGS:
104     ierr = PCSetType(pc, PCSHELL);CHKERRQ(ierr);
105     ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
106     ierr = PCShellSetName(pc, "bfgs");CHKERRQ(ierr);
107     ierr = PCShellSetContext(pc, bnk->M);CHKERRQ(ierr);
108     ierr = PCShellSetApply(pc, MatLMVMSolveShell);CHKERRQ(ierr);
109     break;
110 
111   default:
112     /* Use the pc method set by pc_type */
113     break;
114   }
115 
116   /* Initialize trust-region radius.  The initialization is only performed
117      when we are using Nash, Steihaug-Toint or the Generalized Lanczos method. */
118   if (bnk->is_nash || bnk->is_stcg || bnk->is_gltr) {
119     switch(bnk->init_type) {
120     case BNK_INIT_CONSTANT:
121       /* Use the initial radius specified */
122       break;
123 
124     case BNK_INIT_INTERPOLATION:
125       /* Use the initial radius specified */
126       max_radius = 0.0;
127 
128       for (j = 0; j < j_max; ++j) {
129         fmin = bnk->f;
130         sigma = 0.0;
131 
132         if (needH) {
133           /* Compute the Hessian at the new step, and extract the inactive subsystem */
134           ierr = TaoComputeHessian(tao, tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
135           ierr = ISDestroy(&bnk->inactive_idx);CHKERRQ(ierr);
136           ierr = VecWhichInactive(tao->XL,tao->solution,bnk->unprojected_gradient,tao->XU,PETSC_TRUE,&bnk->inactive_idx);CHKERRQ(ierr);
137           ierr = TaoMatGetSubMat(tao->hessian, bnk->inactive_idx, bnk->Xwork, TAO_SUBSET_MASK, &bnk->H_inactive);CHKERRQ(ierr);
138           needH = 0;
139         }
140 
141         for (i = 0; i < i_max; ++i) {
142           /* Take a steepest descent step and snap it to bounds */
143           ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr);
144           ierr = VecAXPY(tao->solution, -tao->trust/bnk->gnorm, tao->gradient);CHKERRQ(ierr);
145           ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr);
146           /* Adjust the trust radius if the bound snapping changed the step */
147           ierr = VecCopy(tao->solution, bnk->W);CHKERRQ(ierr);
148           ierr = VecAXPY(bnk->W, -1.0, bnk->Xold);CHKERRQ(ierr);
149           ierr = VecNorm(bnk->W, NORM_2, &dnorm);CHKERRQ(ierr);
150           if (dnorm != tao->trust) tao->trust = dnorm;
151           /* Compute the objective at the trial */
152           ierr = TaoComputeObjective(tao, tao->solution, &ftrial);CHKERRQ(ierr);
153           ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr);
154           if (PetscIsInfOrNanReal(ftrial)) {
155             tau = bnk->gamma1_i;
156           } else {
157             if (ftrial < fmin) {
158               fmin = ftrial;
159               sigma = -tao->trust / bnk->gnorm;
160             }
161 
162             ierr = MatMult(bnk->H_inactive, tao->gradient, bnk->W);CHKERRQ(ierr);
163             ierr = VecDot(tao->gradient, bnk->W, &prered);CHKERRQ(ierr);
164 
165             prered = tao->trust * (bnk->gnorm - 0.5 * tao->trust * prered / (bnk->gnorm * bnk->gnorm));
166             actred = bnk->f - ftrial;
167             if ((PetscAbsScalar(actred) <= bnk->epsilon) && (PetscAbsScalar(prered) <= bnk->epsilon)) {
168               kappa = 1.0;
169             } else {
170               kappa = actred / prered;
171             }
172 
173             tau_1 = bnk->theta_i * bnk->gnorm * tao->trust / (bnk->theta_i * bnk->gnorm * tao->trust + (1.0 - bnk->theta_i) * prered - actred);
174             tau_2 = bnk->theta_i * bnk->gnorm * tao->trust / (bnk->theta_i * bnk->gnorm * tao->trust - (1.0 + bnk->theta_i) * prered + actred);
175             tau_min = PetscMin(tau_1, tau_2);
176             tau_max = PetscMax(tau_1, tau_2);
177 
178             if (PetscAbsScalar(kappa - 1.0) <= bnk->mu1_i) {
179               /* Great agreement */
180               max_radius = PetscMax(max_radius, tao->trust);
181 
182               if (tau_max < 1.0) {
183                 tau = bnk->gamma3_i;
184               } else if (tau_max > bnk->gamma4_i) {
185                 tau = bnk->gamma4_i;
186               } else if (tau_1 >= 1.0 && tau_1 <= bnk->gamma4_i && tau_2 < 1.0) {
187                 tau = tau_1;
188               } else if (tau_2 >= 1.0 && tau_2 <= bnk->gamma4_i && tau_1 < 1.0) {
189                 tau = tau_2;
190               } else {
191                 tau = tau_max;
192               }
193             } else if (PetscAbsScalar(kappa - 1.0) <= bnk->mu2_i) {
194               /* Good agreement */
195               max_radius = PetscMax(max_radius, tao->trust);
196 
197               if (tau_max < bnk->gamma2_i) {
198                 tau = bnk->gamma2_i;
199               } else if (tau_max > bnk->gamma3_i) {
200                 tau = bnk->gamma3_i;
201               } else {
202                 tau = tau_max;
203               }
204             } else {
205               /* Not good agreement */
206               if (tau_min > 1.0) {
207                 tau = bnk->gamma2_i;
208               } else if (tau_max < bnk->gamma1_i) {
209                 tau = bnk->gamma1_i;
210               } else if ((tau_min < bnk->gamma1_i) && (tau_max >= 1.0)) {
211                 tau = bnk->gamma1_i;
212               } else if ((tau_1 >= bnk->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < bnk->gamma1_i) || (tau_2 >= 1.0))) {
213                 tau = tau_1;
214               } else if ((tau_2 >= bnk->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < bnk->gamma1_i) || (tau_2 >= 1.0))) {
215                 tau = tau_2;
216               } else {
217                 tau = tau_max;
218               }
219             }
220           }
221           tao->trust = tau * tao->trust;
222         }
223 
224         if (fmin < bnk->f) {
225           bnk->f = fmin;
226           ierr = VecAXPY(tao->solution,sigma,tao->gradient);CHKERRQ(ierr);
227           ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr);
228           ierr = TaoComputeGradient(tao,tao->solution,bnk->unprojected_gradient);CHKERRQ(ierr);
229           ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr);
230 
231           ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr);
232           if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute gradient generated Inf or NaN");
233           needH = 1;
234           ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr);
235 
236           ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
237           ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,step);CHKERRQ(ierr);
238           ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
239           if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
240         }
241       }
242       tao->trust = PetscMax(tao->trust, max_radius);
243 
244       /* Modify the radius if it is too large or small */
245       tao->trust = PetscMax(tao->trust, bnk->min_radius);
246       tao->trust = PetscMin(tao->trust, bnk->max_radius);
247       break;
248 
249     default:
250       /* Norm of the first direction will initialize radius */
251       tao->trust = 0.0;
252       break;
253     }
254   }
255 
256   /* Set initial scaling for the BFGS preconditioner
257      This step is done after computing the initial trust-region radius
258      since the function value may have decreased */
259   if (BNK_PC_BFGS == bnk->pc_type) {
260     if (bnk->f != 0.0) {
261       delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm);
262     } else {
263       delta = 2.0 / (bnk->gnorm*bnk->gnorm);
264     }
265     ierr = MatLMVMSetDelta(bnk->M,delta);CHKERRQ(ierr);
266   }
267 
268   /* Set counter for gradient/reset steps*/
269   bnk->newt = 0;
270   bnk->bfgs = 0;
271   bnk->sgrad = 0;
272   bnk->grad = 0;
273   PetscFunctionReturn(0);
274 }
275 
276 /*------------------------------------------------------------*/
277 
278 /* Routine for computing the Newton step.
279 
280   If the safeguard is enabled, the Newton step is verified to be a
281   descent direction, with fallbacks onto BFGS, scaled gradient, and unscaled
282   gradient steps if/when necessary.
283 
284   The function reports back on which type of step has ultimately been stored
285   under tao->stepdirection.
286 */
287 
288 PetscErrorCode TaoBNKComputeStep(Tao tao, PetscBool safeguard, PetscInt *stepType)
289 {
290   PetscErrorCode               ierr;
291   TAO_BNK                      *bnk = (TAO_BNK *)tao->data;
292   KSPConvergedReason           ksp_reason;
293 
294   Vec                          active_step;
295   PetscReal                    gdx, delta, e_min;
296   PetscInt                     bfgsUpdates = 0;
297   PetscInt                     kspits;
298 
299   PetscFunctionBegin;
300   /* Shift the Hessian matrix */
301   if (bnk->pert > 0) {
302     ierr = MatShift(tao->hessian, bnk->pert);CHKERRQ(ierr);
303     if (tao->hessian != tao->hessian_pre) {
304       ierr = MatShift(tao->hessian_pre, bnk->pert);CHKERRQ(ierr);
305     }
306   }
307 
308   /* Determine the inactive and active sets */
309   ierr = ISDestroy(&bnk->inactive_idx);CHKERRQ(ierr);
310   ierr = ISDestroy(&bnk->active_idx);CHKERRQ(ierr);
311   ierr = VecWhichInactive(tao->XL,tao->solution,bnk->unprojected_gradient,tao->XU,PETSC_TRUE,&bnk->inactive_idx);CHKERRQ(ierr);
312   ierr = ISComplementVec(bnk->inactive_idx, tao->solution, &bnk->active_idx);CHKERRQ(ierr);
313 
314   /* Prepare masked matrices for the inactive set */
315   ierr = MatLMVMSetInactive(bnk->M, bnk->inactive_idx);CHKERRQ(ierr);
316   ierr = VecSet(tao->stepdirection, 0.0);CHKERRQ(ierr);
317   ierr = TaoMatGetSubMat(tao->hessian, bnk->inactive_idx, bnk->Xwork, TAO_SUBSET_MASK, &bnk->H_inactive);CHKERRQ(ierr);
318   if (tao->hessian == tao->hessian_pre) {
319     bnk->Hpre_inactive = bnk->H_inactive;
320   } else {
321     ierr = TaoMatGetSubMat(tao->hessian_pre, bnk->inactive_idx, bnk->Xwork, TAO_SUBSET_MASK, &bnk->Hpre_inactive);CHKERRQ(ierr);
322   }
323 
324   /* Solve the Newton system of equations */
325   ierr = KSPSetOperators(tao->ksp,bnk->H_inactive,bnk->Hpre_inactive);CHKERRQ(ierr);
326   if (bnk->is_nash || bnk->is_stcg || bnk->is_gltr) {
327     ierr = KSPCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
328     ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
329     ierr = KSPGetIterationNumber(tao->ksp,&kspits);CHKERRQ(ierr);
330     tao->ksp_its+=kspits;
331     tao->ksp_tot_its+=kspits;
332     ierr = KSPCGGetNormD(tao->ksp,&bnk->dnorm);CHKERRQ(ierr);
333 
334     if (0.0 == tao->trust) {
335       /* Radius was uninitialized; use the norm of the direction */
336       if (bnk->dnorm > 0.0) {
337         tao->trust = bnk->dnorm;
338 
339         /* Modify the radius if it is too large or small */
340         tao->trust = PetscMax(tao->trust, bnk->min_radius);
341         tao->trust = PetscMin(tao->trust, bnk->max_radius);
342       } else {
343         /* The direction was bad; set radius to default value and re-solve
344            the trust-region subproblem to get a direction */
345         tao->trust = tao->trust0;
346 
347         /* Modify the radius if it is too large or small */
348         tao->trust = PetscMax(tao->trust, bnk->min_radius);
349         tao->trust = PetscMin(tao->trust, bnk->max_radius);
350 
351         ierr = KSPCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
352         ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
353         ierr = KSPGetIterationNumber(tao->ksp,&kspits);CHKERRQ(ierr);
354         tao->ksp_its+=kspits;
355         tao->ksp_tot_its+=kspits;
356         ierr = KSPCGGetNormD(tao->ksp,&bnk->dnorm);CHKERRQ(ierr);
357 
358         if (bnk->dnorm == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
359       }
360     }
361   } else {
362     ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
363     ierr = KSPGetIterationNumber(tao->ksp, &kspits);CHKERRQ(ierr);
364     tao->ksp_its += kspits;
365     tao->ksp_tot_its+=kspits;
366   }
367   ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
368   ierr = KSPGetConvergedReason(tao->ksp, &ksp_reason);CHKERRQ(ierr);
369 
370   /* Make sure the BFGS preconditioner is healthy */
371   if (bnk->pc_type == BNK_PC_BFGS) {
372     ierr = MatLMVMGetUpdates(bnk->M, &bfgsUpdates);CHKERRQ(ierr);
373     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (bfgsUpdates > 1)) {
374       /* Preconditioner is numerically indefinite; reset the approximation. */
375       if (bnk->f != 0.0) {
376         delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm);
377       } else {
378         delta = 2.0 / (bnk->gnorm*bnk->gnorm);
379       }
380       ierr = MatLMVMSetDelta(bnk->M,delta);CHKERRQ(ierr);
381       ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr);
382       ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr);
383       bfgsUpdates = 1;
384     }
385   }
386 
387   if (KSP_CONVERGED_ATOL == ksp_reason) {
388     ++bnk->ksp_atol;
389   } else if (KSP_CONVERGED_RTOL == ksp_reason) {
390     ++bnk->ksp_rtol;
391   } else if (KSP_CONVERGED_CG_CONSTRAINED == ksp_reason) {
392     ++bnk->ksp_ctol;
393   } else if (KSP_CONVERGED_CG_NEG_CURVE == ksp_reason) {
394     ++bnk->ksp_negc;
395   } else if (KSP_DIVERGED_DTOL == ksp_reason) {
396     ++bnk->ksp_dtol;
397   } else if (KSP_DIVERGED_ITS == ksp_reason) {
398     ++bnk->ksp_iter;
399   } else {
400     ++bnk->ksp_othr;
401   }
402 
403   /* Check for success (descent direction) */
404   if (safeguard) {
405     ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr);
406     if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
407       /* Newton step is not descent or direction produced Inf or NaN
408          Update the perturbation for next time */
409       if (bnk->pert <= 0.0) {
410         /* Initialize the perturbation */
411         bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm));
412         if (bnk->is_gltr) {
413           ierr = KSPCGGLTRGetMinEig(tao->ksp,&e_min);CHKERRQ(ierr);
414           bnk->pert = PetscMax(bnk->pert, -e_min);
415         }
416       } else {
417         /* Increase the perturbation */
418         bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm));
419       }
420 
421       if (BNK_PC_BFGS != bnk->pc_type) {
422         /* We don't have the bfgs matrix around and updated
423            Must use gradient direction in this case */
424         ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
425         ++bnk->grad;
426         *stepType = BNK_GRADIENT;
427       } else {
428         /* Attempt to use the BFGS direction */
429         ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr);
430 
431         /* Check for success (descent direction)
432            NOTE: Negative gdx here means not a descent direction because
433            the fall-back step is missing a negative sign. */
434         ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr);
435         if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
436           /* BFGS direction is not descent or direction produced not a number
437              We can assert bfgsUpdates > 1 in this case because
438              the first solve produces the scaled gradient direction,
439              which is guaranteed to be descent */
440 
441           /* Use steepest descent direction (scaled) */
442           if (bnk->f != 0.0) {
443             delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm);
444           } else {
445             delta = 2.0 / (bnk->gnorm*bnk->gnorm);
446           }
447           ierr = MatLMVMSetDelta(bnk->M, delta);CHKERRQ(ierr);
448           ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr);
449           ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr);
450           ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr);
451 
452           bfgsUpdates = 1;
453           ++bnk->sgrad;
454           *stepType = BNK_SCALED_GRADIENT;
455         } else {
456           if (1 == bfgsUpdates) {
457             /* The first BFGS direction is always the scaled gradient */
458             ++bnk->sgrad;
459             *stepType = BNK_SCALED_GRADIENT;
460           } else {
461             ++bnk->bfgs;
462             *stepType = BNK_BFGS;
463           }
464         }
465       }
466       /* Make sure the safeguarded fall-back step is zero for actively bounded variables */
467       ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
468       ierr = VecGetSubVector(tao->stepdirection, bnk->active_idx, &active_step);CHKERRQ(ierr);
469       ierr = VecSet(active_step, 0.0);CHKERRQ(ierr);
470       ierr = VecRestoreSubVector(tao->stepdirection, bnk->active_idx, &active_step);CHKERRQ(ierr);
471     } else {
472       /* Computed Newton step is descent */
473       switch (ksp_reason) {
474       case KSP_DIVERGED_NANORINF:
475       case KSP_DIVERGED_BREAKDOWN:
476       case KSP_DIVERGED_INDEFINITE_MAT:
477       case KSP_DIVERGED_INDEFINITE_PC:
478       case KSP_CONVERGED_CG_NEG_CURVE:
479         /* Matrix or preconditioner is indefinite; increase perturbation */
480         if (bnk->pert <= 0.0) {
481           /* Initialize the perturbation */
482           bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm));
483           if (bnk->is_gltr) {
484             ierr = KSPCGGLTRGetMinEig(tao->ksp, &e_min);CHKERRQ(ierr);
485             bnk->pert = PetscMax(bnk->pert, -e_min);
486           }
487         } else {
488           /* Increase the perturbation */
489           bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm));
490         }
491         break;
492 
493       default:
494         /* Newton step computation is good; decrease perturbation */
495         bnk->pert = PetscMin(bnk->psfac * bnk->pert, bnk->pmsfac * bnk->gnorm);
496         if (bnk->pert < bnk->pmin) {
497           bnk->pert = 0.0;
498         }
499         break;
500       }
501 
502       ++bnk->newt;
503       *stepType = BNK_NEWTON;
504     }
505   } else {
506     ++bnk->newt;
507     bnk->pert = 0.0;
508     *stepType = BNK_NEWTON;
509   }
510   PetscFunctionReturn(0);
511 }
512 
513 /*------------------------------------------------------------*/
514 
515 /* Routine for performing a bound-projected More-Thuente line search.
516 
517   Includes fallbacks to BFGS, scaled gradient, and unscaled gradient steps if the
518   Newton step does not produce a valid step length.
519 */
520 
521 PetscErrorCode TaoBNKPerformLineSearch(Tao tao, PetscInt stepType, PetscReal *steplen, TaoLineSearchConvergedReason *reason)
522 {
523   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
524   PetscErrorCode ierr;
525   TaoLineSearchConvergedReason ls_reason;
526 
527   Vec            active_step;
528   PetscReal      e_min, gdx, delta;
529   PetscInt       bfgsUpdates;
530 
531   PetscFunctionBegin;
532   /* Perform the linesearch */
533   ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, bnk->unprojected_gradient, tao->stepdirection, steplen, &ls_reason);CHKERRQ(ierr);
534   ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr);
535   ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
536 
537   while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != BNK_GRADIENT) {
538     /* Linesearch failed, revert solution */
539     bnk->f = bnk->fold;
540     ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr);
541     ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr);
542     ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr);
543 
544     switch(stepType) {
545     case BNK_NEWTON:
546       /* Failed to obtain acceptable iterate with Newton step
547          Update the perturbation for next time */
548       --bnk->newt;
549       if (bnk->pert <= 0.0) {
550         /* Initialize the perturbation */
551         bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm));
552         if (bnk->is_gltr) {
553           ierr = KSPCGGLTRGetMinEig(tao->ksp,&e_min);CHKERRQ(ierr);
554           bnk->pert = PetscMax(bnk->pert, -e_min);
555         }
556       } else {
557         /* Increase the perturbation */
558         bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm));
559       }
560 
561       if (BNK_PC_BFGS != bnk->pc_type) {
562         /* We don't have the bfgs matrix around and being updated
563            Must use gradient direction in this case */
564         ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
565         ++bnk->grad;
566         stepType = BNK_GRADIENT;
567       } else {
568         /* Attempt to use the BFGS direction */
569         ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr);
570         /* Check for success (descent direction)
571            NOTE: Negative gdx means not a descent direction because the step here is missing a negative sign. */
572         ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr);
573         if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
574           /* BFGS direction is not descent or direction produced not a number
575              We can assert bfgsUpdates > 1 in this case
576              Use steepest descent direction (scaled) */
577           if (bnk->f != 0.0) {
578             delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm);
579           } else {
580             delta = 2.0 / (bnk->gnorm*bnk->gnorm);
581           }
582           ierr = MatLMVMSetDelta(bnk->M, delta);CHKERRQ(ierr);
583           ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr);
584           ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr);
585           ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr);
586 
587           bfgsUpdates = 1;
588           ++bnk->sgrad;
589           stepType = BNK_SCALED_GRADIENT;
590         } else {
591           ierr = MatLMVMGetUpdates(bnk->M, &bfgsUpdates);CHKERRQ(ierr);
592           if (1 == bfgsUpdates) {
593             /* The first BFGS direction is always the scaled gradient */
594             ++bnk->sgrad;
595             stepType = BNK_SCALED_GRADIENT;
596           } else {
597             ++bnk->bfgs;
598             stepType = BNK_BFGS;
599           }
600         }
601       }
602       break;
603 
604     case BNK_BFGS:
605       /* Can only enter if pc_type == BNK_PC_BFGS
606          Failed to obtain acceptable iterate with BFGS step
607          Attempt to use the scaled gradient direction */
608       if (bnk->f != 0.0) {
609         delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm);
610       } else {
611         delta = 2.0 / (bnk->gnorm*bnk->gnorm);
612       }
613       ierr = MatLMVMSetDelta(bnk->M, delta);CHKERRQ(ierr);
614       ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr);
615       ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr);
616       ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr);
617 
618       bfgsUpdates = 1;
619       ++bnk->sgrad;
620       stepType = BNK_SCALED_GRADIENT;
621       break;
622 
623     case BNK_SCALED_GRADIENT:
624       /* Can only enter if pc_type == BNK_PC_BFGS
625          The scaled gradient step did not produce a new iterate;
626          attemp to use the gradient direction.
627          Need to make sure we are not using a different diagonal scaling */
628       ierr = MatLMVMSetScale(bnk->M,0);CHKERRQ(ierr);
629       ierr = MatLMVMSetDelta(bnk->M,1.0);CHKERRQ(ierr);
630       ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr);
631       ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr);
632       ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr);
633 
634       bfgsUpdates = 1;
635       ++bnk->grad;
636       stepType = BNK_GRADIENT;
637       break;
638     }
639     /* Make sure the safeguarded fall-back step is zero for actively bounded variables */
640     ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
641     ierr = VecGetSubVector(tao->stepdirection, bnk->active_idx, &active_step);CHKERRQ(ierr);
642     ierr = VecSet(active_step, 0.0);CHKERRQ(ierr);
643     ierr = VecRestoreSubVector(tao->stepdirection, bnk->active_idx, &active_step);CHKERRQ(ierr);
644 
645     /* Perform one last line search with the fall-back step */
646     ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, bnk->unprojected_gradient, tao->stepdirection, steplen, &ls_reason);CHKERRQ(ierr);
647     ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr);
648     ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
649   }
650   *reason = ls_reason;
651   PetscFunctionReturn(0);
652 }
653 
654 /*------------------------------------------------------------*/
655 
656 /* Routine for updating the trust radius.
657 
658   Function features three different update methods:
659   1) Line-search step length based
660   2) Predicted decrease on the CG quadratic model
661   3) Interpolation
662 */
663 
664 PetscErrorCode TaoBNKUpdateTrustRadius(Tao tao, PetscReal prered, PetscReal actred, PetscInt stepType, PetscBool *accept)
665 {
666   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
667   PetscErrorCode ierr;
668 
669   PetscReal      step, kappa;
670   PetscReal      gdx, tau_1, tau_2, tau_min, tau_max;
671 
672   PetscFunctionBegin;
673   /* Update trust region radius */
674   *accept = PETSC_FALSE;
675   switch(bnk->update_type) {
676   case BNK_UPDATE_STEP:
677     *accept = PETSC_TRUE; /* always accept here because line search succeeded */
678     if (stepType == BNK_NEWTON) {
679       ierr = TaoLineSearchGetStepLength(tao->linesearch, &step);CHKERRQ(ierr);
680       if (step < bnk->nu1) {
681         /* Very bad step taken; reduce radius */
682         tao->trust = bnk->omega1 * PetscMin(bnk->dnorm, tao->trust);
683       } else if (step < bnk->nu2) {
684         /* Reasonably bad step taken; reduce radius */
685         tao->trust = bnk->omega2 * PetscMin(bnk->dnorm, tao->trust);
686       } else if (step < bnk->nu3) {
687         /*  Reasonable step was taken; leave radius alone */
688         if (bnk->omega3 < 1.0) {
689           tao->trust = bnk->omega3 * PetscMin(bnk->dnorm, tao->trust);
690         } else if (bnk->omega3 > 1.0) {
691           tao->trust = PetscMax(bnk->omega3 * bnk->dnorm, tao->trust);
692         }
693       } else if (step < bnk->nu4) {
694         /*  Full step taken; increase the radius */
695         tao->trust = PetscMax(bnk->omega4 * bnk->dnorm, tao->trust);
696       } else {
697         /*  More than full step taken; increase the radius */
698         tao->trust = PetscMax(bnk->omega5 * bnk->dnorm, tao->trust);
699       }
700     } else {
701       /*  Newton step was not good; reduce the radius */
702       tao->trust = bnk->omega1 * PetscMin(bnk->dnorm, tao->trust);
703     }
704     break;
705 
706   case BNK_UPDATE_REDUCTION:
707     if (stepType == BNK_NEWTON) {
708       if (prered < 0.0) {
709         /* The predicted reduction has the wrong sign.  This cannot
710            happen in infinite precision arithmetic.  Step should
711            be rejected! */
712         tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm);
713       }
714       else {
715         if (PetscIsInfOrNanReal(actred)) {
716           tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm);
717         } else {
718           if ((PetscAbsScalar(actred) <= bnk->epsilon) &&
719               (PetscAbsScalar(prered) <= bnk->epsilon)) {
720             kappa = 1.0;
721           }
722           else {
723             kappa = actred / prered;
724           }
725 
726           /* Accept or reject the step and update radius */
727           if (kappa < bnk->eta1) {
728             /* Reject the step */
729             tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm);
730           }
731           else {
732             /* Accept the step */
733             *accept = PETSC_TRUE;
734             /* Update the trust region radius only if the computed step is at the trust radius boundary */
735             if (bnk->dnorm == tao->trust) {
736               if (kappa < bnk->eta2) {
737                 /* Marginal bad step */
738                 tao->trust = bnk->alpha2 * tao->trust;
739               }
740               else if (kappa < bnk->eta3) {
741                 /* Reasonable step */
742                 tao->trust = bnk->alpha3 * tao->trust;
743               }
744               else if (kappa < bnk->eta4) {
745                 /* Good step */
746                 tao->trust = bnk->alpha4 * tao->trust;
747               }
748               else {
749                 /* Very good step */
750                 tao->trust = bnk->alpha5 * tao->trust;
751               }
752             }
753           }
754         }
755       }
756     } else {
757       /*  Newton step was not good; reduce the radius */
758       tao->trust = bnk->alpha1 * PetscMin(bnk->dnorm, tao->trust);
759     }
760     break;
761 
762   default:
763     if (stepType == BNK_NEWTON) {
764       if (prered < 0.0) {
765         /*  The predicted reduction has the wrong sign.  This cannot */
766         /*  happen in infinite precision arithmetic.  Step should */
767         /*  be rejected! */
768         tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm);
769       } else {
770         if (PetscIsInfOrNanReal(actred)) {
771           tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm);
772         } else {
773           if ((PetscAbsScalar(actred) <= bnk->epsilon) && (PetscAbsScalar(prered) <= bnk->epsilon)) {
774             kappa = 1.0;
775           } else {
776             kappa = actred / prered;
777           }
778 
779           ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr);
780           tau_1 = bnk->theta * gdx / (bnk->theta * gdx - (1.0 - bnk->theta) * prered + actred);
781           tau_2 = bnk->theta * gdx / (bnk->theta * gdx + (1.0 + bnk->theta) * prered - actred);
782           tau_min = PetscMin(tau_1, tau_2);
783           tau_max = PetscMax(tau_1, tau_2);
784 
785           if (kappa >= 1.0 - bnk->mu1) {
786             /*  Great agreement */
787             *accept = PETSC_TRUE;
788             if (tau_max < 1.0) {
789               tao->trust = PetscMax(tao->trust, bnk->gamma3 * bnk->dnorm);
790             } else if (tau_max > bnk->gamma4) {
791               tao->trust = PetscMax(tao->trust, bnk->gamma4 * bnk->dnorm);
792             } else {
793               tao->trust = PetscMax(tao->trust, tau_max * bnk->dnorm);
794             }
795           } else if (kappa >= 1.0 - bnk->mu2) {
796             /*  Good agreement */
797             *accept = PETSC_TRUE;
798             if (tau_max < bnk->gamma2) {
799               tao->trust = bnk->gamma2 * PetscMin(tao->trust, bnk->dnorm);
800             } else if (tau_max > bnk->gamma3) {
801               tao->trust = PetscMax(tao->trust, bnk->gamma3 * bnk->dnorm);
802             } else if (tau_max < 1.0) {
803               tao->trust = tau_max * PetscMin(tao->trust, bnk->dnorm);
804             } else {
805               tao->trust = PetscMax(tao->trust, tau_max * bnk->dnorm);
806             }
807           } else {
808             /*  Not good agreement */
809             if (tau_min > 1.0) {
810               tao->trust = bnk->gamma2 * PetscMin(tao->trust, bnk->dnorm);
811             } else if (tau_max < bnk->gamma1) {
812               tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm);
813             } else if ((tau_min < bnk->gamma1) && (tau_max >= 1.0)) {
814               tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm);
815             } else if ((tau_1 >= bnk->gamma1) && (tau_1 < 1.0) && ((tau_2 < bnk->gamma1) || (tau_2 >= 1.0))) {
816               tao->trust = tau_1 * PetscMin(tao->trust, bnk->dnorm);
817             } else if ((tau_2 >= bnk->gamma1) && (tau_2 < 1.0) && ((tau_1 < bnk->gamma1) || (tau_2 >= 1.0))) {
818               tao->trust = tau_2 * PetscMin(tao->trust, bnk->dnorm);
819             } else {
820               tao->trust = tau_max * PetscMin(tao->trust, bnk->dnorm);
821             }
822           }
823         }
824       }
825     } else {
826       /*  Newton step was not good; reduce the radius */
827       tao->trust = bnk->gamma1 * PetscMin(bnk->dnorm, tao->trust);
828     }
829   }
830   /* Make sure the radius does not violate min and max settings */
831   tao->trust = PetscMin(tao->trust, bnk->max_radius);
832   tao->trust = PetscMax(tao->trust, bnk->min_radius);
833   PetscFunctionReturn(0);
834 }
835 
836 /* ---------------------------------------------------------- */
837 
838 static PetscErrorCode TaoSetUp_BNK(Tao tao)
839 {
840   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
841   PetscErrorCode ierr;
842 
843   PetscFunctionBegin;
844   if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);}
845   if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr);}
846   if (!bnk->W) {ierr = VecDuplicate(tao->solution,&bnk->W);CHKERRQ(ierr);}
847   if (!bnk->Xold) {ierr = VecDuplicate(tao->solution,&bnk->Xold);CHKERRQ(ierr);}
848   if (!bnk->Gold) {ierr = VecDuplicate(tao->solution,&bnk->Gold);CHKERRQ(ierr);}
849   if (!bnk->Xwork) {ierr = VecDuplicate(tao->solution,&bnk->Xwork);CHKERRQ(ierr);}
850   if (!bnk->Gwork) {ierr = VecDuplicate(tao->solution,&bnk->Gwork);CHKERRQ(ierr);}
851   if (!bnk->unprojected_gradient) {ierr = VecDuplicate(tao->solution,&bnk->unprojected_gradient);CHKERRQ(ierr);}
852   if (!bnk->unprojected_gradient_old) {ierr = VecDuplicate(tao->solution,&bnk->unprojected_gradient_old);CHKERRQ(ierr);}
853   bnk->Diag = 0;
854   bnk->M = 0;
855   bnk->H_inactive = 0;
856   bnk->Hpre_inactive = 0;
857   PetscFunctionReturn(0);
858 }
859 
860 /*------------------------------------------------------------*/
861 
862 static PetscErrorCode TaoDestroy_BNK(Tao tao)
863 {
864   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
865   PetscErrorCode ierr;
866 
867   PetscFunctionBegin;
868   if (tao->setupcalled) {
869     ierr = VecDestroy(&bnk->W);CHKERRQ(ierr);
870     ierr = VecDestroy(&bnk->Xold);CHKERRQ(ierr);
871     ierr = VecDestroy(&bnk->Gold);CHKERRQ(ierr);
872     ierr = VecDestroy(&bnk->Xwork);CHKERRQ(ierr);
873     ierr = VecDestroy(&bnk->Gwork);CHKERRQ(ierr);
874     ierr = VecDestroy(&bnk->unprojected_gradient);CHKERRQ(ierr);
875     ierr = VecDestroy(&bnk->unprojected_gradient_old);CHKERRQ(ierr);
876   }
877   ierr = VecDestroy(&bnk->Diag);CHKERRQ(ierr);
878   ierr = MatDestroy(&bnk->M);CHKERRQ(ierr);
879   ierr = PetscFree(tao->data);CHKERRQ(ierr);
880   PetscFunctionReturn(0);
881 }
882 
883 /*------------------------------------------------------------*/
884 
885 static PetscErrorCode TaoSetFromOptions_BNK(PetscOptionItems *PetscOptionsObject,Tao tao)
886 {
887   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
888   PetscErrorCode ierr;
889 
890   PetscFunctionBegin;
891   ierr = PetscOptionsHead(PetscOptionsObject,"Newton line search method for unconstrained optimization");CHKERRQ(ierr);
892   ierr = PetscOptionsEList("-tao_bnk_pc_type", "pc type", "", BNK_PC, BNK_PC_TYPES, BNK_PC[bnk->pc_type], &bnk->pc_type, 0);CHKERRQ(ierr);
893   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);
894   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);
895   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);
896   ierr = PetscOptionsReal("-tao_bnk_sval", "perturbation starting value", "", bnk->sval, &bnk->sval,NULL);CHKERRQ(ierr);
897   ierr = PetscOptionsReal("-tao_bnk_imin", "minimum initial perturbation", "", bnk->imin, &bnk->imin,NULL);CHKERRQ(ierr);
898   ierr = PetscOptionsReal("-tao_bnk_imax", "maximum initial perturbation", "", bnk->imax, &bnk->imax,NULL);CHKERRQ(ierr);
899   ierr = PetscOptionsReal("-tao_bnk_imfac", "initial merit factor", "", bnk->imfac, &bnk->imfac,NULL);CHKERRQ(ierr);
900   ierr = PetscOptionsReal("-tao_bnk_pmin", "minimum perturbation", "", bnk->pmin, &bnk->pmin,NULL);CHKERRQ(ierr);
901   ierr = PetscOptionsReal("-tao_bnk_pmax", "maximum perturbation", "", bnk->pmax, &bnk->pmax,NULL);CHKERRQ(ierr);
902   ierr = PetscOptionsReal("-tao_bnk_pgfac", "growth factor", "", bnk->pgfac, &bnk->pgfac,NULL);CHKERRQ(ierr);
903   ierr = PetscOptionsReal("-tao_bnk_psfac", "shrink factor", "", bnk->psfac, &bnk->psfac,NULL);CHKERRQ(ierr);
904   ierr = PetscOptionsReal("-tao_bnk_pmgfac", "merit growth factor", "", bnk->pmgfac, &bnk->pmgfac,NULL);CHKERRQ(ierr);
905   ierr = PetscOptionsReal("-tao_bnk_pmsfac", "merit shrink factor", "", bnk->pmsfac, &bnk->pmsfac,NULL);CHKERRQ(ierr);
906   ierr = PetscOptionsReal("-tao_bnk_eta1", "poor steplength; reduce radius", "", bnk->eta1, &bnk->eta1,NULL);CHKERRQ(ierr);
907   ierr = PetscOptionsReal("-tao_bnk_eta2", "reasonable steplength; leave radius alone", "", bnk->eta2, &bnk->eta2,NULL);CHKERRQ(ierr);
908   ierr = PetscOptionsReal("-tao_bnk_eta3", "good steplength; increase radius", "", bnk->eta3, &bnk->eta3,NULL);CHKERRQ(ierr);
909   ierr = PetscOptionsReal("-tao_bnk_eta4", "excellent steplength; greatly increase radius", "", bnk->eta4, &bnk->eta4,NULL);CHKERRQ(ierr);
910   ierr = PetscOptionsReal("-tao_bnk_alpha1", "", "", bnk->alpha1, &bnk->alpha1,NULL);CHKERRQ(ierr);
911   ierr = PetscOptionsReal("-tao_bnk_alpha2", "", "", bnk->alpha2, &bnk->alpha2,NULL);CHKERRQ(ierr);
912   ierr = PetscOptionsReal("-tao_bnk_alpha3", "", "", bnk->alpha3, &bnk->alpha3,NULL);CHKERRQ(ierr);
913   ierr = PetscOptionsReal("-tao_bnk_alpha4", "", "", bnk->alpha4, &bnk->alpha4,NULL);CHKERRQ(ierr);
914   ierr = PetscOptionsReal("-tao_bnk_alpha5", "", "", bnk->alpha5, &bnk->alpha5,NULL);CHKERRQ(ierr);
915   ierr = PetscOptionsReal("-tao_bnk_nu1", "poor steplength; reduce radius", "", bnk->nu1, &bnk->nu1,NULL);CHKERRQ(ierr);
916   ierr = PetscOptionsReal("-tao_bnk_nu2", "reasonable steplength; leave radius alone", "", bnk->nu2, &bnk->nu2,NULL);CHKERRQ(ierr);
917   ierr = PetscOptionsReal("-tao_bnk_nu3", "good steplength; increase radius", "", bnk->nu3, &bnk->nu3,NULL);CHKERRQ(ierr);
918   ierr = PetscOptionsReal("-tao_bnk_nu4", "excellent steplength; greatly increase radius", "", bnk->nu4, &bnk->nu4,NULL);CHKERRQ(ierr);
919   ierr = PetscOptionsReal("-tao_bnk_omega1", "", "", bnk->omega1, &bnk->omega1,NULL);CHKERRQ(ierr);
920   ierr = PetscOptionsReal("-tao_bnk_omega2", "", "", bnk->omega2, &bnk->omega2,NULL);CHKERRQ(ierr);
921   ierr = PetscOptionsReal("-tao_bnk_omega3", "", "", bnk->omega3, &bnk->omega3,NULL);CHKERRQ(ierr);
922   ierr = PetscOptionsReal("-tao_bnk_omega4", "", "", bnk->omega4, &bnk->omega4,NULL);CHKERRQ(ierr);
923   ierr = PetscOptionsReal("-tao_bnk_omega5", "", "", bnk->omega5, &bnk->omega5,NULL);CHKERRQ(ierr);
924   ierr = PetscOptionsReal("-tao_bnk_mu1_i", "", "", bnk->mu1_i, &bnk->mu1_i,NULL);CHKERRQ(ierr);
925   ierr = PetscOptionsReal("-tao_bnk_mu2_i", "", "", bnk->mu2_i, &bnk->mu2_i,NULL);CHKERRQ(ierr);
926   ierr = PetscOptionsReal("-tao_bnk_gamma1_i", "", "", bnk->gamma1_i, &bnk->gamma1_i,NULL);CHKERRQ(ierr);
927   ierr = PetscOptionsReal("-tao_bnk_gamma2_i", "", "", bnk->gamma2_i, &bnk->gamma2_i,NULL);CHKERRQ(ierr);
928   ierr = PetscOptionsReal("-tao_bnk_gamma3_i", "", "", bnk->gamma3_i, &bnk->gamma3_i,NULL);CHKERRQ(ierr);
929   ierr = PetscOptionsReal("-tao_bnk_gamma4_i", "", "", bnk->gamma4_i, &bnk->gamma4_i,NULL);CHKERRQ(ierr);
930   ierr = PetscOptionsReal("-tao_bnk_theta_i", "", "", bnk->theta_i, &bnk->theta_i,NULL);CHKERRQ(ierr);
931   ierr = PetscOptionsReal("-tao_bnk_mu1", "", "", bnk->mu1, &bnk->mu1,NULL);CHKERRQ(ierr);
932   ierr = PetscOptionsReal("-tao_bnk_mu2", "", "", bnk->mu2, &bnk->mu2,NULL);CHKERRQ(ierr);
933   ierr = PetscOptionsReal("-tao_bnk_gamma1", "", "", bnk->gamma1, &bnk->gamma1,NULL);CHKERRQ(ierr);
934   ierr = PetscOptionsReal("-tao_bnk_gamma2", "", "", bnk->gamma2, &bnk->gamma2,NULL);CHKERRQ(ierr);
935   ierr = PetscOptionsReal("-tao_bnk_gamma3", "", "", bnk->gamma3, &bnk->gamma3,NULL);CHKERRQ(ierr);
936   ierr = PetscOptionsReal("-tao_bnk_gamma4", "", "", bnk->gamma4, &bnk->gamma4,NULL);CHKERRQ(ierr);
937   ierr = PetscOptionsReal("-tao_bnk_theta", "", "", bnk->theta, &bnk->theta,NULL);CHKERRQ(ierr);
938   ierr = PetscOptionsReal("-tao_bnk_min_radius", "lower bound on initial radius", "", bnk->min_radius, &bnk->min_radius,NULL);CHKERRQ(ierr);
939   ierr = PetscOptionsReal("-tao_bnk_max_radius", "upper bound on radius", "", bnk->max_radius, &bnk->max_radius,NULL);CHKERRQ(ierr);
940   ierr = PetscOptionsReal("-tao_bnk_epsilon", "tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL);CHKERRQ(ierr);
941   ierr = PetscOptionsTail();CHKERRQ(ierr);
942   ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
943   ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
944   PetscFunctionReturn(0);
945 }
946 
947 /*------------------------------------------------------------*/
948 
949 static PetscErrorCode TaoView_BNK(Tao tao, PetscViewer viewer)
950 {
951   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
952   PetscInt       nrejects;
953   PetscBool      isascii;
954   PetscErrorCode ierr;
955 
956   PetscFunctionBegin;
957   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
958   if (isascii) {
959     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
960     if (BNK_PC_BFGS == bnk->pc_type && bnk->M) {
961       ierr = MatLMVMGetRejects(bnk->M,&nrejects);CHKERRQ(ierr);
962       ierr = PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n",nrejects);CHKERRQ(ierr);
963     }
964     ierr = PetscViewerASCIIPrintf(viewer, "Newton steps: %D\n", bnk->newt);CHKERRQ(ierr);
965     ierr = PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", bnk->bfgs);CHKERRQ(ierr);
966     ierr = PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", bnk->sgrad);CHKERRQ(ierr);
967     ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", bnk->grad);CHKERRQ(ierr);
968     ierr = PetscViewerASCIIPrintf(viewer, "KSP termination reasons:\n");CHKERRQ(ierr);
969     ierr = PetscViewerASCIIPrintf(viewer, "  atol: %D\n", bnk->ksp_atol);CHKERRQ(ierr);
970     ierr = PetscViewerASCIIPrintf(viewer, "  rtol: %D\n", bnk->ksp_rtol);CHKERRQ(ierr);
971     ierr = PetscViewerASCIIPrintf(viewer, "  ctol: %D\n", bnk->ksp_ctol);CHKERRQ(ierr);
972     ierr = PetscViewerASCIIPrintf(viewer, "  negc: %D\n", bnk->ksp_negc);CHKERRQ(ierr);
973     ierr = PetscViewerASCIIPrintf(viewer, "  dtol: %D\n", bnk->ksp_dtol);CHKERRQ(ierr);
974     ierr = PetscViewerASCIIPrintf(viewer, "  iter: %D\n", bnk->ksp_iter);CHKERRQ(ierr);
975     ierr = PetscViewerASCIIPrintf(viewer, "  othr: %D\n", bnk->ksp_othr);CHKERRQ(ierr);
976     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
977   }
978   PetscFunctionReturn(0);
979 }
980 
981 /* ---------------------------------------------------------- */
982 
983 /*MC
984   TAOBNK - Shared base-type for Bounded Newton-Krylov type algorithms.
985   At each iteration, the BNK methods solve the symmetric
986   system of equations to obtain the step diretion dk:
987               Hk dk = -gk
988   for free variables only. The step can be globalized either through
989   trust-region methods, or a line search, or a heuristic mixture of both.
990 
991     Options Database Keys:
992 + -tao_bnk_pc_type - "none","ahess","bfgs","petsc"
993 . -tao_bnk_bfgs_scale_type - "ahess","phess","bfgs"
994 . -tao_bnk_init_type - "constant","direction","interpolation"
995 . -tao_bnk_update_type - "step","direction","interpolation"
996 . -tao_bnk_sval - perturbation starting value
997 . -tao_bnk_imin - minimum initial perturbation
998 . -tao_bnk_imax - maximum initial perturbation
999 . -tao_bnk_pmin - minimum perturbation
1000 . -tao_bnk_pmax - maximum perturbation
1001 . -tao_bnk_pgfac - growth factor
1002 . -tao_bnk_psfac - shrink factor
1003 . -tao_bnk_imfac - initial merit factor
1004 . -tao_bnk_pmgfac - merit growth factor
1005 . -tao_bnk_pmsfac - merit shrink factor
1006 . -tao_bnk_eta1 - poor steplength; reduce radius
1007 . -tao_bnk_eta2 - reasonable steplength; leave radius
1008 . -tao_bnk_eta3 - good steplength; increase readius
1009 . -tao_bnk_eta4 - excellent steplength; greatly increase radius
1010 . -tao_bnk_alpha1 - alpha1 reduction
1011 . -tao_bnk_alpha2 - alpha2 reduction
1012 . -tao_bnk_alpha3 - alpha3 reduction
1013 . -tao_bnk_alpha4 - alpha4 reduction
1014 . -tao_bnk_alpha - alpha5 reduction
1015 . -tao_bnk_mu1 - mu1 interpolation update
1016 . -tao_bnk_mu2 - mu2 interpolation update
1017 . -tao_bnk_gamma1 - gamma1 interpolation update
1018 . -tao_bnk_gamma2 - gamma2 interpolation update
1019 . -tao_bnk_gamma3 - gamma3 interpolation update
1020 . -tao_bnk_gamma4 - gamma4 interpolation update
1021 . -tao_bnk_theta - theta interpolation update
1022 . -tao_bnk_omega1 - omega1 step update
1023 . -tao_bnk_omega2 - omega2 step update
1024 . -tao_bnk_omega3 - omega3 step update
1025 . -tao_bnk_omega4 - omega4 step update
1026 . -tao_bnk_omega5 - omega5 step update
1027 . -tao_bnk_mu1_i -  mu1 interpolation init factor
1028 . -tao_bnk_mu2_i -  mu2 interpolation init factor
1029 . -tao_bnk_gamma1_i -  gamma1 interpolation init factor
1030 . -tao_bnk_gamma2_i -  gamma2 interpolation init factor
1031 . -tao_bnk_gamma3_i -  gamma3 interpolation init factor
1032 . -tao_bnk_gamma4_i -  gamma4 interpolation init factor
1033 - -tao_bnk_theta_i -  theta interpolation init factor
1034 
1035   Level: beginner
1036 M*/
1037 
1038 PetscErrorCode TaoCreate_BNK(Tao tao)
1039 {
1040   TAO_BNK        *bnk;
1041   const char     *morethuente_type = TAOLINESEARCHMT;
1042   PetscErrorCode ierr;
1043 
1044   PetscFunctionBegin;
1045   ierr = PetscNewLog(tao,&bnk);CHKERRQ(ierr);
1046 
1047   tao->ops->setup = TaoSetUp_BNK;
1048   tao->ops->view = TaoView_BNK;
1049   tao->ops->setfromoptions = TaoSetFromOptions_BNK;
1050   tao->ops->destroy = TaoDestroy_BNK;
1051 
1052   /*  Override default settings (unless already changed) */
1053   if (!tao->max_it_changed) tao->max_it = 50;
1054   if (!tao->trust0_changed) tao->trust0 = 100.0;
1055 
1056   tao->data = (void*)bnk;
1057 
1058   /*  Hessian shifting parameters */
1059   bnk->sval   = 0.0;
1060   bnk->imin   = 1.0e-4;
1061   bnk->imax   = 1.0e+2;
1062   bnk->imfac  = 1.0e-1;
1063 
1064   bnk->pmin   = 1.0e-12;
1065   bnk->pmax   = 1.0e+2;
1066   bnk->pgfac  = 1.0e+1;
1067   bnk->psfac  = 4.0e-1;
1068   bnk->pmgfac = 1.0e-1;
1069   bnk->pmsfac = 1.0e-1;
1070 
1071   /*  Default values for trust-region radius update based on steplength */
1072   bnk->nu1 = 0.25;
1073   bnk->nu2 = 0.50;
1074   bnk->nu3 = 1.00;
1075   bnk->nu4 = 1.25;
1076 
1077   bnk->omega1 = 0.25;
1078   bnk->omega2 = 0.50;
1079   bnk->omega3 = 1.00;
1080   bnk->omega4 = 2.00;
1081   bnk->omega5 = 4.00;
1082 
1083   /*  Default values for trust-region radius update based on reduction */
1084   bnk->eta1 = 1.0e-4;
1085   bnk->eta2 = 0.25;
1086   bnk->eta3 = 0.50;
1087   bnk->eta4 = 0.90;
1088 
1089   bnk->alpha1 = 0.25;
1090   bnk->alpha2 = 0.50;
1091   bnk->alpha3 = 1.00;
1092   bnk->alpha4 = 2.00;
1093   bnk->alpha5 = 4.00;
1094 
1095   /*  Default values for trust-region radius update based on interpolation */
1096   bnk->mu1 = 0.10;
1097   bnk->mu2 = 0.50;
1098 
1099   bnk->gamma1 = 0.25;
1100   bnk->gamma2 = 0.50;
1101   bnk->gamma3 = 2.00;
1102   bnk->gamma4 = 4.00;
1103 
1104   bnk->theta = 0.05;
1105 
1106   /*  Default values for trust region initialization based on interpolation */
1107   bnk->mu1_i = 0.35;
1108   bnk->mu2_i = 0.50;
1109 
1110   bnk->gamma1_i = 0.0625;
1111   bnk->gamma2_i = 0.5;
1112   bnk->gamma3_i = 2.0;
1113   bnk->gamma4_i = 5.0;
1114 
1115   bnk->theta_i = 0.25;
1116 
1117   /*  Remaining parameters */
1118   bnk->min_radius = 1.0e-10;
1119   bnk->max_radius = 1.0e10;
1120   bnk->epsilon = 1.0e-6;
1121 
1122   bnk->pc_type         = BNK_PC_BFGS;
1123   bnk->bfgs_scale_type = BFGS_SCALE_PHESS;
1124   bnk->init_type       = BNK_INIT_INTERPOLATION;
1125   bnk->update_type     = BNK_UPDATE_INTERPOLATION;
1126 
1127   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr);
1128   ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr);
1129   ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr);
1130   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr);
1131   ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
1132 
1133   /*  Set linear solver to default for symmetric matrices */
1134   ierr = KSPCreate(((PetscObject)tao)->comm,&tao->ksp);CHKERRQ(ierr);
1135   ierr = PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);CHKERRQ(ierr);
1136   ierr = KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);CHKERRQ(ierr);
1137   ierr = KSPSetType(tao->ksp,KSPCGSTCG);CHKERRQ(ierr);
1138   PetscFunctionReturn(0);
1139 }
1140