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