xref: /petsc/src/tao/unconstrained/impls/ntr/ntr.c (revision 65f8aed5f7eaa1e2ef2ddeffe666264e0669c876) 
1 #include <../src/tao/unconstrained/impls/ntr/ntrimpl.h>
2 
3 #include <petscksp.h>
4 
5 #define NTR_INIT_CONSTANT         0
6 #define NTR_INIT_DIRECTION        1
7 #define NTR_INIT_INTERPOLATION    2
8 #define NTR_INIT_TYPES            3
9 
10 #define NTR_UPDATE_REDUCTION      0
11 #define NTR_UPDATE_INTERPOLATION  1
12 #define NTR_UPDATE_TYPES          2
13 
14 static const char *NTR_INIT[64] = {"constant","direction","interpolation"};
15 
16 static const char *NTR_UPDATE[64] = {"reduction","interpolation"};
17 
18 /*
19    TaoSolve_NTR - Implements Newton's Method with a trust region approach
20    for solving unconstrained minimization problems.
21 
22    The basic algorithm is taken from MINPACK-2 (dstrn).
23 
24    TaoSolve_NTR computes a local minimizer of a twice differentiable function
25    f by applying a trust region variant of Newton's method.  At each stage
26    of the algorithm, we use the prconditioned conjugate gradient method to
27    determine an approximate minimizer of the quadratic equation
28 
29         q(s) = <s, Hs + g>
30 
31    subject to the trust region constraint
32 
33         || s ||_M <= radius,
34 
35    where radius is the trust region radius and M is a symmetric positive
36    definite matrix (the preconditioner).  Here g is the gradient and H
37    is the Hessian matrix.
38 
39    Note:  TaoSolve_NTR MUST use the iterative solver KSPCGNASH, KSPCGSTCG,
40           or KSPCGGLTR.  Thus, we set KSPCGNASH, KSPCGSTCG, or KSPCGGLTR in this
41           routine regardless of what the user may have previously specified.
42 */
43 static PetscErrorCode TaoSolve_NTR(Tao tao)
44 {
45   TAO_NTR            *tr = (TAO_NTR *)tao->data;
46   KSPType            ksp_type;
47   PetscBool          is_nash,is_stcg,is_gltr,is_bfgs,is_jacobi,is_symmetric,sym_set;
48   KSPConvergedReason ksp_reason;
49   PC                 pc;
50   PetscReal          fmin, ftrial, prered, actred, kappa, sigma, beta;
51   PetscReal          tau, tau_1, tau_2, tau_max, tau_min, max_radius;
52   PetscReal          f, gnorm;
53 
54   PetscReal          norm_d;
55   PetscErrorCode     ierr;
56   PetscInt           bfgsUpdates = 0;
57   PetscInt           needH;
58 
59   PetscInt           i_max = 5;
60   PetscInt           j_max = 1;
61   PetscInt           i, j, N, n, its;
62 
63   PetscFunctionBegin;
64   if (tao->XL || tao->XU || tao->ops->computebounds) {
65     ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");CHKERRQ(ierr);
66   }
67 
68   ierr = KSPGetType(tao->ksp,&ksp_type);CHKERRQ(ierr);
69   ierr = PetscStrcmp(ksp_type,KSPCGNASH,&is_nash);CHKERRQ(ierr);
70   ierr = PetscStrcmp(ksp_type,KSPCGSTCG,&is_stcg);CHKERRQ(ierr);
71   ierr = PetscStrcmp(ksp_type,KSPCGGLTR,&is_gltr);CHKERRQ(ierr);
72   if (!is_nash && !is_stcg && !is_gltr) {
73     SETERRQ(PETSC_COMM_SELF,1,"TAO_NTR requires nash, stcg, or gltr for the KSP");
74   }
75 
76   /* Initialize the radius and modify if it is too large or small */
77   tao->trust = tao->trust0;
78   tao->trust = PetscMax(tao->trust, tr->min_radius);
79   tao->trust = PetscMin(tao->trust, tr->max_radius);
80 
81 /* Allocate the vectors needed for the BFGS approximation */
82 ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr);
83 ierr = PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);CHKERRQ(ierr);
84 ierr = PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);CHKERRQ(ierr);
85 if (is_bfgs) {
86   tr->bfgs_pre = pc;
87   ierr = PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M);CHKERRQ(ierr);
88   ierr = VecGetLocalSize(tao->solution, &n);CHKERRQ(ierr);
89   ierr = VecGetSize(tao->solution, &N);CHKERRQ(ierr);
90   ierr = MatSetSizes(tr->M, n, n, N, N);CHKERRQ(ierr);
91   ierr = MatLMVMAllocate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr);
92   ierr = MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric);CHKERRQ(ierr);
93   if (!sym_set || !is_symmetric) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
94 } else if (is_jacobi) {
95   ierr = PCJacobiSetUseAbs(pc,PETSC_TRUE);CHKERRQ(ierr);
96 }
97 
98   /* Check convergence criteria */
99   ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
100   ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
101   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Inf or NaN");
102   needH = 1;
103 
104   tao->reason = TAO_CONTINUE_ITERATING;
105   ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
106   ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);CHKERRQ(ierr);
107   ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
108   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
109 
110   /*  Initialize trust-region radius */
111   switch(tr->init_type) {
112   case NTR_INIT_CONSTANT:
113     /*  Use the initial radius specified */
114     break;
115 
116   case NTR_INIT_INTERPOLATION:
117     /*  Use the initial radius specified */
118     max_radius = 0.0;
119 
120     for (j = 0; j < j_max; ++j) {
121       fmin = f;
122       sigma = 0.0;
123 
124       if (needH) {
125         ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
126         needH = 0;
127       }
128 
129       for (i = 0; i < i_max; ++i) {
130 
131         ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr);
132         ierr = VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);CHKERRQ(ierr);
133         ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);
134 
135         if (PetscIsInfOrNanReal(ftrial)) {
136           tau = tr->gamma1_i;
137         }
138         else {
139           if (ftrial < fmin) {
140             fmin = ftrial;
141             sigma = -tao->trust / gnorm;
142           }
143 
144           ierr = MatMult(tao->hessian, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
145           ierr = VecDot(tao->gradient, tao->stepdirection, &prered);CHKERRQ(ierr);
146 
147           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
148           actred = f - ftrial;
149           if ((PetscAbsScalar(actred) <= tr->epsilon) &&
150               (PetscAbsScalar(prered) <= tr->epsilon)) {
151             kappa = 1.0;
152           }
153           else {
154             kappa = actred / prered;
155           }
156 
157           tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
158           tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
159           tau_min = PetscMin(tau_1, tau_2);
160           tau_max = PetscMax(tau_1, tau_2);
161 
162           if (PetscAbsScalar(kappa - 1.0) <= tr->mu1_i) {
163             /*  Great agreement */
164             max_radius = PetscMax(max_radius, tao->trust);
165 
166             if (tau_max < 1.0) {
167               tau = tr->gamma3_i;
168             }
169             else if (tau_max > tr->gamma4_i) {
170               tau = tr->gamma4_i;
171             }
172             else {
173               tau = tau_max;
174             }
175           }
176           else if (PetscAbsScalar(kappa - 1.0) <= tr->mu2_i) {
177             /*  Good agreement */
178             max_radius = PetscMax(max_radius, tao->trust);
179 
180             if (tau_max < tr->gamma2_i) {
181               tau = tr->gamma2_i;
182             }
183             else if (tau_max > tr->gamma3_i) {
184               tau = tr->gamma3_i;
185             }
186             else {
187               tau = tau_max;
188             }
189           }
190           else {
191             /*  Not good agreement */
192             if (tau_min > 1.0) {
193               tau = tr->gamma2_i;
194             }
195             else if (tau_max < tr->gamma1_i) {
196               tau = tr->gamma1_i;
197             }
198             else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
199               tau = tr->gamma1_i;
200             }
201             else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) &&
202                      ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
203               tau = tau_1;
204             }
205             else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) &&
206                      ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
207               tau = tau_2;
208             }
209             else {
210               tau = tau_max;
211             }
212           }
213         }
214         tao->trust = tau * tao->trust;
215       }
216 
217       if (fmin < f) {
218         f = fmin;
219         ierr = VecAXPY(tao->solution, sigma, tao->gradient);CHKERRQ(ierr);
220         ierr = TaoComputeGradient(tao,tao->solution, tao->gradient);CHKERRQ(ierr);
221 
222         ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
223 
224         if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
225         needH = 1;
226 
227         ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
228         ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);CHKERRQ(ierr);
229         ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
230         if (tao->reason != TAO_CONTINUE_ITERATING) {
231           PetscFunctionReturn(0);
232         }
233       }
234     }
235     tao->trust = PetscMax(tao->trust, max_radius);
236 
237     /*  Modify the radius if it is too large or small */
238     tao->trust = PetscMax(tao->trust, tr->min_radius);
239     tao->trust = PetscMin(tao->trust, tr->max_radius);
240     break;
241 
242   default:
243     /*  Norm of the first direction will initialize radius */
244     tao->trust = 0.0;
245     break;
246   }
247 
248   /* Have not converged; continue with Newton method */
249   while (tao->reason == TAO_CONTINUE_ITERATING) {
250     ++tao->niter;
251     tao->ksp_its=0;
252     /* Compute the Hessian */
253     if (needH) {
254       ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
255       needH = 0;
256     }
257 
258     if (tr->bfgs_pre) {
259       /* Update the limited memory preconditioner */
260       ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr);
261       ++bfgsUpdates;
262     }
263 
264     while (tao->reason == TAO_CONTINUE_ITERATING) {
265       ierr = KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);CHKERRQ(ierr);
266 
267       /* Solve the trust region subproblem */
268       ierr = KSPCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
269       ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
270       ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
271       tao->ksp_its+=its;
272       tao->ksp_tot_its+=its;
273       ierr = KSPCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
274 
275       if (0.0 == tao->trust) {
276         /* Radius was uninitialized; use the norm of the direction */
277         if (norm_d > 0.0) {
278           tao->trust = norm_d;
279 
280           /* Modify the radius if it is too large or small */
281           tao->trust = PetscMax(tao->trust, tr->min_radius);
282           tao->trust = PetscMin(tao->trust, tr->max_radius);
283         }
284         else {
285           /* The direction was bad; set radius to default value and re-solve
286              the trust-region subproblem to get a direction */
287           tao->trust = tao->trust0;
288 
289           /* Modify the radius if it is too large or small */
290           tao->trust = PetscMax(tao->trust, tr->min_radius);
291           tao->trust = PetscMin(tao->trust, tr->max_radius);
292 
293           ierr = KSPCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
294           ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
295           ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
296           tao->ksp_its+=its;
297           tao->ksp_tot_its+=its;
298           ierr = KSPCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
299 
300           if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
301         }
302       }
303       ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
304       ierr = KSPGetConvergedReason(tao->ksp, &ksp_reason);CHKERRQ(ierr);
305       if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) {
306         /* Preconditioner is numerically indefinite; reset the
307            approximate if using BFGS preconditioning. */
308         ierr = MatLMVMReset(tr->M, PETSC_FALSE);CHKERRQ(ierr);
309         ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr);
310         bfgsUpdates = 1;
311       }
312 
313       if (NTR_UPDATE_REDUCTION == tr->update_type) {
314         /* Get predicted reduction */
315         ierr = KSPCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
316         if (prered >= 0.0) {
317           /* The predicted reduction has the wrong sign.  This cannot
318              happen in infinite precision arithmetic.  Step should
319              be rejected! */
320           tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
321         }
322         else {
323           /* Compute trial step and function value */
324           ierr = VecCopy(tao->solution,tr->W);CHKERRQ(ierr);
325           ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr);
326           ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);
327 
328           if (PetscIsInfOrNanReal(ftrial)) {
329             tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
330           } else {
331             /* Compute and actual reduction */
332             actred = f - ftrial;
333             prered = -prered;
334             if ((PetscAbsScalar(actred) <= tr->epsilon) &&
335                 (PetscAbsScalar(prered) <= tr->epsilon)) {
336               kappa = 1.0;
337             }
338             else {
339               kappa = actred / prered;
340             }
341 
342             /* Accept or reject the step and update radius */
343             if (kappa < tr->eta1) {
344               /* Reject the step */
345               tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
346             }
347             else {
348               /* Accept the step */
349               if (kappa < tr->eta2) {
350                 /* Marginal bad step */
351                 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
352               }
353               else if (kappa < tr->eta3) {
354                 /* Reasonable step */
355                 tao->trust = tr->alpha3 * tao->trust;
356               }
357               else if (kappa < tr->eta4) {
358                 /* Good step */
359                 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
360               }
361               else {
362                 /* Very good step */
363                 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
364               }
365               break;
366             }
367           }
368         }
369       }
370       else {
371         /* Get predicted reduction */
372         ierr = KSPCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
373         if (prered >= 0.0) {
374           /* The predicted reduction has the wrong sign.  This cannot
375              happen in infinite precision arithmetic.  Step should
376              be rejected! */
377           tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
378         }
379         else {
380           ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr);
381           ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr);
382           ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);
383           if (PetscIsInfOrNanReal(ftrial)) {
384             tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
385           }
386           else {
387             ierr = VecDot(tao->gradient, tao->stepdirection, &beta);CHKERRQ(ierr);
388             actred = f - ftrial;
389             prered = -prered;
390             if ((PetscAbsScalar(actred) <= tr->epsilon) &&
391                 (PetscAbsScalar(prered) <= tr->epsilon)) {
392               kappa = 1.0;
393             }
394             else {
395               kappa = actred / prered;
396             }
397 
398             tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
399             tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
400             tau_min = PetscMin(tau_1, tau_2);
401             tau_max = PetscMax(tau_1, tau_2);
402 
403             if (kappa >= 1.0 - tr->mu1) {
404               /* Great agreement; accept step and update radius */
405               if (tau_max < 1.0) {
406                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
407               }
408               else if (tau_max > tr->gamma4) {
409                 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
410               }
411               else {
412                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
413               }
414               break;
415             }
416             else if (kappa >= 1.0 - tr->mu2) {
417               /* Good agreement */
418 
419               if (tau_max < tr->gamma2) {
420                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
421               }
422               else if (tau_max > tr->gamma3) {
423                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
424               }
425               else if (tau_max < 1.0) {
426                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
427               }
428               else {
429                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
430               }
431               break;
432             }
433             else {
434               /* Not good agreement */
435               if (tau_min > 1.0) {
436                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
437               }
438               else if (tau_max < tr->gamma1) {
439                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
440               }
441               else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
442                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
443               }
444               else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) &&
445                        ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
446                 tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
447               }
448               else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) &&
449                        ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
450                 tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
451               }
452               else {
453                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
454               }
455             }
456           }
457         }
458       }
459 
460       /* The step computed was not good and the radius was decreased.
461          Monitor the radius to terminate. */
462       ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
463       ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);CHKERRQ(ierr);
464       ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
465     }
466 
467     /* The radius may have been increased; modify if it is too large */
468     tao->trust = PetscMin(tao->trust, tr->max_radius);
469 
470     if (tao->reason == TAO_CONTINUE_ITERATING) {
471       ierr = VecCopy(tr->W, tao->solution);CHKERRQ(ierr);
472       f = ftrial;
473       ierr = TaoComputeGradient(tao, tao->solution, tao->gradient);CHKERRQ(ierr);
474       ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
475       if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
476       needH = 1;
477       ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
478       ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);CHKERRQ(ierr);
479       ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
480     }
481   }
482   PetscFunctionReturn(0);
483 }
484 
485 /*------------------------------------------------------------*/
486 static PetscErrorCode TaoSetUp_NTR(Tao tao)
487 {
488   TAO_NTR *tr = (TAO_NTR *)tao->data;
489   PetscErrorCode ierr;
490 
491   PetscFunctionBegin;
492   if (!tao->gradient) {ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr);}
493   if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);}
494   if (!tr->W) {ierr = VecDuplicate(tao->solution, &tr->W);CHKERRQ(ierr);}
495 
496   tr->bfgs_pre = 0;
497   tr->M = 0;
498   PetscFunctionReturn(0);
499 }
500 
501 /*------------------------------------------------------------*/
502 static PetscErrorCode TaoDestroy_NTR(Tao tao)
503 {
504   TAO_NTR        *tr = (TAO_NTR *)tao->data;
505   PetscErrorCode ierr;
506 
507   PetscFunctionBegin;
508   if (tao->setupcalled) {
509     ierr = VecDestroy(&tr->W);CHKERRQ(ierr);
510   }
511   ierr = PetscFree(tao->data);CHKERRQ(ierr);
512   PetscFunctionReturn(0);
513 }
514 
515 /*------------------------------------------------------------*/
516 static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao)
517 {
518   TAO_NTR        *tr = (TAO_NTR *)tao->data;
519   PetscErrorCode ierr;
520 
521   PetscFunctionBegin;
522   ierr = PetscOptionsHead(PetscOptionsObject,"Newton trust region method for unconstrained optimization");CHKERRQ(ierr);
523   ierr = PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type,NULL);CHKERRQ(ierr);
524   ierr = PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type,NULL);CHKERRQ(ierr);
525   ierr = PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL);CHKERRQ(ierr);
526   ierr = PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL);CHKERRQ(ierr);
527   ierr = PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL);CHKERRQ(ierr);
528   ierr = PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL);CHKERRQ(ierr);
529   ierr = PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL);CHKERRQ(ierr);
530   ierr = PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL);CHKERRQ(ierr);
531   ierr = PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL);CHKERRQ(ierr);
532   ierr = PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL);CHKERRQ(ierr);
533   ierr = PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL);CHKERRQ(ierr);
534   ierr = PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL);CHKERRQ(ierr);
535   ierr = PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL);CHKERRQ(ierr);
536   ierr = PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL);CHKERRQ(ierr);
537   ierr = PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL);CHKERRQ(ierr);
538   ierr = PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL);CHKERRQ(ierr);
539   ierr = PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL);CHKERRQ(ierr);
540   ierr = PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL);CHKERRQ(ierr);
541   ierr = PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL);CHKERRQ(ierr);
542   ierr = PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL);CHKERRQ(ierr);
543   ierr = PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL);CHKERRQ(ierr);
544   ierr = PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL);CHKERRQ(ierr);
545   ierr = PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL);CHKERRQ(ierr);
546   ierr = PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL);CHKERRQ(ierr);
547   ierr = PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL);CHKERRQ(ierr);
548   ierr = PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL);CHKERRQ(ierr);
549   ierr = PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL);CHKERRQ(ierr);
550   ierr = PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL);CHKERRQ(ierr);
551   ierr = PetscOptionsTail();CHKERRQ(ierr);
552   ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
553   PetscFunctionReturn(0);
554 }
555 
556 /*------------------------------------------------------------*/
557 /*MC
558   TAONTR - Newton's method with trust region for unconstrained minimization.
559   At each iteration, the Newton trust region method solves the system.
560   NTR expects a KSP solver with a trust region radius.
561             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k
562 
563   Options Database Keys:
564 + -tao_ntr_pc_type - "none","ahess","bfgs","petsc"
565 . -tao_ntr_bfgs_scale_type - type of scaling with bfgs pc, "ahess" or "bfgs"
566 . -tao_ntr_init_type - "constant","direction","interpolation"
567 . -tao_ntr_update_type - "reduction","interpolation"
568 . -tao_ntr_min_radius - lower bound on trust region radius
569 . -tao_ntr_max_radius - upper bound on trust region radius
570 . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
571 . -tao_ntr_mu1_i - mu1 interpolation init factor
572 . -tao_ntr_mu2_i - mu2 interpolation init factor
573 . -tao_ntr_gamma1_i - gamma1 interpolation init factor
574 . -tao_ntr_gamma2_i - gamma2 interpolation init factor
575 . -tao_ntr_gamma3_i - gamma3 interpolation init factor
576 . -tao_ntr_gamma4_i - gamma4 interpolation init factor
577 . -tao_ntr_theta_i - thetha1 interpolation init factor
578 . -tao_ntr_eta1 - eta1 reduction update factor
579 . -tao_ntr_eta2 - eta2 reduction update factor
580 . -tao_ntr_eta3 - eta3 reduction update factor
581 . -tao_ntr_eta4 - eta4 reduction update factor
582 . -tao_ntr_alpha1 - alpha1 reduction update factor
583 . -tao_ntr_alpha2 - alpha2 reduction update factor
584 . -tao_ntr_alpha3 - alpha3 reduction update factor
585 . -tao_ntr_alpha4 - alpha4 reduction update factor
586 . -tao_ntr_alpha4 - alpha4 reduction update factor
587 . -tao_ntr_mu1 - mu1 interpolation update
588 . -tao_ntr_mu2 - mu2 interpolation update
589 . -tao_ntr_gamma1 - gamma1 interpolcation update
590 . -tao_ntr_gamma2 - gamma2 interpolcation update
591 . -tao_ntr_gamma3 - gamma3 interpolcation update
592 . -tao_ntr_gamma4 - gamma4 interpolation update
593 - -tao_ntr_theta - theta interpolation update
594 
595   Level: beginner
596 M*/
597 
598 PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
599 {
600   TAO_NTR *tr;
601   PetscErrorCode ierr;
602 
603   PetscFunctionBegin;
604 
605   ierr = PetscNewLog(tao,&tr);CHKERRQ(ierr);
606 
607   tao->ops->setup = TaoSetUp_NTR;
608   tao->ops->solve = TaoSolve_NTR;
609   tao->ops->setfromoptions = TaoSetFromOptions_NTR;
610   tao->ops->destroy = TaoDestroy_NTR;
611 
612   /* Override default settings (unless already changed) */
613   if (!tao->max_it_changed) tao->max_it = 50;
614   if (!tao->trust0_changed) tao->trust0 = 100.0;
615   tao->data = (void*)tr;
616 
617   /*  Standard trust region update parameters */
618   tr->eta1 = 1.0e-4;
619   tr->eta2 = 0.25;
620   tr->eta3 = 0.50;
621   tr->eta4 = 0.90;
622 
623   tr->alpha1 = 0.25;
624   tr->alpha2 = 0.50;
625   tr->alpha3 = 1.00;
626   tr->alpha4 = 2.00;
627   tr->alpha5 = 4.00;
628 
629   /*  Interpolation trust region update parameters */
630   tr->mu1 = 0.10;
631   tr->mu2 = 0.50;
632 
633   tr->gamma1 = 0.25;
634   tr->gamma2 = 0.50;
635   tr->gamma3 = 2.00;
636   tr->gamma4 = 4.00;
637 
638   tr->theta = 0.05;
639 
640   /*  Interpolation parameters for initialization */
641   tr->mu1_i = 0.35;
642   tr->mu2_i = 0.50;
643 
644   tr->gamma1_i = 0.0625;
645   tr->gamma2_i = 0.50;
646   tr->gamma3_i = 2.00;
647   tr->gamma4_i = 5.00;
648 
649   tr->theta_i = 0.25;
650 
651   tr->min_radius = 1.0e-10;
652   tr->max_radius = 1.0e10;
653   tr->epsilon    = 1.0e-6;
654 
655   tr->init_type       = NTR_INIT_INTERPOLATION;
656   tr->update_type     = NTR_UPDATE_REDUCTION;
657 
658   /* Set linear solver to default for trust region */
659   ierr = KSPCreate(((PetscObject)tao)->comm,&tao->ksp);CHKERRQ(ierr);
660   ierr = PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);CHKERRQ(ierr);
661   ierr = KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);CHKERRQ(ierr);
662   ierr = KSPSetType(tao->ksp,KSPCGSTCG);CHKERRQ(ierr);
663   PetscFunctionReturn(0);
664 }
665