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