Actual source code: ntr.c
petsc-3.10.5 2019-03-28
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 KSPCGNASH, KSPCGSTCG,
40: or KSPCGGLTR. Thus, we set KSPCGNASH, KSPCGSTCG, or KSPCGGLTR in this
41: routine regardless of what the user may have previously specified.
42: */
43: static PetscErrorCode TaoSolve_NTR(Tao tao)
44: {
45: TAO_NTR *tr = (TAO_NTR *)tao->data;
46: KSPType ksp_type;
47: PetscBool is_nash,is_stcg,is_gltr,is_bfgs,is_jacobi,is_symmetric,sym_set;
48: KSPConvergedReason ksp_reason;
49: PC pc;
50: PetscReal fmin, ftrial, prered, actred, kappa, sigma, beta;
51: PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
52: PetscReal f, gnorm;
54: PetscReal norm_d;
55: PetscErrorCode ierr;
56: PetscInt bfgsUpdates = 0;
57: PetscInt needH;
59: PetscInt i_max = 5;
60: PetscInt j_max = 1;
61: PetscInt i, j, N, n, its;
64: if (tao->XL || tao->XU || tao->ops->computebounds) {
65: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");
66: }
68: KSPGetType(tao->ksp,&ksp_type);
69: PetscStrcmp(ksp_type,KSPCGNASH,&is_nash);
70: PetscStrcmp(ksp_type,KSPCGSTCG,&is_stcg);
71: PetscStrcmp(ksp_type,KSPCGGLTR,&is_gltr);
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: }
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);
81: /* Allocate the vectors needed for the BFGS approximation */
82: KSPGetPC(tao->ksp, &pc);
83: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);
84: PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);
85: if (is_bfgs) {
86: tr->bfgs_pre = pc;
87: PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M);
88: VecGetLocalSize(tao->solution, &n);
89: VecGetSize(tao->solution, &N);
90: MatSetSizes(tr->M, n, n, N, N);
91: MatLMVMAllocate(tr->M, tao->solution, tao->gradient);
92: MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric);
93: if (!sym_set || !is_symmetric) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
94: } else if (is_jacobi) {
95: PCJacobiSetUseAbs(pc,PETSC_TRUE);
96: }
98: /* Check convergence criteria */
99: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
100: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
101: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Inf or NaN");
102: needH = 1;
104: tao->reason = TAO_CONTINUE_ITERATING;
105: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
106: TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);
107: (*tao->ops->convergencetest)(tao,tao->cnvP);
108: if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
110: /* Initialize trust-region radius */
111: switch(tr->init_type) {
112: case NTR_INIT_CONSTANT:
113: /* Use the initial radius specified */
114: break;
116: case NTR_INIT_INTERPOLATION:
117: /* Use the initial radius specified */
118: max_radius = 0.0;
120: for (j = 0; j < j_max; ++j) {
121: fmin = f;
122: sigma = 0.0;
124: if (needH) {
125: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
126: needH = 0;
127: }
129: for (i = 0; i < i_max; ++i) {
131: VecCopy(tao->solution, tr->W);
132: VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);
133: TaoComputeObjective(tao, tr->W, &ftrial);
135: if (PetscIsInfOrNanReal(ftrial)) {
136: tau = tr->gamma1_i;
137: }
138: else {
139: if (ftrial < fmin) {
140: fmin = ftrial;
141: sigma = -tao->trust / gnorm;
142: }
144: MatMult(tao->hessian, tao->gradient, tao->stepdirection);
145: VecDot(tao->gradient, tao->stepdirection, &prered);
147: prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
148: actred = f - ftrial;
149: if ((PetscAbsScalar(actred) <= tr->epsilon) &&
150: (PetscAbsScalar(prered) <= tr->epsilon)) {
151: kappa = 1.0;
152: }
153: else {
154: kappa = actred / prered;
155: }
157: tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
158: tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
159: tau_min = PetscMin(tau_1, tau_2);
160: tau_max = PetscMax(tau_1, tau_2);
162: if (PetscAbsScalar(kappa - 1.0) <= tr->mu1_i) {
163: /* Great agreement */
164: max_radius = PetscMax(max_radius, tao->trust);
166: if (tau_max < 1.0) {
167: tau = tr->gamma3_i;
168: }
169: else if (tau_max > tr->gamma4_i) {
170: tau = tr->gamma4_i;
171: }
172: else {
173: tau = tau_max;
174: }
175: }
176: else if (PetscAbsScalar(kappa - 1.0) <= tr->mu2_i) {
177: /* Good agreement */
178: max_radius = PetscMax(max_radius, tao->trust);
180: if (tau_max < tr->gamma2_i) {
181: tau = tr->gamma2_i;
182: }
183: else if (tau_max > tr->gamma3_i) {
184: tau = tr->gamma3_i;
185: }
186: else {
187: tau = tau_max;
188: }
189: }
190: else {
191: /* Not good agreement */
192: if (tau_min > 1.0) {
193: tau = tr->gamma2_i;
194: }
195: else if (tau_max < tr->gamma1_i) {
196: tau = tr->gamma1_i;
197: }
198: else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
199: tau = tr->gamma1_i;
200: }
201: else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) &&
202: ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
203: tau = tau_1;
204: }
205: else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) &&
206: ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
207: tau = tau_2;
208: }
209: else {
210: tau = tau_max;
211: }
212: }
213: }
214: tao->trust = tau * tao->trust;
215: }
217: if (fmin < f) {
218: f = fmin;
219: VecAXPY(tao->solution, sigma, tao->gradient);
220: TaoComputeGradient(tao,tao->solution, tao->gradient);
222: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
224: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
225: needH = 1;
227: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
228: TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);
229: (*tao->ops->convergencetest)(tao,tao->cnvP);
230: if (tao->reason != TAO_CONTINUE_ITERATING) {
231: return(0);
232: }
233: }
234: }
235: tao->trust = PetscMax(tao->trust, max_radius);
237: /* Modify the radius if it is too large or small */
238: tao->trust = PetscMax(tao->trust, tr->min_radius);
239: tao->trust = PetscMin(tao->trust, tr->max_radius);
240: break;
242: default:
243: /* Norm of the first direction will initialize radius */
244: tao->trust = 0.0;
245: break;
246: }
248: /* Have not converged; continue with Newton method */
249: while (tao->reason == TAO_CONTINUE_ITERATING) {
250: ++tao->niter;
251: tao->ksp_its=0;
252: /* Compute the Hessian */
253: if (needH) {
254: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
255: needH = 0;
256: }
258: if (tr->bfgs_pre) {
259: /* Update the limited memory preconditioner */
260: MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
261: ++bfgsUpdates;
262: }
264: while (tao->reason == TAO_CONTINUE_ITERATING) {
265: KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);
267: /* Solve the trust region subproblem */
268: KSPCGSetRadius(tao->ksp,tao->trust);
269: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
270: KSPGetIterationNumber(tao->ksp,&its);
271: tao->ksp_its+=its;
272: tao->ksp_tot_its+=its;
273: KSPCGGetNormD(tao->ksp, &norm_d);
275: if (0.0 == tao->trust) {
276: /* Radius was uninitialized; use the norm of the direction */
277: if (norm_d > 0.0) {
278: tao->trust = norm_d;
280: /* Modify the radius if it is too large or small */
281: tao->trust = PetscMax(tao->trust, tr->min_radius);
282: tao->trust = PetscMin(tao->trust, tr->max_radius);
283: }
284: else {
285: /* The direction was bad; set radius to default value and re-solve
286: the trust-region subproblem to get a direction */
287: tao->trust = tao->trust0;
289: /* Modify the radius if it is too large or small */
290: tao->trust = PetscMax(tao->trust, tr->min_radius);
291: tao->trust = PetscMin(tao->trust, tr->max_radius);
293: KSPCGSetRadius(tao->ksp,tao->trust);
294: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
295: KSPGetIterationNumber(tao->ksp,&its);
296: tao->ksp_its+=its;
297: tao->ksp_tot_its+=its;
298: KSPCGGetNormD(tao->ksp, &norm_d);
300: if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
301: }
302: }
303: VecScale(tao->stepdirection, -1.0);
304: KSPGetConvergedReason(tao->ksp, &ksp_reason);
305: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) {
306: /* Preconditioner is numerically indefinite; reset the
307: approximate if using BFGS preconditioning. */
308: MatLMVMReset(tr->M, PETSC_FALSE);
309: MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
310: bfgsUpdates = 1;
311: }
313: if (NTR_UPDATE_REDUCTION == tr->update_type) {
314: /* Get predicted reduction */
315: KSPCGGetObjFcn(tao->ksp,&prered);
316: if (prered >= 0.0) {
317: /* The predicted reduction has the wrong sign. This cannot
318: happen in infinite precision arithmetic. Step should
319: be rejected! */
320: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
321: }
322: else {
323: /* Compute trial step and function value */
324: VecCopy(tao->solution,tr->W);
325: VecAXPY(tr->W, 1.0, tao->stepdirection);
326: TaoComputeObjective(tao, tr->W, &ftrial);
328: if (PetscIsInfOrNanReal(ftrial)) {
329: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
330: } else {
331: /* Compute and actual reduction */
332: actred = f - ftrial;
333: prered = -prered;
334: if ((PetscAbsScalar(actred) <= tr->epsilon) &&
335: (PetscAbsScalar(prered) <= tr->epsilon)) {
336: kappa = 1.0;
337: }
338: else {
339: kappa = actred / prered;
340: }
342: /* Accept or reject the step and update radius */
343: if (kappa < tr->eta1) {
344: /* Reject the step */
345: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
346: }
347: else {
348: /* Accept the step */
349: if (kappa < tr->eta2) {
350: /* Marginal bad step */
351: tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
352: }
353: else if (kappa < tr->eta3) {
354: /* Reasonable step */
355: tao->trust = tr->alpha3 * tao->trust;
356: }
357: else if (kappa < tr->eta4) {
358: /* Good step */
359: tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
360: }
361: else {
362: /* Very good step */
363: tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
364: }
365: break;
366: }
367: }
368: }
369: }
370: else {
371: /* Get predicted reduction */
372: KSPCGGetObjFcn(tao->ksp,&prered);
373: if (prered >= 0.0) {
374: /* The predicted reduction has the wrong sign. This cannot
375: happen in infinite precision arithmetic. Step should
376: be rejected! */
377: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
378: }
379: else {
380: VecCopy(tao->solution, tr->W);
381: VecAXPY(tr->W, 1.0, tao->stepdirection);
382: TaoComputeObjective(tao, tr->W, &ftrial);
383: if (PetscIsInfOrNanReal(ftrial)) {
384: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
385: }
386: else {
387: VecDot(tao->gradient, tao->stepdirection, &beta);
388: actred = f - ftrial;
389: prered = -prered;
390: if ((PetscAbsScalar(actred) <= tr->epsilon) &&
391: (PetscAbsScalar(prered) <= tr->epsilon)) {
392: kappa = 1.0;
393: }
394: else {
395: kappa = actred / prered;
396: }
398: tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
399: tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
400: tau_min = PetscMin(tau_1, tau_2);
401: tau_max = PetscMax(tau_1, tau_2);
403: if (kappa >= 1.0 - tr->mu1) {
404: /* Great agreement; accept step and update radius */
405: if (tau_max < 1.0) {
406: tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
407: }
408: else if (tau_max > tr->gamma4) {
409: tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
410: }
411: else {
412: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
413: }
414: break;
415: }
416: else if (kappa >= 1.0 - tr->mu2) {
417: /* Good agreement */
419: if (tau_max < tr->gamma2) {
420: tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
421: }
422: else if (tau_max > tr->gamma3) {
423: tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
424: }
425: else if (tau_max < 1.0) {
426: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
427: }
428: else {
429: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
430: }
431: break;
432: }
433: else {
434: /* Not good agreement */
435: if (tau_min > 1.0) {
436: tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
437: }
438: else if (tau_max < tr->gamma1) {
439: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
440: }
441: else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
442: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
443: }
444: else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) &&
445: ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
446: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
447: }
448: else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) &&
449: ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
450: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
451: }
452: else {
453: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
454: }
455: }
456: }
457: }
458: }
460: /* The step computed was not good and the radius was decreased.
461: Monitor the radius to terminate. */
462: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
463: TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);
464: (*tao->ops->convergencetest)(tao,tao->cnvP);
465: }
467: /* The radius may have been increased; modify if it is too large */
468: tao->trust = PetscMin(tao->trust, tr->max_radius);
470: if (tao->reason == TAO_CONTINUE_ITERATING) {
471: VecCopy(tr->W, tao->solution);
472: f = ftrial;
473: TaoComputeGradient(tao, tao->solution, tao->gradient);
474: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
475: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
476: needH = 1;
477: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
478: TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);
479: (*tao->ops->convergencetest)(tao,tao->cnvP);
480: }
481: }
482: return(0);
483: }
485: /*------------------------------------------------------------*/
486: static PetscErrorCode TaoSetUp_NTR(Tao tao)
487: {
488: TAO_NTR *tr = (TAO_NTR *)tao->data;
492: if (!tao->gradient) {VecDuplicate(tao->solution, &tao->gradient);}
493: if (!tao->stepdirection) {VecDuplicate(tao->solution, &tao->stepdirection);}
494: if (!tr->W) {VecDuplicate(tao->solution, &tr->W);}
496: tr->bfgs_pre = 0;
497: tr->M = 0;
498: return(0);
499: }
501: /*------------------------------------------------------------*/
502: static PetscErrorCode TaoDestroy_NTR(Tao tao)
503: {
504: TAO_NTR *tr = (TAO_NTR *)tao->data;
508: if (tao->setupcalled) {
509: VecDestroy(&tr->W);
510: }
511: PetscFree(tao->data);
512: return(0);
513: }
515: /*------------------------------------------------------------*/
516: static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao)
517: {
518: TAO_NTR *tr = (TAO_NTR *)tao->data;
522: PetscOptionsHead(PetscOptionsObject,"Newton trust region method for unconstrained optimization");
523: PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type,NULL);
524: PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type,NULL);
525: PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL);
526: PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL);
527: PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL);
528: PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL);
529: PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL);
530: PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL);
531: PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL);
532: PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL);
533: PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL);
534: PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL);
535: PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL);
536: PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL);
537: PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL);
538: PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL);
539: PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL);
540: PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL);
541: PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL);
542: PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL);
543: PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL);
544: PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL);
545: PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL);
546: PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL);
547: PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL);
548: PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL);
549: PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL);
550: PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL);
551: PetscOptionsTail();
552: KSPSetFromOptions(tao->ksp);
553: return(0);
554: }
556: /*------------------------------------------------------------*/
557: /*MC
558: TAONTR - Newton's method with trust region for unconstrained minimization.
559: At each iteration, the Newton trust region method solves the system.
560: NTR expects a KSP solver with a trust region radius.
561: min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k
563: Options Database Keys:
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
593: Level: beginner
594: M*/
596: PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
597: {
598: TAO_NTR *tr;
603: PetscNewLog(tao,&tr);
605: tao->ops->setup = TaoSetUp_NTR;
606: tao->ops->solve = TaoSolve_NTR;
607: tao->ops->setfromoptions = TaoSetFromOptions_NTR;
608: tao->ops->destroy = TaoDestroy_NTR;
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;
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;
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;
627: /* Interpolation trust region update parameters */
628: tr->mu1 = 0.10;
629: tr->mu2 = 0.50;
631: tr->gamma1 = 0.25;
632: tr->gamma2 = 0.50;
633: tr->gamma3 = 2.00;
634: tr->gamma4 = 4.00;
636: tr->theta = 0.05;
638: /* Interpolation parameters for initialization */
639: tr->mu1_i = 0.35;
640: tr->mu2_i = 0.50;
642: tr->gamma1_i = 0.0625;
643: tr->gamma2_i = 0.50;
644: tr->gamma3_i = 2.00;
645: tr->gamma4_i = 5.00;
647: tr->theta_i = 0.25;
649: tr->min_radius = 1.0e-10;
650: tr->max_radius = 1.0e10;
651: tr->epsilon = 1.0e-6;
653: tr->init_type = NTR_INIT_INTERPOLATION;
654: tr->update_type = NTR_UPDATE_REDUCTION;
656: /* Set linear solver to default for trust region */
657: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
658: PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1);
659: KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
660: KSPAppendOptionsPrefix(tao->ksp,"tao_ntr_");
661: KSPSetType(tao->ksp,KSPCGSTCG);
662: return(0);
663: }