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