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