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