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