xref: /petsc/src/tao/unconstrained/impls/ntr/ntr.c (revision bebe2cf65d55febe21a5af8db2bd2e168caaa2e7) 
1 #include <../src/tao/matrix/lmvmmat.h>
2 #include <../src/tao/unconstrained/impls/ntr/ntr.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 NTR_KSP_NASH    0
10 #define NTR_KSP_STCG    1
11 #define NTR_KSP_GLTR    2
12 #define NTR_KSP_TYPES   3
13 
14 #define NTR_PC_NONE     0
15 #define NTR_PC_AHESS    1
16 #define NTR_PC_BFGS     2
17 #define NTR_PC_PETSC    3
18 #define NTR_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 NTR_INIT_CONSTANT         0
25 #define NTR_INIT_DIRECTION        1
26 #define NTR_INIT_INTERPOLATION    2
27 #define NTR_INIT_TYPES            3
28 
29 #define NTR_UPDATE_REDUCTION      0
30 #define NTR_UPDATE_INTERPOLATION  1
31 #define NTR_UPDATE_TYPES          2
32 
33 static const char *NTR_KSP[64] = {  "nash", "stcg", "gltr"};
34 
35 static const char *NTR_PC[64] = {  "none", "ahess", "bfgs", "petsc"};
36 
37 static const char *BFGS_SCALE[64] = {  "ahess", "bfgs"};
38 
39 static const char *NTR_INIT[64] = {  "constant", "direction", "interpolation"};
40 
41 static const char *NTR_UPDATE[64] = {  "reduction", "interpolation"};
42 
43 /*  Routine for BFGS preconditioner */
44 static PetscErrorCode MatLMVMSolveShell(PC pc, Vec xin, Vec xout);
45 
46 /*
47    TaoSolve_NTR - Implements Newton's Method with a trust region approach
48    for solving unconstrained minimization problems.
49 
50    The basic algorithm is taken from MINPACK-2 (dstrn).
51 
52    TaoSolve_NTR computes a local minimizer of a twice differentiable function
53    f by applying a trust region variant of Newton's method.  At each stage
54    of the algorithm, we use the prconditioned conjugate gradient method to
55    determine an approximate minimizer of the quadratic equation
56 
57         q(s) = <s, Hs + g>
58 
59    subject to the trust region constraint
60 
61         || s ||_M <= radius,
62 
63    where radius is the trust region radius and M is a symmetric positive
64    definite matrix (the preconditioner).  Here g is the gradient and H
65    is the Hessian matrix.
66 
67    Note:  TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG,
68           or KSPGLTR.  Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this
69           routine regardless of what the user may have previously specified.
70 */
71 #undef __FUNCT__
72 #define __FUNCT__ "TaoSolve_NTR"
73 static PetscErrorCode TaoSolve_NTR(Tao tao)
74 {
75   TAO_NTR            *tr = (TAO_NTR *)tao->data;
76   PC                 pc;
77   KSPConvergedReason ksp_reason;
78   TaoConvergedReason reason;
79   PetscReal          fmin, ftrial, prered, actred, kappa, sigma, beta;
80   PetscReal          tau, tau_1, tau_2, tau_max, tau_min, max_radius;
81   PetscReal          f, gnorm;
82 
83   PetscReal          delta;
84   PetscReal          norm_d;
85   PetscErrorCode     ierr;
86   PetscInt           bfgsUpdates = 0;
87   PetscInt           needH;
88 
89   PetscInt           i_max = 5;
90   PetscInt           j_max = 1;
91   PetscInt           i, j, N, n, its;
92 
93   PetscFunctionBegin;
94   if (tao->XL || tao->XU || tao->ops->computebounds) {
95     ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");CHKERRQ(ierr);
96   }
97 
98   tao->trust = tao->trust0;
99 
100   /* Modify the radius if it is too large or small */
101   tao->trust = PetscMax(tao->trust, tr->min_radius);
102   tao->trust = PetscMin(tao->trust, tr->max_radius);
103 
104 
105   if (NTR_PC_BFGS == tr->pc_type && !tr->M) {
106     ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr);
107     ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);
108     ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tr->M);CHKERRQ(ierr);
109     ierr = MatLMVMAllocateVectors(tr->M,tao->solution);CHKERRQ(ierr);
110   }
111 
112   /* Check convergence criteria */
113   ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
114   ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
115   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
116   needH = 1;
117 
118   ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr);
119   if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
120 
121   /* Create vectors for the limited memory preconditioner */
122   if ((NTR_PC_BFGS == tr->pc_type) &&
123       (BFGS_SCALE_BFGS != tr->bfgs_scale_type)) {
124     if (!tr->Diag) {
125         ierr = VecDuplicate(tao->solution, &tr->Diag);CHKERRQ(ierr);
126     }
127   }
128 
129   switch(tr->ksp_type) {
130   case NTR_KSP_NASH:
131     ierr = KSPSetType(tao->ksp, KSPNASH);CHKERRQ(ierr);
132     ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
133     break;
134 
135   case NTR_KSP_STCG:
136     ierr = KSPSetType(tao->ksp, KSPSTCG);CHKERRQ(ierr);
137     ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
138     break;
139 
140   default:
141     ierr = KSPSetType(tao->ksp, KSPGLTR);CHKERRQ(ierr);
142     ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
143     break;
144   }
145 
146   /*  Modify the preconditioner to use the bfgs approximation */
147   ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr);
148   switch(tr->pc_type) {
149   case NTR_PC_NONE:
150     ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr);
151     ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
152     break;
153 
154   case NTR_PC_AHESS:
155     ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr);
156     ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
157     ierr = PCJacobiSetUseAbs(pc,PETSC_TRUE);CHKERRQ(ierr);
158     break;
159 
160   case NTR_PC_BFGS:
161     ierr = PCSetType(pc, PCSHELL);CHKERRQ(ierr);
162     ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
163     ierr = PCShellSetName(pc, "bfgs");CHKERRQ(ierr);
164     ierr = PCShellSetContext(pc, tr->M);CHKERRQ(ierr);
165     ierr = PCShellSetApply(pc, MatLMVMSolveShell);CHKERRQ(ierr);
166     break;
167 
168   default:
169     /*  Use the pc method set by pc_type */
170     break;
171   }
172 
173   /*  Initialize trust-region radius */
174   switch(tr->init_type) {
175   case NTR_INIT_CONSTANT:
176     /*  Use the initial radius specified */
177     break;
178 
179   case NTR_INIT_INTERPOLATION:
180     /*  Use the initial radius specified */
181     max_radius = 0.0;
182 
183     for (j = 0; j < j_max; ++j) {
184       fmin = f;
185       sigma = 0.0;
186 
187       if (needH) {
188         ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
189         needH = 0;
190       }
191 
192       for (i = 0; i < i_max; ++i) {
193 
194         ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr);
195         ierr = VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);CHKERRQ(ierr);
196         ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);
197 
198         if (PetscIsInfOrNanReal(ftrial)) {
199           tau = tr->gamma1_i;
200         }
201         else {
202           if (ftrial < fmin) {
203             fmin = ftrial;
204             sigma = -tao->trust / gnorm;
205           }
206 
207           ierr = MatMult(tao->hessian, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
208           ierr = VecDot(tao->gradient, tao->stepdirection, &prered);CHKERRQ(ierr);
209 
210           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
211           actred = f - ftrial;
212           if ((PetscAbsScalar(actred) <= tr->epsilon) &&
213               (PetscAbsScalar(prered) <= tr->epsilon)) {
214             kappa = 1.0;
215           }
216           else {
217             kappa = actred / prered;
218           }
219 
220           tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
221           tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
222           tau_min = PetscMin(tau_1, tau_2);
223           tau_max = PetscMax(tau_1, tau_2);
224 
225           if (PetscAbsScalar(kappa - 1.0) <= tr->mu1_i) {
226             /*  Great agreement */
227             max_radius = PetscMax(max_radius, tao->trust);
228 
229             if (tau_max < 1.0) {
230               tau = tr->gamma3_i;
231             }
232             else if (tau_max > tr->gamma4_i) {
233               tau = tr->gamma4_i;
234             }
235             else {
236               tau = tau_max;
237             }
238           }
239           else if (PetscAbsScalar(kappa - 1.0) <= tr->mu2_i) {
240             /*  Good agreement */
241             max_radius = PetscMax(max_radius, tao->trust);
242 
243             if (tau_max < tr->gamma2_i) {
244               tau = tr->gamma2_i;
245             }
246             else if (tau_max > tr->gamma3_i) {
247               tau = tr->gamma3_i;
248             }
249             else {
250               tau = tau_max;
251             }
252           }
253           else {
254             /*  Not good agreement */
255             if (tau_min > 1.0) {
256               tau = tr->gamma2_i;
257             }
258             else if (tau_max < tr->gamma1_i) {
259               tau = tr->gamma1_i;
260             }
261             else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
262               tau = tr->gamma1_i;
263             }
264             else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) &&
265                      ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
266               tau = tau_1;
267             }
268             else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) &&
269                      ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
270               tau = tau_2;
271             }
272             else {
273               tau = tau_max;
274             }
275           }
276         }
277         tao->trust = tau * tao->trust;
278       }
279 
280       if (fmin < f) {
281         f = fmin;
282         ierr = VecAXPY(tao->solution, sigma, tao->gradient);CHKERRQ(ierr);
283         ierr = TaoComputeGradient(tao,tao->solution, tao->gradient);CHKERRQ(ierr);
284 
285         ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr);
286 
287         if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
288         needH = 1;
289 
290         ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr);
291         if (reason != TAO_CONTINUE_ITERATING) {
292           PetscFunctionReturn(0);
293         }
294       }
295     }
296     tao->trust = PetscMax(tao->trust, max_radius);
297 
298     /*  Modify the radius if it is too large or small */
299     tao->trust = PetscMax(tao->trust, tr->min_radius);
300     tao->trust = PetscMin(tao->trust, tr->max_radius);
301     break;
302 
303   default:
304     /*  Norm of the first direction will initialize radius */
305     tao->trust = 0.0;
306     break;
307   }
308 
309   /* Set initial scaling for the BFGS preconditioner
310      This step is done after computing the initial trust-region radius
311      since the function value may have decreased */
312   if (NTR_PC_BFGS == tr->pc_type) {
313     if (f != 0.0) {
314       delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
315     }
316     else {
317       delta = 2.0 / (gnorm*gnorm);
318     }
319     ierr = MatLMVMSetDelta(tr->M,delta);CHKERRQ(ierr);
320   }
321 
322   /* Have not converged; continue with Newton method */
323   while (reason == TAO_CONTINUE_ITERATING) {
324     ++tao->niter;
325     tao->ksp_its=0;
326     /* Compute the Hessian */
327     if (needH) {
328       ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
329       needH = 0;
330     }
331 
332     if (NTR_PC_BFGS == tr->pc_type) {
333       if (BFGS_SCALE_AHESS == tr->bfgs_scale_type) {
334         /* Obtain diagonal for the bfgs preconditioner */
335         ierr = MatGetDiagonal(tao->hessian, tr->Diag);CHKERRQ(ierr);
336         ierr = VecAbs(tr->Diag);CHKERRQ(ierr);
337         ierr = VecReciprocal(tr->Diag);CHKERRQ(ierr);
338         ierr = MatLMVMSetScale(tr->M,tr->Diag);CHKERRQ(ierr);
339       }
340 
341       /* Update the limited memory preconditioner */
342       ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr);
343       ++bfgsUpdates;
344     }
345 
346     while (reason == TAO_CONTINUE_ITERATING) {
347       ierr = KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);CHKERRQ(ierr);
348 
349       /* Solve the trust region subproblem */
350       if (NTR_KSP_NASH == tr->ksp_type) {
351         ierr = KSPNASHSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
352         ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
353         ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
354         tao->ksp_its+=its;
355         tao->ksp_tot_its+=its;
356         ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
357       } else if (NTR_KSP_STCG == tr->ksp_type) {
358         ierr = KSPSTCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
359         ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
360         ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
361         tao->ksp_its+=its;
362         tao->ksp_tot_its+=its;
363         ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
364       } else { /* NTR_KSP_GLTR */
365         ierr = KSPGLTRSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
366         ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
367         ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
368         tao->ksp_its+=its;
369         tao->ksp_tot_its+=its;
370         ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
371       }
372 
373       if (0.0 == tao->trust) {
374         /* Radius was uninitialized; use the norm of the direction */
375         if (norm_d > 0.0) {
376           tao->trust = norm_d;
377 
378           /* Modify the radius if it is too large or small */
379           tao->trust = PetscMax(tao->trust, tr->min_radius);
380           tao->trust = PetscMin(tao->trust, tr->max_radius);
381         }
382         else {
383           /* The direction was bad; set radius to default value and re-solve
384              the trust-region subproblem to get a direction */
385           tao->trust = tao->trust0;
386 
387           /* Modify the radius if it is too large or small */
388           tao->trust = PetscMax(tao->trust, tr->min_radius);
389           tao->trust = PetscMin(tao->trust, tr->max_radius);
390 
391           if (NTR_KSP_NASH == tr->ksp_type) {
392             ierr = KSPNASHSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
393             ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
394             ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
395             tao->ksp_its+=its;
396             tao->ksp_tot_its+=its;
397             ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
398           } else if (NTR_KSP_STCG == tr->ksp_type) {
399             ierr = KSPSTCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
400             ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
401             ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
402             tao->ksp_its+=its;
403             tao->ksp_tot_its+=its;
404             ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
405           } else { /* NTR_KSP_GLTR */
406             ierr = KSPGLTRSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
407             ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
408             ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
409             tao->ksp_its+=its;
410             tao->ksp_tot_its+=its;
411             ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
412           }
413 
414           if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
415         }
416       }
417       ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
418       ierr = KSPGetConvergedReason(tao->ksp, &ksp_reason);CHKERRQ(ierr);
419       if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
420           (NTR_PC_BFGS == tr->pc_type) && (bfgsUpdates > 1)) {
421         /* Preconditioner is numerically indefinite; reset the
422            approximate if using BFGS preconditioning. */
423 
424         if (f != 0.0) {
425           delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
426         }
427         else {
428           delta = 2.0 / (gnorm*gnorm);
429         }
430         ierr = MatLMVMSetDelta(tr->M, delta);CHKERRQ(ierr);
431         ierr = MatLMVMReset(tr->M);CHKERRQ(ierr);
432         ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr);
433         bfgsUpdates = 1;
434       }
435 
436       if (NTR_UPDATE_REDUCTION == tr->update_type) {
437         /* Get predicted reduction */
438         if (NTR_KSP_NASH == tr->ksp_type) {
439           ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
440         } else if (NTR_KSP_STCG == tr->ksp_type) {
441           ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
442         } else { /* gltr */
443           ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
444         }
445 
446         if (prered >= 0.0) {
447           /* The predicted reduction has the wrong sign.  This cannot
448              happen in infinite precision arithmetic.  Step should
449              be rejected! */
450           tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
451         }
452         else {
453           /* Compute trial step and function value */
454           ierr = VecCopy(tao->solution,tr->W);CHKERRQ(ierr);
455           ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr);
456           ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);
457 
458           if (PetscIsInfOrNanReal(ftrial)) {
459             tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
460           } else {
461             /* Compute and actual reduction */
462             actred = f - ftrial;
463             prered = -prered;
464             if ((PetscAbsScalar(actred) <= tr->epsilon) &&
465                 (PetscAbsScalar(prered) <= tr->epsilon)) {
466               kappa = 1.0;
467             }
468             else {
469               kappa = actred / prered;
470             }
471 
472             /* Accept or reject the step and update radius */
473             if (kappa < tr->eta1) {
474               /* Reject the step */
475               tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
476             }
477             else {
478               /* Accept the step */
479               if (kappa < tr->eta2) {
480                 /* Marginal bad step */
481                 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
482               }
483               else if (kappa < tr->eta3) {
484                 /* Reasonable step */
485                 tao->trust = tr->alpha3 * tao->trust;
486               }
487               else if (kappa < tr->eta4) {
488                 /* Good step */
489                 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
490               }
491               else {
492                 /* Very good step */
493                 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
494               }
495               break;
496             }
497           }
498         }
499       }
500       else {
501         /* Get predicted reduction */
502         if (NTR_KSP_NASH == tr->ksp_type) {
503           ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
504         } else if (NTR_KSP_STCG == tr->ksp_type) {
505           ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
506         } else { /* gltr */
507           ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
508         }
509 
510         if (prered >= 0.0) {
511           /* The predicted reduction has the wrong sign.  This cannot
512              happen in infinite precision arithmetic.  Step should
513              be rejected! */
514           tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
515         }
516         else {
517           ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr);
518           ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr);
519           ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);
520           if (PetscIsInfOrNanReal(ftrial)) {
521             tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
522           }
523           else {
524             ierr = VecDot(tao->gradient, tao->stepdirection, &beta);CHKERRQ(ierr);
525             actred = f - ftrial;
526             prered = -prered;
527             if ((PetscAbsScalar(actred) <= tr->epsilon) &&
528                 (PetscAbsScalar(prered) <= tr->epsilon)) {
529               kappa = 1.0;
530             }
531             else {
532               kappa = actred / prered;
533             }
534 
535             tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
536             tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
537             tau_min = PetscMin(tau_1, tau_2);
538             tau_max = PetscMax(tau_1, tau_2);
539 
540             if (kappa >= 1.0 - tr->mu1) {
541               /* Great agreement; accept step and update radius */
542               if (tau_max < 1.0) {
543                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
544               }
545               else if (tau_max > tr->gamma4) {
546                 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
547               }
548               else {
549                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
550               }
551               break;
552             }
553             else if (kappa >= 1.0 - tr->mu2) {
554               /* Good agreement */
555 
556               if (tau_max < tr->gamma2) {
557                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
558               }
559               else if (tau_max > tr->gamma3) {
560                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
561               }
562               else if (tau_max < 1.0) {
563                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
564               }
565               else {
566                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
567               }
568               break;
569             }
570             else {
571               /* Not good agreement */
572               if (tau_min > 1.0) {
573                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
574               }
575               else if (tau_max < tr->gamma1) {
576                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
577               }
578               else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
579                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
580               }
581               else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) &&
582                        ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
583                 tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
584               }
585               else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) &&
586                        ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
587                 tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
588               }
589               else {
590                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
591               }
592             }
593           }
594         }
595       }
596 
597       /* The step computed was not good and the radius was decreased.
598          Monitor the radius to terminate. */
599       ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust, &reason);CHKERRQ(ierr);
600     }
601 
602     /* The radius may have been increased; modify if it is too large */
603     tao->trust = PetscMin(tao->trust, tr->max_radius);
604 
605     if (reason == TAO_CONTINUE_ITERATING) {
606       ierr = VecCopy(tr->W, tao->solution);CHKERRQ(ierr);
607       f = ftrial;
608       ierr = TaoComputeGradient(tao, tao->solution, tao->gradient);CHKERRQ(ierr);
609       ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr);
610       if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
611       needH = 1;
612       ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust, &reason);CHKERRQ(ierr);
613     }
614   }
615   PetscFunctionReturn(0);
616 }
617 
618 /*------------------------------------------------------------*/
619 #undef __FUNCT__
620 #define __FUNCT__ "TaoSetUp_NTR"
621 static PetscErrorCode TaoSetUp_NTR(Tao tao)
622 {
623   TAO_NTR *tr = (TAO_NTR *)tao->data;
624   PetscErrorCode ierr;
625 
626   PetscFunctionBegin;
627 
628   if (!tao->gradient) {ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr);}
629   if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);}
630   if (!tr->W) {ierr = VecDuplicate(tao->solution, &tr->W);CHKERRQ(ierr);}
631 
632   tr->Diag = 0;
633   tr->M = 0;
634 
635 
636   PetscFunctionReturn(0);
637 }
638 
639 /*------------------------------------------------------------*/
640 #undef __FUNCT__
641 #define __FUNCT__ "TaoDestroy_NTR"
642 static PetscErrorCode TaoDestroy_NTR(Tao tao)
643 {
644   TAO_NTR        *tr = (TAO_NTR *)tao->data;
645   PetscErrorCode ierr;
646 
647   PetscFunctionBegin;
648   if (tao->setupcalled) {
649     ierr = VecDestroy(&tr->W);CHKERRQ(ierr);
650   }
651   ierr = MatDestroy(&tr->M);CHKERRQ(ierr);
652   ierr = VecDestroy(&tr->Diag);CHKERRQ(ierr);
653   ierr = PetscFree(tao->data);CHKERRQ(ierr);
654   PetscFunctionReturn(0);
655 }
656 
657 /*------------------------------------------------------------*/
658 #undef __FUNCT__
659 #define __FUNCT__ "TaoSetFromOptions_NTR"
660 static PetscErrorCode TaoSetFromOptions_NTR(PetscOptions *PetscOptionsObject,Tao tao)
661 {
662   TAO_NTR        *tr = (TAO_NTR *)tao->data;
663   PetscErrorCode ierr;
664 
665   PetscFunctionBegin;
666   ierr = PetscOptionsHead(PetscOptionsObject,"Newton trust region method for unconstrained optimization");CHKERRQ(ierr);
667   ierr = PetscOptionsEList("-tao_ntr_ksp_type", "ksp type", "", NTR_KSP, NTR_KSP_TYPES, NTR_KSP[tr->ksp_type], &tr->ksp_type,NULL);CHKERRQ(ierr);
668   ierr = PetscOptionsEList("-tao_ntr_pc_type", "pc type", "", NTR_PC, NTR_PC_TYPES, NTR_PC[tr->pc_type], &tr->pc_type,NULL);CHKERRQ(ierr);
669   ierr = PetscOptionsEList("-tao_ntr_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tr->bfgs_scale_type], &tr->bfgs_scale_type,NULL);CHKERRQ(ierr);
670   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);
671   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);
672   ierr = PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL);CHKERRQ(ierr);
673   ierr = PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL);CHKERRQ(ierr);
674   ierr = PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL);CHKERRQ(ierr);
675   ierr = PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL);CHKERRQ(ierr);
676   ierr = PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL);CHKERRQ(ierr);
677   ierr = PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL);CHKERRQ(ierr);
678   ierr = PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL);CHKERRQ(ierr);
679   ierr = PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL);CHKERRQ(ierr);
680   ierr = PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL);CHKERRQ(ierr);
681   ierr = PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL);CHKERRQ(ierr);
682   ierr = PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL);CHKERRQ(ierr);
683   ierr = PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL);CHKERRQ(ierr);
684   ierr = PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL);CHKERRQ(ierr);
685   ierr = PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL);CHKERRQ(ierr);
686   ierr = PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL);CHKERRQ(ierr);
687   ierr = PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL);CHKERRQ(ierr);
688   ierr = PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL);CHKERRQ(ierr);
689   ierr = PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL);CHKERRQ(ierr);
690   ierr = PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL);CHKERRQ(ierr);
691   ierr = PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL);CHKERRQ(ierr);
692   ierr = PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL);CHKERRQ(ierr);
693   ierr = PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL);CHKERRQ(ierr);
694   ierr = PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL);CHKERRQ(ierr);
695   ierr = PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL);CHKERRQ(ierr);
696   ierr = PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL);CHKERRQ(ierr);
697   ierr = PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL);CHKERRQ(ierr);
698   ierr = PetscOptionsTail();CHKERRQ(ierr);
699   ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
700   PetscFunctionReturn(0);
701 }
702 
703 /*------------------------------------------------------------*/
704 #undef __FUNCT__
705 #define __FUNCT__ "TaoView_NTR"
706 static PetscErrorCode TaoView_NTR(Tao tao, PetscViewer viewer)
707 {
708   TAO_NTR        *tr = (TAO_NTR *)tao->data;
709   PetscErrorCode ierr;
710   PetscInt       nrejects;
711   PetscBool      isascii;
712 
713   PetscFunctionBegin;
714   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
715   if (isascii) {
716     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
717     if (NTR_PC_BFGS == tr->pc_type && tr->M) {
718       ierr = MatLMVMGetRejects(tr->M, &nrejects);CHKERRQ(ierr);
719       ierr = PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n", nrejects);CHKERRQ(ierr);
720     }
721     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
722   }
723   PetscFunctionReturn(0);
724 }
725 
726 /*------------------------------------------------------------*/
727 /*MC
728   TAONTR - Newton's method with trust region for unconstrained minimization.
729   At each iteration, the Newton trust region method solves the system.
730   NTR expects a KSP solver with a trust region radius.
731             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k
732 
733   Options Database Keys:
734 + -tao_ntr_ksp_type - "nash","stcg","gltr"
735 . -tao_ntr_pc_type - "none","ahess","bfgs","petsc"
736 . -tao_ntr_bfgs_scale_type - type of scaling with bfgs pc, "ahess" or "bfgs"
737 . -tao_ntr_init_type - "constant","direction","interpolation"
738 . -tao_ntr_update_type - "reduction","interpolation"
739 . -tao_ntr_min_radius - lower bound on trust region radius
740 . -tao_ntr_max_radius - upper bound on trust region radius
741 . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
742 . -tao_ntr_mu1_i - mu1 interpolation init factor
743 . -tao_ntr_mu2_i - mu2 interpolation init factor
744 . -tao_ntr_gamma1_i - gamma1 interpolation init factor
745 . -tao_ntr_gamma2_i - gamma2 interpolation init factor
746 . -tao_ntr_gamma3_i - gamma3 interpolation init factor
747 . -tao_ntr_gamma4_i - gamma4 interpolation init factor
748 . -tao_ntr_theta_i - thetha1 interpolation init factor
749 . -tao_ntr_eta1 - eta1 reduction update factor
750 . -tao_ntr_eta2 - eta2 reduction update factor
751 . -tao_ntr_eta3 - eta3 reduction update factor
752 . -tao_ntr_eta4 - eta4 reduction update factor
753 . -tao_ntr_alpha1 - alpha1 reduction update factor
754 . -tao_ntr_alpha2 - alpha2 reduction update factor
755 . -tao_ntr_alpha3 - alpha3 reduction update factor
756 . -tao_ntr_alpha4 - alpha4 reduction update factor
757 . -tao_ntr_alpha4 - alpha4 reduction update factor
758 . -tao_ntr_mu1 - mu1 interpolation update
759 . -tao_ntr_mu2 - mu2 interpolation update
760 . -tao_ntr_gamma1 - gamma1 interpolcation update
761 . -tao_ntr_gamma2 - gamma2 interpolcation update
762 . -tao_ntr_gamma3 - gamma3 interpolcation update
763 . -tao_ntr_gamma4 - gamma4 interpolation update
764 - -tao_ntr_theta - theta interpolation update
765 
766   Level: beginner
767 M*/
768 
769 #undef __FUNCT__
770 #define __FUNCT__ "TaoCreate_NTR"
771 PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
772 {
773   TAO_NTR *tr;
774   PetscErrorCode ierr;
775 
776   PetscFunctionBegin;
777 
778   ierr = PetscNewLog(tao,&tr);CHKERRQ(ierr);
779 
780   tao->ops->setup = TaoSetUp_NTR;
781   tao->ops->solve = TaoSolve_NTR;
782   tao->ops->view = TaoView_NTR;
783   tao->ops->setfromoptions = TaoSetFromOptions_NTR;
784   tao->ops->destroy = TaoDestroy_NTR;
785 
786   /* Override default settings (unless already changed) */
787   if (!tao->max_it_changed) tao->max_it = 50;
788   if (!tao->trust0_changed) tao->trust0 = 100.0;
789 #if defined(PETSC_USE_REAL_SINGLE)
790   if (!tao->fatol_changed) tao->fatol = 1.0e-5;
791   if (!tao->frtol_changed) tao->frtol = 1.0e-5;
792 #else
793   if (!tao->fatol_changed) tao->fatol = 1.0e-10;
794   if (!tao->frtol_changed) tao->frtol = 1.0e-10;
795 #endif
796   tao->data = (void*)tr;
797 
798   /*  Standard trust region update parameters */
799   tr->eta1 = 1.0e-4;
800   tr->eta2 = 0.25;
801   tr->eta3 = 0.50;
802   tr->eta4 = 0.90;
803 
804   tr->alpha1 = 0.25;
805   tr->alpha2 = 0.50;
806   tr->alpha3 = 1.00;
807   tr->alpha4 = 2.00;
808   tr->alpha5 = 4.00;
809 
810   /*  Interpolation parameters */
811   tr->mu1_i = 0.35;
812   tr->mu2_i = 0.50;
813 
814   tr->gamma1_i = 0.0625;
815   tr->gamma2_i = 0.50;
816   tr->gamma3_i = 2.00;
817   tr->gamma4_i = 5.00;
818 
819   tr->theta_i = 0.25;
820 
821   /*  Interpolation trust region update parameters */
822   tr->mu1 = 0.10;
823   tr->mu2 = 0.50;
824 
825   tr->gamma1 = 0.25;
826   tr->gamma2 = 0.50;
827   tr->gamma3 = 2.00;
828   tr->gamma4 = 4.00;
829 
830   tr->theta = 0.05;
831 
832   tr->min_radius = 1.0e-10;
833   tr->max_radius = 1.0e10;
834   tr->epsilon = 1.0e-6;
835 
836   tr->ksp_type        = NTR_KSP_STCG;
837   tr->pc_type         = NTR_PC_BFGS;
838   tr->bfgs_scale_type = BFGS_SCALE_AHESS;
839   tr->init_type       = NTR_INIT_INTERPOLATION;
840   tr->update_type     = NTR_UPDATE_REDUCTION;
841 
842 
843   /* Set linear solver to default for trust region */
844   ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr);
845   ierr = KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);CHKERRQ(ierr);
846   PetscFunctionReturn(0);
847 }
848 
849 
850 #undef __FUNCT__
851 #define __FUNCT__ "MatLMVMSolveShell"
852 static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x)
853 {
854     PetscErrorCode ierr;
855     Mat M;
856     PetscFunctionBegin;
857     PetscValidHeaderSpecific(pc,PC_CLASSID,1);
858     PetscValidHeaderSpecific(b,VEC_CLASSID,2);
859     PetscValidHeaderSpecific(x,VEC_CLASSID,3);
860     ierr = PCShellGetContext(pc,(void**)&M);CHKERRQ(ierr);
861     ierr = MatLMVMSolve(M, b, x);CHKERRQ(ierr);
862     PetscFunctionReturn(0);
863 }
864