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