xref: /petsc/src/tao/unconstrained/impls/ntl/ntl.c (revision 4fc747eaadbeca11629f314a99edccbc2ed7b3d3) 
1 #include <../src/tao/matrix/lmvmmat.h>
2 #include <../src/tao/unconstrained/impls/ntl/ntl.h>
3 
4 #include <petscksp.h>
5 
6 #define NTL_KSP_NASH    0
7 #define NTL_KSP_STCG    1
8 #define NTL_KSP_GLTR    2
9 #define NTL_KSP_TYPES   3
10 
11 #define NTL_PC_NONE     0
12 #define NTL_PC_AHESS    1
13 #define NTL_PC_BFGS     2
14 #define NTL_PC_PETSC    3
15 #define NTL_PC_TYPES    4
16 
17 #define BFGS_SCALE_AHESS        0
18 #define BFGS_SCALE_BFGS         1
19 #define BFGS_SCALE_TYPES        2
20 
21 #define NTL_INIT_CONSTANT         0
22 #define NTL_INIT_DIRECTION        1
23 #define NTL_INIT_INTERPOLATION    2
24 #define NTL_INIT_TYPES            3
25 
26 #define NTL_UPDATE_REDUCTION      0
27 #define NTL_UPDATE_INTERPOLATION  1
28 #define NTL_UPDATE_TYPES          2
29 
30 static const char *NTL_KSP[64] = {"nash", "stcg", "gltr"};
31 
32 static const char *NTL_PC[64] = {"none", "ahess", "bfgs", "petsc"};
33 
34 static const char *BFGS_SCALE[64] = {"ahess", "bfgs"};
35 
36 static const char *NTL_INIT[64] = {"constant", "direction", "interpolation"};
37 
38 static const char *NTL_UPDATE[64] = {"reduction", "interpolation"};
39 
40 /* Routine for BFGS preconditioner */
41 
42 #undef __FUNCT__
43 #define __FUNCT__ "MatLMVMSolveShell"
44 static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x)
45 {
46   PetscErrorCode ierr;
47   Mat            M;
48 
49   PetscFunctionBegin;
50   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
51   PetscValidHeaderSpecific(b,VEC_CLASSID,2);
52   PetscValidHeaderSpecific(x,VEC_CLASSID,3);
53   ierr = PCShellGetContext(pc,(void**)&M);CHKERRQ(ierr);
54   ierr = MatLMVMSolve(M, b, x);CHKERRQ(ierr);
55   PetscFunctionReturn(0);
56 }
57 
58 /* Implements Newton's Method with a trust-region, line-search approach for
59    solving unconstrained minimization problems.  A More'-Thuente line search
60    is used to guarantee that the bfgs preconditioner remains positive
61    definite. */
62 
63 #define NTL_NEWTON              0
64 #define NTL_BFGS                1
65 #define NTL_SCALED_GRADIENT     2
66 #define NTL_GRADIENT            3
67 
68 #undef __FUNCT__
69 #define __FUNCT__ "TaoSolve_NTL"
70 static PetscErrorCode TaoSolve_NTL(Tao tao)
71 {
72   TAO_NTL                      *tl = (TAO_NTL *)tao->data;
73   PC                           pc;
74   KSPConvergedReason           ksp_reason;
75   TaoConvergedReason           reason;
76   TaoLineSearchConvergedReason ls_reason;
77 
78   PetscReal                    fmin, ftrial, prered, actred, kappa, sigma;
79   PetscReal                    tau, tau_1, tau_2, tau_max, tau_min, max_radius;
80   PetscReal                    f, fold, gdx, gnorm;
81   PetscReal                    step = 1.0;
82 
83   PetscReal                    delta;
84   PetscReal                    norm_d = 0.0;
85   PetscErrorCode               ierr;
86   PetscInt                     stepType;
87   PetscInt                     its;
88 
89   PetscInt                     bfgsUpdates = 0;
90   PetscInt                     needH;
91 
92   PetscInt                     i_max = 5;
93   PetscInt                     j_max = 1;
94   PetscInt                     i, j, n, N;
95 
96   PetscInt                     tr_reject;
97 
98   PetscFunctionBegin;
99   if (tao->XL || tao->XU || tao->ops->computebounds) {
100     ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n");CHKERRQ(ierr);
101   }
102 
103   /* Initialize trust-region radius */
104   tao->trust = tao->trust0;
105 
106   /* Modify the radius if it is too large or small */
107   tao->trust = PetscMax(tao->trust, tl->min_radius);
108   tao->trust = PetscMin(tao->trust, tl->max_radius);
109 
110   if (NTL_PC_BFGS == tl->pc_type && !tl->M) {
111     ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr);
112     ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);
113     ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tl->M);CHKERRQ(ierr);
114     ierr = MatLMVMAllocateVectors(tl->M,tao->solution);CHKERRQ(ierr);
115   }
116 
117   /* Check convergence criteria */
118   ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
119   ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr);
120   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
121   needH = 1;
122 
123   ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr);
124   if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
125 
126   /* Create vectors for the limited memory preconditioner */
127   if ((NTL_PC_BFGS == tl->pc_type) && (BFGS_SCALE_BFGS != tl->bfgs_scale_type)) {
128     if (!tl->Diag) {
129       ierr = VecDuplicate(tao->solution, &tl->Diag);CHKERRQ(ierr);
130     }
131   }
132 
133   /* Modify the linear solver to a conjugate gradient method */
134   switch(tl->ksp_type) {
135   case NTL_KSP_NASH:
136     ierr = KSPSetType(tao->ksp, KSPNASH);CHKERRQ(ierr);
137     ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
138     break;
139 
140   case NTL_KSP_STCG:
141     ierr = KSPSetType(tao->ksp, KSPSTCG);CHKERRQ(ierr);
142     ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
143     break;
144 
145   default:
146     ierr = KSPSetType(tao->ksp, KSPGLTR);CHKERRQ(ierr);
147     ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
148     break;
149   }
150 
151   /* Modify the preconditioner to use the bfgs approximation */
152   ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr);
153   switch(tl->pc_type) {
154   case NTL_PC_NONE:
155     ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr);
156     ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
157     break;
158 
159   case NTL_PC_AHESS:
160     ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr);
161     ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
162     ierr = PCJacobiSetUseAbs(pc,PETSC_TRUE);CHKERRQ(ierr);
163     break;
164 
165   case NTL_PC_BFGS:
166     ierr = PCSetType(pc, PCSHELL);CHKERRQ(ierr);
167     ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
168     ierr = PCShellSetName(pc, "bfgs");CHKERRQ(ierr);
169     ierr = PCShellSetContext(pc, tl->M);CHKERRQ(ierr);
170     ierr = PCShellSetApply(pc, MatLMVMSolveShell);CHKERRQ(ierr);
171     break;
172 
173   default:
174     /* Use the pc method set by pc_type */
175     break;
176   }
177 
178   /* Initialize trust-region radius.  The initialization is only performed
179      when we are using Steihaug-Toint or the Generalized Lanczos method. */
180   switch(tl->init_type) {
181   case NTL_INIT_CONSTANT:
182     /* Use the initial radius specified */
183     break;
184 
185   case NTL_INIT_INTERPOLATION:
186     /* Use the initial radius specified */
187     max_radius = 0.0;
188 
189     for (j = 0; j < j_max; ++j) {
190       fmin = f;
191       sigma = 0.0;
192 
193       if (needH) {
194         ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
195         needH = 0;
196       }
197 
198       for (i = 0; i < i_max; ++i) {
199         ierr = VecCopy(tao->solution, tl->W);CHKERRQ(ierr);
200         ierr = VecAXPY(tl->W, -tao->trust/gnorm, tao->gradient);CHKERRQ(ierr);
201 
202         ierr = TaoComputeObjective(tao, tl->W, &ftrial);CHKERRQ(ierr);
203         if (PetscIsInfOrNanReal(ftrial)) {
204           tau = tl->gamma1_i;
205         } else {
206           if (ftrial < fmin) {
207             fmin = ftrial;
208             sigma = -tao->trust / gnorm;
209           }
210 
211           ierr = MatMult(tao->hessian, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
212           ierr = VecDot(tao->gradient, tao->stepdirection, &prered);CHKERRQ(ierr);
213 
214           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
215           actred = f - ftrial;
216           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
217             kappa = 1.0;
218           } else {
219             kappa = actred / prered;
220           }
221 
222           tau_1 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred);
223           tau_2 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred);
224           tau_min = PetscMin(tau_1, tau_2);
225           tau_max = PetscMax(tau_1, tau_2);
226 
227           if (PetscAbsScalar(kappa - 1.0) <= tl->mu1_i) {
228             /* Great agreement */
229             max_radius = PetscMax(max_radius, tao->trust);
230 
231             if (tau_max < 1.0) {
232               tau = tl->gamma3_i;
233             } else if (tau_max > tl->gamma4_i) {
234               tau = tl->gamma4_i;
235             } else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) {
236               tau = tau_1;
237             } else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) {
238               tau = tau_2;
239             } else {
240               tau = tau_max;
241             }
242           } else if (PetscAbsScalar(kappa - 1.0) <= tl->mu2_i) {
243             /* Good agreement */
244             max_radius = PetscMax(max_radius, tao->trust);
245 
246             if (tau_max < tl->gamma2_i) {
247               tau = tl->gamma2_i;
248             } else if (tau_max > tl->gamma3_i) {
249               tau = tl->gamma3_i;
250             } else {
251               tau = tau_max;
252             }
253           } else {
254             /* Not good agreement */
255             if (tau_min > 1.0) {
256               tau = tl->gamma2_i;
257             } else if (tau_max < tl->gamma1_i) {
258               tau = tl->gamma1_i;
259             } else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) {
260               tau = tl->gamma1_i;
261             } else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) &&  ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) {
262               tau = tau_1;
263             } else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) &&  ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) {
264               tau = tau_2;
265             } else {
266               tau = tau_max;
267             }
268           }
269         }
270         tao->trust = tau * tao->trust;
271       }
272 
273       if (fmin < f) {
274         f = fmin;
275         ierr = VecAXPY(tao->solution, sigma, tao->gradient);CHKERRQ(ierr);
276         ierr = TaoComputeGradient(tao, tao->solution, tao->gradient);CHKERRQ(ierr);
277 
278         ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr);
279         if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
280         needH = 1;
281 
282         ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr);
283         if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
284       }
285     }
286     tao->trust = PetscMax(tao->trust, max_radius);
287 
288     /* Modify the radius if it is too large or small */
289     tao->trust = PetscMax(tao->trust, tl->min_radius);
290     tao->trust = PetscMin(tao->trust, tl->max_radius);
291     break;
292 
293   default:
294     /* Norm of the first direction will initialize radius */
295     tao->trust = 0.0;
296     break;
297   }
298 
299   /* Set initial scaling for the BFGS preconditioner
300      This step is done after computing the initial trust-region radius
301      since the function value may have decreased */
302   if (NTL_PC_BFGS == tl->pc_type) {
303     if (f != 0.0) {
304       delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
305     } else {
306       delta = 2.0 / (gnorm*gnorm);
307     }
308     ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr);
309   }
310 
311   /* Set counter for gradient/reset steps */
312   tl->ntrust = 0;
313   tl->newt = 0;
314   tl->bfgs = 0;
315   tl->sgrad = 0;
316   tl->grad = 0;
317 
318   /* Have not converged; continue with Newton method */
319   while (reason == TAO_CONTINUE_ITERATING) {
320     ++tao->niter;
321     tao->ksp_its=0;
322     /* Compute the Hessian */
323     if (needH) {
324       ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
325     }
326 
327     if (NTL_PC_BFGS == tl->pc_type) {
328       if (BFGS_SCALE_AHESS == tl->bfgs_scale_type) {
329         /* Obtain diagonal for the bfgs preconditioner */
330         ierr = MatGetDiagonal(tao->hessian, tl->Diag);CHKERRQ(ierr);
331         ierr = VecAbs(tl->Diag);CHKERRQ(ierr);
332         ierr = VecReciprocal(tl->Diag);CHKERRQ(ierr);
333         ierr = MatLMVMSetScale(tl->M, tl->Diag);CHKERRQ(ierr);
334       }
335 
336       /* Update the limited memory preconditioner */
337       ierr = MatLMVMUpdate(tl->M,tao->solution, tao->gradient);CHKERRQ(ierr);
338       ++bfgsUpdates;
339     }
340     ierr = KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);CHKERRQ(ierr);
341     /* Solve the Newton system of equations */
342     if (NTL_KSP_NASH == tl->ksp_type) {
343       ierr = KSPNASHSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr);
344       ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
345       ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
346       tao->ksp_its+=its;
347       tao->ksp_tot_its+=its;
348       ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
349     } else if (NTL_KSP_STCG == tl->ksp_type) {
350       ierr = KSPSTCGSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr);
351       ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
352       ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
353       tao->ksp_its+=its;
354       tao->ksp_tot_its+=its;
355       ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
356     } else { /* NTL_KSP_GLTR */
357       ierr = KSPGLTRSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr);
358       ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
359       ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
360       tao->ksp_its+=its;
361       tao->ksp_tot_its+=its;
362       ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
363     }
364 
365     if (0.0 == tao->trust) {
366       /* Radius was uninitialized; use the norm of the direction */
367       if (norm_d > 0.0) {
368         tao->trust = norm_d;
369 
370         /* Modify the radius if it is too large or small */
371         tao->trust = PetscMax(tao->trust, tl->min_radius);
372         tao->trust = PetscMin(tao->trust, tl->max_radius);
373       } else {
374         /* The direction was bad; set radius to default value and re-solve
375            the trust-region subproblem to get a direction */
376         tao->trust = tao->trust0;
377 
378         /* Modify the radius if it is too large or small */
379         tao->trust = PetscMax(tao->trust, tl->min_radius);
380         tao->trust = PetscMin(tao->trust, tl->max_radius);
381 
382         if (NTL_KSP_NASH == tl->ksp_type) {
383           ierr = KSPNASHSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr);
384           ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
385           ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
386           tao->ksp_its+=its;
387           tao->ksp_tot_its+=its;
388           ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
389         } else if (NTL_KSP_STCG == tl->ksp_type) {
390           ierr = KSPSTCGSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr);
391           ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
392           ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
393           tao->ksp_its+=its;
394           tao->ksp_tot_its+=its;
395           ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
396         } else { /* NTL_KSP_GLTR */
397           ierr = KSPGLTRSetRadius(tao->ksp,tl->max_radius);CHKERRQ(ierr);
398           ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
399           ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
400           tao->ksp_its+=its;
401           tao->ksp_tot_its+=its;
402           ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
403         }
404 
405 
406         if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
407       }
408     }
409 
410     ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
411     ierr = KSPGetConvergedReason(tao->ksp, &ksp_reason);CHKERRQ(ierr);
412     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (NTL_PC_BFGS == tl->pc_type) && (bfgsUpdates > 1)) {
413       /* Preconditioner is numerically indefinite; reset the
414          approximate if using BFGS preconditioning. */
415 
416       if (f != 0.0) {
417         delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
418       } else {
419         delta = 2.0 / (gnorm*gnorm);
420       }
421       ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr);
422       ierr = MatLMVMReset(tl->M);CHKERRQ(ierr);
423       ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr);
424       bfgsUpdates = 1;
425     }
426 
427     /* Check trust-region reduction conditions */
428     tr_reject = 0;
429     if (NTL_UPDATE_REDUCTION == tl->update_type) {
430       /* Get predicted reduction */
431       if (NTL_KSP_NASH == tl->ksp_type) {
432         ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
433       } else if (NTL_KSP_STCG == tl->ksp_type) {
434         ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
435       } else { /* gltr */
436         ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
437       }
438 
439       if (prered >= 0.0) {
440         /* The predicted reduction has the wrong sign.  This cannot
441            happen in infinite precision arithmetic.  Step should
442            be rejected! */
443         tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
444         tr_reject = 1;
445       } else {
446         /* Compute trial step and function value */
447         ierr = VecCopy(tao->solution, tl->W);CHKERRQ(ierr);
448         ierr = VecAXPY(tl->W, 1.0, tao->stepdirection);CHKERRQ(ierr);
449         ierr = TaoComputeObjective(tao, tl->W, &ftrial);CHKERRQ(ierr);
450 
451         if (PetscIsInfOrNanReal(ftrial)) {
452           tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
453           tr_reject = 1;
454         } else {
455           /* Compute and actual reduction */
456           actred = f - ftrial;
457           prered = -prered;
458           if ((PetscAbsScalar(actred) <= tl->epsilon) &&
459               (PetscAbsScalar(prered) <= tl->epsilon)) {
460             kappa = 1.0;
461           } else {
462             kappa = actred / prered;
463           }
464 
465           /* Accept of reject the step and update radius */
466           if (kappa < tl->eta1) {
467             /* Reject the step */
468             tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
469             tr_reject = 1;
470           } else {
471             /* Accept the step */
472             if (kappa < tl->eta2) {
473               /* Marginal bad step */
474               tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
475             } else if (kappa < tl->eta3) {
476               /* Reasonable step */
477               tao->trust = tl->alpha3 * tao->trust;
478             } else if (kappa < tl->eta4) {
479               /* Good step */
480               tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
481             } else {
482               /* Very good step */
483               tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
484             }
485           }
486         }
487       }
488     } else {
489       /* Get predicted reduction */
490       if (NTL_KSP_NASH == tl->ksp_type) {
491         ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
492       } else if (NTL_KSP_STCG == tl->ksp_type) {
493         ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
494       } else { /* gltr */
495         ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
496       }
497 
498       if (prered >= 0.0) {
499         /* The predicted reduction has the wrong sign.  This cannot
500            happen in infinite precision arithmetic.  Step should
501            be rejected! */
502         tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
503         tr_reject = 1;
504       } else {
505         ierr = VecCopy(tao->solution, tl->W);CHKERRQ(ierr);
506         ierr = VecAXPY(tl->W, 1.0, tao->stepdirection);CHKERRQ(ierr);
507         ierr = TaoComputeObjective(tao, tl->W, &ftrial);CHKERRQ(ierr);
508         if (PetscIsInfOrNanReal(ftrial)) {
509           tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
510           tr_reject = 1;
511         } else {
512           ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr);
513 
514           actred = f - ftrial;
515           prered = -prered;
516           if ((PetscAbsScalar(actred) <= tl->epsilon) &&
517               (PetscAbsScalar(prered) <= tl->epsilon)) {
518             kappa = 1.0;
519           } else {
520             kappa = actred / prered;
521           }
522 
523           tau_1 = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
524           tau_2 = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
525           tau_min = PetscMin(tau_1, tau_2);
526           tau_max = PetscMax(tau_1, tau_2);
527 
528           if (kappa >= 1.0 - tl->mu1) {
529             /* Great agreement; accept step and update radius */
530             if (tau_max < 1.0) {
531               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
532             } else if (tau_max > tl->gamma4) {
533               tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
534             } else {
535               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
536             }
537           } else if (kappa >= 1.0 - tl->mu2) {
538             /* Good agreement */
539 
540             if (tau_max < tl->gamma2) {
541               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
542             } else if (tau_max > tl->gamma3) {
543               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
544             } else if (tau_max < 1.0) {
545               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
546             } else {
547               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
548             }
549           } else {
550             /* Not good agreement */
551             if (tau_min > 1.0) {
552               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
553             } else if (tau_max < tl->gamma1) {
554               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
555             } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
556               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
557             } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
558               tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
559             } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
560               tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
561             } else {
562               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
563             }
564             tr_reject = 1;
565           }
566         }
567       }
568     }
569 
570     if (tr_reject) {
571       /* The trust-region constraints rejected the step.  Apply a linesearch.
572          Check for descent direction. */
573       ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr);
574       if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
575         /* Newton step is not descent or direction produced Inf or NaN */
576 
577         if (NTL_PC_BFGS != tl->pc_type) {
578           /* We don't have the bfgs matrix around and updated
579              Must use gradient direction in this case */
580           ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
581           ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
582           ++tl->grad;
583           stepType = NTL_GRADIENT;
584         } else {
585           /* Attempt to use the BFGS direction */
586           ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
587           ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
588 
589           /* Check for success (descent direction) */
590           ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr);
591           if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
592             /* BFGS direction is not descent or direction produced not a number
593                We can assert bfgsUpdates > 1 in this case because
594                the first solve produces the scaled gradient direction,
595                which is guaranteed to be descent */
596 
597             /* Use steepest descent direction (scaled) */
598             if (f != 0.0) {
599               delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
600             } else {
601               delta = 2.0 / (gnorm*gnorm);
602             }
603             ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr);
604             ierr = MatLMVMReset(tl->M);CHKERRQ(ierr);
605             ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr);
606             ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
607             ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
608 
609             bfgsUpdates = 1;
610             ++tl->sgrad;
611             stepType = NTL_SCALED_GRADIENT;
612           } else {
613             if (1 == bfgsUpdates) {
614               /* The first BFGS direction is always the scaled gradient */
615               ++tl->sgrad;
616               stepType = NTL_SCALED_GRADIENT;
617             } else {
618               ++tl->bfgs;
619               stepType = NTL_BFGS;
620             }
621           }
622         }
623       } else {
624         /* Computed Newton step is descent */
625         ++tl->newt;
626         stepType = NTL_NEWTON;
627       }
628 
629       /* Perform the linesearch */
630       fold = f;
631       ierr = VecCopy(tao->solution, tl->Xold);CHKERRQ(ierr);
632       ierr = VecCopy(tao->gradient, tl->Gold);CHKERRQ(ierr);
633 
634       step = 1.0;
635       ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);CHKERRQ(ierr);
636       ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
637 
638       while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) {      /* Linesearch failed */
639         /* Linesearch failed */
640         f = fold;
641         ierr = VecCopy(tl->Xold, tao->solution);CHKERRQ(ierr);
642         ierr = VecCopy(tl->Gold, tao->gradient);CHKERRQ(ierr);
643 
644         switch(stepType) {
645         case NTL_NEWTON:
646           /* Failed to obtain acceptable iterate with Newton step */
647 
648           if (NTL_PC_BFGS != tl->pc_type) {
649             /* We don't have the bfgs matrix around and being updated
650                Must use gradient direction in this case */
651             ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
652             ++tl->grad;
653             stepType = NTL_GRADIENT;
654           } else {
655             /* Attempt to use the BFGS direction */
656             ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
657 
658 
659             /* Check for success (descent direction) */
660             ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr);
661             if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
662               /* BFGS direction is not descent or direction produced
663                  not a number.  We can assert bfgsUpdates > 1 in this case
664                  Use steepest descent direction (scaled) */
665 
666               if (f != 0.0) {
667                 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
668               } else {
669                 delta = 2.0 / (gnorm*gnorm);
670               }
671               ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr);
672               ierr = MatLMVMReset(tl->M);CHKERRQ(ierr);
673               ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr);
674               ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
675 
676               bfgsUpdates = 1;
677               ++tl->sgrad;
678               stepType = NTL_SCALED_GRADIENT;
679             } else {
680               if (1 == bfgsUpdates) {
681                 /* The first BFGS direction is always the scaled gradient */
682                 ++tl->sgrad;
683                 stepType = NTL_SCALED_GRADIENT;
684               } else {
685                 ++tl->bfgs;
686                 stepType = NTL_BFGS;
687               }
688             }
689           }
690           break;
691 
692         case NTL_BFGS:
693           /* Can only enter if pc_type == NTL_PC_BFGS
694              Failed to obtain acceptable iterate with BFGS step
695              Attempt to use the scaled gradient direction */
696 
697           if (f != 0.0) {
698             delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
699           } else {
700             delta = 2.0 / (gnorm*gnorm);
701           }
702           ierr = MatLMVMSetDelta(tl->M, delta);CHKERRQ(ierr);
703           ierr = MatLMVMReset(tl->M);CHKERRQ(ierr);
704           ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr);
705           ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
706 
707           bfgsUpdates = 1;
708           ++tl->sgrad;
709           stepType = NTL_SCALED_GRADIENT;
710           break;
711 
712         case NTL_SCALED_GRADIENT:
713           /* Can only enter if pc_type == NTL_PC_BFGS
714              The scaled gradient step did not produce a new iterate;
715              attemp to use the gradient direction.
716              Need to make sure we are not using a different diagonal scaling */
717           ierr = MatLMVMSetScale(tl->M, tl->Diag);CHKERRQ(ierr);
718           ierr = MatLMVMSetDelta(tl->M, 1.0);CHKERRQ(ierr);
719           ierr = MatLMVMReset(tl->M);CHKERRQ(ierr);
720           ierr = MatLMVMUpdate(tl->M, tao->solution, tao->gradient);CHKERRQ(ierr);
721           ierr = MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
722 
723           bfgsUpdates = 1;
724           ++tl->grad;
725           stepType = NTL_GRADIENT;
726           break;
727         }
728         ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
729 
730         /* This may be incorrect; linesearch has values for stepmax and stepmin
731            that should be reset. */
732         step = 1.0;
733         ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);CHKERRQ(ierr);
734         ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
735       }
736 
737       if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
738         /* Failed to find an improving point */
739         f = fold;
740         ierr = VecCopy(tl->Xold, tao->solution);CHKERRQ(ierr);
741         ierr = VecCopy(tl->Gold, tao->gradient);CHKERRQ(ierr);
742         tao->trust = 0.0;
743         step = 0.0;
744         reason = TAO_DIVERGED_LS_FAILURE;
745         tao->reason = TAO_DIVERGED_LS_FAILURE;
746         break;
747       } else if (stepType == NTL_NEWTON) {
748         if (step < tl->nu1) {
749           /* Very bad step taken; reduce radius */
750           tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
751         } else if (step < tl->nu2) {
752           /* Reasonably bad step taken; reduce radius */
753           tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
754         } else if (step < tl->nu3) {
755           /* Reasonable step was taken; leave radius alone */
756           if (tl->omega3 < 1.0) {
757             tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
758           } else if (tl->omega3 > 1.0) {
759             tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
760           }
761         } else if (step < tl->nu4) {
762           /* Full step taken; increase the radius */
763           tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
764         } else {
765           /* More than full step taken; increase the radius */
766           tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
767         }
768       } else {
769         /* Newton step was not good; reduce the radius */
770         tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
771       }
772     } else {
773       /* Trust-region step is accepted */
774       ierr = VecCopy(tl->W, tao->solution);CHKERRQ(ierr);
775       f = ftrial;
776       ierr = TaoComputeGradient(tao, tao->solution, tao->gradient);CHKERRQ(ierr);
777       ++tl->ntrust;
778     }
779 
780     /* The radius may have been increased; modify if it is too large */
781     tao->trust = PetscMin(tao->trust, tl->max_radius);
782 
783     /* Check for converged */
784     ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr);
785     if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number");
786     needH = 1;
787 
788     ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust, &reason);CHKERRQ(ierr);
789   }
790   PetscFunctionReturn(0);
791 }
792 
793 /* ---------------------------------------------------------- */
794 #undef __FUNCT__
795 #define __FUNCT__ "TaoSetUp_NTL"
796 static PetscErrorCode TaoSetUp_NTL(Tao tao)
797 {
798   TAO_NTL        *tl = (TAO_NTL *)tao->data;
799   PetscErrorCode ierr;
800 
801   PetscFunctionBegin;
802   if (!tao->gradient) {ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr); }
803   if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);}
804   if (!tl->W) { ierr = VecDuplicate(tao->solution, &tl->W);CHKERRQ(ierr);}
805   if (!tl->Xold) { ierr = VecDuplicate(tao->solution, &tl->Xold);CHKERRQ(ierr);}
806   if (!tl->Gold) { ierr = VecDuplicate(tao->solution, &tl->Gold);CHKERRQ(ierr);}
807   tl->Diag = 0;
808   tl->M = 0;
809   PetscFunctionReturn(0);
810 }
811 
812 /*------------------------------------------------------------*/
813 #undef __FUNCT__
814 #define __FUNCT__ "TaoDestroy_NTL"
815 static PetscErrorCode TaoDestroy_NTL(Tao tao)
816 {
817   TAO_NTL        *tl = (TAO_NTL *)tao->data;
818   PetscErrorCode ierr;
819 
820   PetscFunctionBegin;
821   if (tao->setupcalled) {
822     ierr = VecDestroy(&tl->W);CHKERRQ(ierr);
823     ierr = VecDestroy(&tl->Xold);CHKERRQ(ierr);
824     ierr = VecDestroy(&tl->Gold);CHKERRQ(ierr);
825   }
826   ierr = VecDestroy(&tl->Diag);CHKERRQ(ierr);
827   ierr = MatDestroy(&tl->M);CHKERRQ(ierr);
828   ierr = PetscFree(tao->data);CHKERRQ(ierr);
829   PetscFunctionReturn(0);
830 }
831 
832 /*------------------------------------------------------------*/
833 #undef __FUNCT__
834 #define __FUNCT__ "TaoSetFromOptions_NTL"
835 static PetscErrorCode TaoSetFromOptions_NTL(PetscOptionItems *PetscOptionsObject,Tao tao)
836 {
837   TAO_NTL        *tl = (TAO_NTL *)tao->data;
838   PetscErrorCode ierr;
839 
840   PetscFunctionBegin;
841   ierr = PetscOptionsHead(PetscOptionsObject,"Newton trust region with line search method for unconstrained optimization");CHKERRQ(ierr);
842   ierr = PetscOptionsEList("-tao_ntl_ksp_type", "ksp type", "", NTL_KSP, NTL_KSP_TYPES, NTL_KSP[tl->ksp_type], &tl->ksp_type,NULL);CHKERRQ(ierr);
843   ierr = PetscOptionsEList("-tao_ntl_pc_type", "pc type", "", NTL_PC, NTL_PC_TYPES, NTL_PC[tl->pc_type], &tl->pc_type,NULL);CHKERRQ(ierr);
844   ierr = PetscOptionsEList("-tao_ntl_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tl->bfgs_scale_type], &tl->bfgs_scale_type,NULL);CHKERRQ(ierr);
845   ierr = PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type,NULL);CHKERRQ(ierr);
846   ierr = PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type,NULL);CHKERRQ(ierr);
847   ierr = PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1,NULL);CHKERRQ(ierr);
848   ierr = PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2,NULL);CHKERRQ(ierr);
849   ierr = PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3,NULL);CHKERRQ(ierr);
850   ierr = PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4,NULL);CHKERRQ(ierr);
851   ierr = PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1,NULL);CHKERRQ(ierr);
852   ierr = PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2,NULL);CHKERRQ(ierr);
853   ierr = PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3,NULL);CHKERRQ(ierr);
854   ierr = PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4,NULL);CHKERRQ(ierr);
855   ierr = PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5,NULL);CHKERRQ(ierr);
856   ierr = PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1,NULL);CHKERRQ(ierr);
857   ierr = PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2,NULL);CHKERRQ(ierr);
858   ierr = PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3,NULL);CHKERRQ(ierr);
859   ierr = PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4,NULL);CHKERRQ(ierr);
860   ierr = PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1,NULL);CHKERRQ(ierr);
861   ierr = PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2,NULL);CHKERRQ(ierr);
862   ierr = PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3,NULL);CHKERRQ(ierr);
863   ierr = PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4,NULL);CHKERRQ(ierr);
864   ierr = PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5,NULL);CHKERRQ(ierr);
865   ierr = PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i,NULL);CHKERRQ(ierr);
866   ierr = PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i,NULL);CHKERRQ(ierr);
867   ierr = PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i,NULL);CHKERRQ(ierr);
868   ierr = PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i,NULL);CHKERRQ(ierr);
869   ierr = PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i,NULL);CHKERRQ(ierr);
870   ierr = PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i,NULL);CHKERRQ(ierr);
871   ierr = PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i,NULL);CHKERRQ(ierr);
872   ierr = PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1,NULL);CHKERRQ(ierr);
873   ierr = PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2,NULL);CHKERRQ(ierr);
874   ierr = PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1,NULL);CHKERRQ(ierr);
875   ierr = PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2,NULL);CHKERRQ(ierr);
876   ierr = PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3,NULL);CHKERRQ(ierr);
877   ierr = PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4,NULL);CHKERRQ(ierr);
878   ierr = PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta,NULL);CHKERRQ(ierr);
879   ierr = PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius,NULL);CHKERRQ(ierr);
880   ierr = PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius,NULL);CHKERRQ(ierr);
881   ierr = PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon,NULL);CHKERRQ(ierr);
882   ierr = PetscOptionsTail();CHKERRQ(ierr);
883   ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
884   ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
885   PetscFunctionReturn(0);
886 }
887 
888 /*------------------------------------------------------------*/
889 #undef __FUNCT__
890 #define __FUNCT__ "TaoView_NTL"
891 static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer)
892 {
893   TAO_NTL        *tl = (TAO_NTL *)tao->data;
894   PetscInt       nrejects;
895   PetscBool      isascii;
896   PetscErrorCode ierr;
897 
898   PetscFunctionBegin;
899   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
900   if (isascii) {
901     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
902     if (NTL_PC_BFGS == tl->pc_type && tl->M) {
903       ierr = MatLMVMGetRejects(tl->M, &nrejects);CHKERRQ(ierr);
904       ierr = PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n", nrejects);CHKERRQ(ierr);
905     }
906     ierr = PetscViewerASCIIPrintf(viewer, "Trust-region steps: %D\n", tl->ntrust);CHKERRQ(ierr);
907     ierr = PetscViewerASCIIPrintf(viewer, "Newton search steps: %D\n", tl->newt);CHKERRQ(ierr);
908     ierr = PetscViewerASCIIPrintf(viewer, "BFGS search steps: %D\n", tl->bfgs);CHKERRQ(ierr);
909     ierr = PetscViewerASCIIPrintf(viewer, "Scaled gradient search steps: %D\n", tl->sgrad);CHKERRQ(ierr);
910     ierr = PetscViewerASCIIPrintf(viewer, "Gradient search steps: %D\n", tl->grad);CHKERRQ(ierr);
911     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
912   }
913   PetscFunctionReturn(0);
914 }
915 
916 /* ---------------------------------------------------------- */
917 /*MC
918   TAONTR - Newton's method with trust region and linesearch
919   for unconstrained minimization.
920   At each iteration, the Newton trust region method solves the system for d
921   and performs a line search in the d direction:
922 
923             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k
924 
925   Options Database Keys:
926 + -tao_ntl_ksp_type - "nash","stcg","gltr"
927 . -tao_ntl_pc_type - "none","ahess","bfgs","petsc"
928 . -tao_ntl_bfgs_scale_type - type of scaling with bfgs pc, "ahess" or "bfgs"
929 . -tao_ntl_init_type - "constant","direction","interpolation"
930 . -tao_ntl_update_type - "reduction","interpolation"
931 . -tao_ntl_min_radius - lower bound on trust region radius
932 . -tao_ntl_max_radius - upper bound on trust region radius
933 . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction
934 . -tao_ntl_mu1_i - mu1 interpolation init factor
935 . -tao_ntl_mu2_i - mu2 interpolation init factor
936 . -tao_ntl_gamma1_i - gamma1 interpolation init factor
937 . -tao_ntl_gamma2_i - gamma2 interpolation init factor
938 . -tao_ntl_gamma3_i - gamma3 interpolation init factor
939 . -tao_ntl_gamma4_i - gamma4 interpolation init factor
940 . -tao_ntl_theta_i - thetha1 interpolation init factor
941 . -tao_ntl_eta1 - eta1 reduction update factor
942 . -tao_ntl_eta2 - eta2 reduction update factor
943 . -tao_ntl_eta3 - eta3 reduction update factor
944 . -tao_ntl_eta4 - eta4 reduction update factor
945 . -tao_ntl_alpha1 - alpha1 reduction update factor
946 . -tao_ntl_alpha2 - alpha2 reduction update factor
947 . -tao_ntl_alpha3 - alpha3 reduction update factor
948 . -tao_ntl_alpha4 - alpha4 reduction update factor
949 . -tao_ntl_alpha4 - alpha4 reduction update factor
950 . -tao_ntl_mu1 - mu1 interpolation update
951 . -tao_ntl_mu2 - mu2 interpolation update
952 . -tao_ntl_gamma1 - gamma1 interpolcation update
953 . -tao_ntl_gamma2 - gamma2 interpolcation update
954 . -tao_ntl_gamma3 - gamma3 interpolcation update
955 . -tao_ntl_gamma4 - gamma4 interpolation update
956 - -tao_ntl_theta - theta1 interpolation update
957 
958   Level: beginner
959 M*/
960 
961 #undef __FUNCT__
962 #define __FUNCT__ "TaoCreate_NTL"
963 PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao)
964 {
965   TAO_NTL        *tl;
966   PetscErrorCode ierr;
967   const char     *morethuente_type = TAOLINESEARCHMT;
968 
969   PetscFunctionBegin;
970   ierr = PetscNewLog(tao,&tl);CHKERRQ(ierr);
971   tao->ops->setup = TaoSetUp_NTL;
972   tao->ops->solve = TaoSolve_NTL;
973   tao->ops->view = TaoView_NTL;
974   tao->ops->setfromoptions = TaoSetFromOptions_NTL;
975   tao->ops->destroy = TaoDestroy_NTL;
976 
977   /* Override default settings (unless already changed) */
978   if (!tao->max_it_changed) tao->max_it = 50;
979   if (!tao->trust0_changed) tao->trust0 = 100.0;
980 
981   tao->data = (void*)tl;
982 
983   /* Default values for trust-region radius update based on steplength */
984   tl->nu1 = 0.25;
985   tl->nu2 = 0.50;
986   tl->nu3 = 1.00;
987   tl->nu4 = 1.25;
988 
989   tl->omega1 = 0.25;
990   tl->omega2 = 0.50;
991   tl->omega3 = 1.00;
992   tl->omega4 = 2.00;
993   tl->omega5 = 4.00;
994 
995   /* Default values for trust-region radius update based on reduction */
996   tl->eta1 = 1.0e-4;
997   tl->eta2 = 0.25;
998   tl->eta3 = 0.50;
999   tl->eta4 = 0.90;
1000 
1001   tl->alpha1 = 0.25;
1002   tl->alpha2 = 0.50;
1003   tl->alpha3 = 1.00;
1004   tl->alpha4 = 2.00;
1005   tl->alpha5 = 4.00;
1006 
1007   /* Default values for trust-region radius update based on interpolation */
1008   tl->mu1 = 0.10;
1009   tl->mu2 = 0.50;
1010 
1011   tl->gamma1 = 0.25;
1012   tl->gamma2 = 0.50;
1013   tl->gamma3 = 2.00;
1014   tl->gamma4 = 4.00;
1015 
1016   tl->theta = 0.05;
1017 
1018   /* Default values for trust region initialization based on interpolation */
1019   tl->mu1_i = 0.35;
1020   tl->mu2_i = 0.50;
1021 
1022   tl->gamma1_i = 0.0625;
1023   tl->gamma2_i = 0.5;
1024   tl->gamma3_i = 2.0;
1025   tl->gamma4_i = 5.0;
1026 
1027   tl->theta_i = 0.25;
1028 
1029   /* Remaining parameters */
1030   tl->min_radius = 1.0e-10;
1031   tl->max_radius = 1.0e10;
1032   tl->epsilon = 1.0e-6;
1033 
1034   tl->ksp_type        = NTL_KSP_STCG;
1035   tl->pc_type         = NTL_PC_BFGS;
1036   tl->bfgs_scale_type = BFGS_SCALE_AHESS;
1037   tl->init_type       = NTL_INIT_INTERPOLATION;
1038   tl->update_type     = NTL_UPDATE_REDUCTION;
1039 
1040   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
1041   ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type);CHKERRQ(ierr);
1042   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch, tao);CHKERRQ(ierr);
1043   ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
1044   ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr);
1045   ierr = KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);CHKERRQ(ierr);
1046   PetscFunctionReturn(0);
1047 }
1048 
1049 
1050 
1051 
1052