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