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 needH;
57: PetscInt i_max = 5;
58: PetscInt j_max = 1;
59: PetscInt i, j, N, n, its;
61: if (tao->XL || tao->XU || tao->ops->computebounds) PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");
63: KSPGetType(tao->ksp, &ksp_type);
64: PetscStrcmp(ksp_type, KSPNASH, &is_nash);
65: PetscStrcmp(ksp_type, KSPSTCG, &is_stcg);
66: PetscStrcmp(ksp_type, KSPGLTR, &is_gltr);
69: /* Initialize the radius and modify if it is too large or small */
70: tao->trust = tao->trust0;
71: tao->trust = PetscMax(tao->trust, tr->min_radius);
72: tao->trust = PetscMin(tao->trust, tr->max_radius);
74: /* Allocate the vectors needed for the BFGS approximation */
75: KSPGetPC(tao->ksp, &pc);
76: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);
77: PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);
78: if (is_bfgs) {
79: tr->bfgs_pre = pc;
80: PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M);
81: VecGetLocalSize(tao->solution, &n);
82: VecGetSize(tao->solution, &N);
83: MatSetSizes(tr->M, n, n, N, N);
84: MatLMVMAllocate(tr->M, tao->solution, tao->gradient);
85: MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric);
87: } else if (is_jacobi) PCJacobiSetUseAbs(pc, PETSC_TRUE);
89: /* Check convergence criteria */
90: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
91: TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);
93: needH = 1;
95: tao->reason = TAO_CONTINUE_ITERATING;
96: TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
97: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0);
98: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
99: if (tao->reason != TAO_CONTINUE_ITERATING) return 0;
101: /* Initialize trust-region radius */
102: switch (tr->init_type) {
103: case NTR_INIT_CONSTANT:
104: /* Use the initial radius specified */
105: break;
107: case NTR_INIT_INTERPOLATION:
108: /* Use the initial radius specified */
109: max_radius = 0.0;
111: for (j = 0; j < j_max; ++j) {
112: fmin = f;
113: sigma = 0.0;
115: if (needH) {
116: TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);
117: needH = 0;
118: }
120: for (i = 0; i < i_max; ++i) {
121: VecCopy(tao->solution, tr->W);
122: VecAXPY(tr->W, -tao->trust / gnorm, tao->gradient);
123: TaoComputeObjective(tao, tr->W, &ftrial);
125: if (PetscIsInfOrNanReal(ftrial)) {
126: tau = tr->gamma1_i;
127: } else {
128: if (ftrial < fmin) {
129: fmin = ftrial;
130: sigma = -tao->trust / gnorm;
131: }
133: MatMult(tao->hessian, tao->gradient, tao->stepdirection);
134: VecDot(tao->gradient, tao->stepdirection, &prered);
136: prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
137: actred = f - ftrial;
138: if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
139: kappa = 1.0;
140: } else {
141: kappa = actred / prered;
142: }
144: tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
145: tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
146: tau_min = PetscMin(tau_1, tau_2);
147: tau_max = PetscMax(tau_1, tau_2);
149: if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu1_i) {
150: /* Great agreement */
151: max_radius = PetscMax(max_radius, tao->trust);
153: if (tau_max < 1.0) {
154: tau = tr->gamma3_i;
155: } else if (tau_max > tr->gamma4_i) {
156: tau = tr->gamma4_i;
157: } else {
158: tau = tau_max;
159: }
160: } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu2_i) {
161: /* Good agreement */
162: max_radius = PetscMax(max_radius, tao->trust);
164: if (tau_max < tr->gamma2_i) {
165: tau = tr->gamma2_i;
166: } else if (tau_max > tr->gamma3_i) {
167: tau = tr->gamma3_i;
168: } else {
169: tau = tau_max;
170: }
171: } else {
172: /* Not good agreement */
173: if (tau_min > 1.0) {
174: tau = tr->gamma2_i;
175: } else if (tau_max < tr->gamma1_i) {
176: tau = tr->gamma1_i;
177: } else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
178: tau = tr->gamma1_i;
179: } else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
180: tau = tau_1;
181: } else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
182: tau = tau_2;
183: } else {
184: tau = tau_max;
185: }
186: }
187: }
188: tao->trust = tau * tao->trust;
189: }
191: if (fmin < f) {
192: f = fmin;
193: VecAXPY(tao->solution, sigma, tao->gradient);
194: TaoComputeGradient(tao, tao->solution, tao->gradient);
196: TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);
198: needH = 1;
200: TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
201: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0);
202: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
203: if (tao->reason != TAO_CONTINUE_ITERATING) return 0;
204: }
205: }
206: tao->trust = PetscMax(tao->trust, max_radius);
208: /* Modify the radius if it is too large or small */
209: tao->trust = PetscMax(tao->trust, tr->min_radius);
210: tao->trust = PetscMin(tao->trust, tr->max_radius);
211: break;
213: default:
214: /* Norm of the first direction will initialize radius */
215: tao->trust = 0.0;
216: break;
217: }
219: /* Have not converged; continue with Newton method */
220: while (tao->reason == TAO_CONTINUE_ITERATING) {
221: /* Call general purpose update function */
222: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
223: ++tao->niter;
224: tao->ksp_its = 0;
225: /* Compute the Hessian */
226: if (needH) {
227: TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);
228: needH = 0;
229: }
231: if (tr->bfgs_pre) {
232: /* Update the limited memory preconditioner */
233: MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
234: }
236: while (tao->reason == TAO_CONTINUE_ITERATING) {
237: KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);
239: /* Solve the trust region subproblem */
240: KSPCGSetRadius(tao->ksp, tao->trust);
241: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
242: KSPGetIterationNumber(tao->ksp, &its);
243: tao->ksp_its += its;
244: tao->ksp_tot_its += its;
245: KSPCGGetNormD(tao->ksp, &norm_d);
247: if (0.0 == tao->trust) {
248: /* Radius was uninitialized; use the norm of the direction */
249: if (norm_d > 0.0) {
250: tao->trust = norm_d;
252: /* Modify the radius if it is too large or small */
253: tao->trust = PetscMax(tao->trust, tr->min_radius);
254: tao->trust = PetscMin(tao->trust, tr->max_radius);
255: } else {
256: /* The direction was bad; set radius to default value and re-solve
257: the trust-region subproblem to get a direction */
258: tao->trust = tao->trust0;
260: /* Modify the radius if it is too large or small */
261: tao->trust = PetscMax(tao->trust, tr->min_radius);
262: tao->trust = PetscMin(tao->trust, tr->max_radius);
264: KSPCGSetRadius(tao->ksp, tao->trust);
265: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
266: KSPGetIterationNumber(tao->ksp, &its);
267: tao->ksp_its += its;
268: tao->ksp_tot_its += its;
269: KSPCGGetNormD(tao->ksp, &norm_d);
272: }
273: }
274: VecScale(tao->stepdirection, -1.0);
275: KSPGetConvergedReason(tao->ksp, &ksp_reason);
276: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) {
277: /* Preconditioner is numerically indefinite; reset the
278: approximate if using BFGS preconditioning. */
279: MatLMVMReset(tr->M, PETSC_FALSE);
280: MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
281: }
283: if (NTR_UPDATE_REDUCTION == tr->update_type) {
284: /* Get predicted reduction */
285: KSPCGGetObjFcn(tao->ksp, &prered);
286: if (prered >= 0.0) {
287: /* The predicted reduction has the wrong sign. This cannot
288: happen in infinite precision arithmetic. Step should
289: be rejected! */
290: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
291: } else {
292: /* Compute trial step and function value */
293: VecCopy(tao->solution, tr->W);
294: VecAXPY(tr->W, 1.0, tao->stepdirection);
295: TaoComputeObjective(tao, tr->W, &ftrial);
297: if (PetscIsInfOrNanReal(ftrial)) {
298: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
299: } else {
300: /* Compute and actual reduction */
301: actred = f - ftrial;
302: prered = -prered;
303: if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
304: kappa = 1.0;
305: } else {
306: kappa = actred / prered;
307: }
309: /* Accept or reject the step and update radius */
310: if (kappa < tr->eta1) {
311: /* Reject the step */
312: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
313: } else {
314: /* Accept the step */
315: if (kappa < tr->eta2) {
316: /* Marginal bad step */
317: tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
318: } else if (kappa < tr->eta3) {
319: /* Reasonable step */
320: tao->trust = tr->alpha3 * tao->trust;
321: } else if (kappa < tr->eta4) {
322: /* Good step */
323: tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
324: } else {
325: /* Very good step */
326: tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
327: }
328: break;
329: }
330: }
331: }
332: } else {
333: /* Get predicted reduction */
334: KSPCGGetObjFcn(tao->ksp, &prered);
335: if (prered >= 0.0) {
336: /* The predicted reduction has the wrong sign. This cannot
337: happen in infinite precision arithmetic. Step should
338: be rejected! */
339: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
340: } else {
341: VecCopy(tao->solution, tr->W);
342: VecAXPY(tr->W, 1.0, tao->stepdirection);
343: TaoComputeObjective(tao, tr->W, &ftrial);
344: if (PetscIsInfOrNanReal(ftrial)) {
345: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
346: } else {
347: VecDot(tao->gradient, tao->stepdirection, &beta);
348: actred = f - ftrial;
349: prered = -prered;
350: if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
351: kappa = 1.0;
352: } else {
353: kappa = actred / prered;
354: }
356: tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
357: tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
358: tau_min = PetscMin(tau_1, tau_2);
359: tau_max = PetscMax(tau_1, tau_2);
361: if (kappa >= 1.0 - tr->mu1) {
362: /* Great agreement; accept step and update radius */
363: if (tau_max < 1.0) {
364: tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
365: } else if (tau_max > tr->gamma4) {
366: tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
367: } else {
368: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
369: }
370: break;
371: } else if (kappa >= 1.0 - tr->mu2) {
372: /* Good agreement */
374: if (tau_max < tr->gamma2) {
375: tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
376: } else if (tau_max > tr->gamma3) {
377: tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
378: } else if (tau_max < 1.0) {
379: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
380: } else {
381: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
382: }
383: break;
384: } else {
385: /* Not good agreement */
386: if (tau_min > 1.0) {
387: tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
388: } else if (tau_max < tr->gamma1) {
389: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
390: } else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
391: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
392: } else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) && ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
393: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
394: } else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) && ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
395: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
396: } else {
397: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
398: }
399: }
400: }
401: }
402: }
404: /* The step computed was not good and the radius was decreased.
405: Monitor the radius to terminate. */
406: TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
407: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust);
408: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
409: }
411: /* The radius may have been increased; modify if it is too large */
412: tao->trust = PetscMin(tao->trust, tr->max_radius);
414: if (tao->reason == TAO_CONTINUE_ITERATING) {
415: VecCopy(tr->W, tao->solution);
416: f = ftrial;
417: TaoComputeGradient(tao, tao->solution, tao->gradient);
418: TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);
420: needH = 1;
421: TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
422: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust);
423: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
424: }
425: }
426: return 0;
427: }
429: /*------------------------------------------------------------*/
430: static PetscErrorCode TaoSetUp_NTR(Tao tao)
431: {
432: TAO_NTR *tr = (TAO_NTR *)tao->data;
434: if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
435: if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
436: if (!tr->W) VecDuplicate(tao->solution, &tr->W);
438: tr->bfgs_pre = NULL;
439: tr->M = NULL;
440: return 0;
441: }
443: /*------------------------------------------------------------*/
444: static PetscErrorCode TaoDestroy_NTR(Tao tao)
445: {
446: TAO_NTR *tr = (TAO_NTR *)tao->data;
448: if (tao->setupcalled) VecDestroy(&tr->W);
449: KSPDestroy(&tao->ksp);
450: PetscFree(tao->data);
451: return 0;
452: }
454: /*------------------------------------------------------------*/
455: static PetscErrorCode TaoSetFromOptions_NTR(Tao tao, PetscOptionItems *PetscOptionsObject)
456: {
457: TAO_NTR *tr = (TAO_NTR *)tao->data;
459: PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region method for unconstrained optimization");
460: PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type, NULL);
461: PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type, NULL);
462: PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1, NULL);
463: PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2, NULL);
464: PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3, NULL);
465: PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4, NULL);
466: PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1, NULL);
467: PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2, NULL);
468: PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3, NULL);
469: PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4, NULL);
470: PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5, NULL);
471: PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1, NULL);
472: PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2, NULL);
473: PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1, NULL);
474: PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2, NULL);
475: PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3, NULL);
476: PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4, NULL);
477: PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta, NULL);
478: PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i, NULL);
479: PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i, NULL);
480: PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i, NULL);
481: PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i, NULL);
482: PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i, NULL);
483: PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i, NULL);
484: PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i, NULL);
485: PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius, NULL);
486: PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius, NULL);
487: PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon, NULL);
488: PetscOptionsHeadEnd();
489: KSPSetFromOptions(tao->ksp);
490: return 0;
491: }
493: /*------------------------------------------------------------*/
494: /*MC
495: TAONTR - Newton's method with trust region for unconstrained minimization.
496: At each iteration, the Newton trust region method solves the system.
497: NTR expects a KSP solver with a trust region radius.
498: min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k
500: Options Database Keys:
501: + -tao_ntr_init_type - "constant","direction","interpolation"
502: . -tao_ntr_update_type - "reduction","interpolation"
503: . -tao_ntr_min_radius - lower bound on trust region radius
504: . -tao_ntr_max_radius - upper bound on trust region radius
505: . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
506: . -tao_ntr_mu1_i - mu1 interpolation init factor
507: . -tao_ntr_mu2_i - mu2 interpolation init factor
508: . -tao_ntr_gamma1_i - gamma1 interpolation init factor
509: . -tao_ntr_gamma2_i - gamma2 interpolation init factor
510: . -tao_ntr_gamma3_i - gamma3 interpolation init factor
511: . -tao_ntr_gamma4_i - gamma4 interpolation init factor
512: . -tao_ntr_theta_i - theta1 interpolation init factor
513: . -tao_ntr_eta1 - eta1 reduction update factor
514: . -tao_ntr_eta2 - eta2 reduction update factor
515: . -tao_ntr_eta3 - eta3 reduction update factor
516: . -tao_ntr_eta4 - eta4 reduction update factor
517: . -tao_ntr_alpha1 - alpha1 reduction update factor
518: . -tao_ntr_alpha2 - alpha2 reduction update factor
519: . -tao_ntr_alpha3 - alpha3 reduction update factor
520: . -tao_ntr_alpha4 - alpha4 reduction update factor
521: . -tao_ntr_alpha4 - alpha4 reduction update factor
522: . -tao_ntr_mu1 - mu1 interpolation update
523: . -tao_ntr_mu2 - mu2 interpolation update
524: . -tao_ntr_gamma1 - gamma1 interpolcation update
525: . -tao_ntr_gamma2 - gamma2 interpolcation update
526: . -tao_ntr_gamma3 - gamma3 interpolcation update
527: . -tao_ntr_gamma4 - gamma4 interpolation update
528: - -tao_ntr_theta - theta interpolation update
530: Level: beginner
531: M*/
533: PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
534: {
535: TAO_NTR *tr;
538: PetscNew(&tr);
540: tao->ops->setup = TaoSetUp_NTR;
541: tao->ops->solve = TaoSolve_NTR;
542: tao->ops->setfromoptions = TaoSetFromOptions_NTR;
543: tao->ops->destroy = TaoDestroy_NTR;
545: /* Override default settings (unless already changed) */
546: if (!tao->max_it_changed) tao->max_it = 50;
547: if (!tao->trust0_changed) tao->trust0 = 100.0;
548: tao->data = (void *)tr;
550: /* Standard trust region update parameters */
551: tr->eta1 = 1.0e-4;
552: tr->eta2 = 0.25;
553: tr->eta3 = 0.50;
554: tr->eta4 = 0.90;
556: tr->alpha1 = 0.25;
557: tr->alpha2 = 0.50;
558: tr->alpha3 = 1.00;
559: tr->alpha4 = 2.00;
560: tr->alpha5 = 4.00;
562: /* Interpolation trust region update parameters */
563: tr->mu1 = 0.10;
564: tr->mu2 = 0.50;
566: tr->gamma1 = 0.25;
567: tr->gamma2 = 0.50;
568: tr->gamma3 = 2.00;
569: tr->gamma4 = 4.00;
571: tr->theta = 0.05;
573: /* Interpolation parameters for initialization */
574: tr->mu1_i = 0.35;
575: tr->mu2_i = 0.50;
577: tr->gamma1_i = 0.0625;
578: tr->gamma2_i = 0.50;
579: tr->gamma3_i = 2.00;
580: tr->gamma4_i = 5.00;
582: tr->theta_i = 0.25;
584: tr->min_radius = 1.0e-10;
585: tr->max_radius = 1.0e10;
586: tr->epsilon = 1.0e-6;
588: tr->init_type = NTR_INIT_INTERPOLATION;
589: tr->update_type = NTR_UPDATE_REDUCTION;
591: /* Set linear solver to default for trust region */
592: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
593: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
594: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
595: KSPAppendOptionsPrefix(tao->ksp, "tao_ntr_");
596: KSPSetType(tao->ksp, KSPSTCG);
597: return 0;
598: }