Actual source code: ntr.c
petsc-3.12.5 2020-03-29
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: 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: PetscInfo(tao,"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,KSPNASH,&is_nash);
70: PetscStrcmp(ksp_type,KSPSTCG,&is_stcg);
71: PetscStrcmp(ksp_type,KSPGLTR,&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 - (PetscReal)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 - (PetscReal)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: /* Call general purpose update function */
251: if (tao->ops->update) {
252: (*tao->ops->update)(tao, tao->niter, tao->user_update);
253: }
254: ++tao->niter;
255: tao->ksp_its=0;
256: /* Compute the Hessian */
257: if (needH) {
258: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
259: needH = 0;
260: }
262: if (tr->bfgs_pre) {
263: /* Update the limited memory preconditioner */
264: MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
265: ++bfgsUpdates;
266: }
268: while (tao->reason == TAO_CONTINUE_ITERATING) {
269: KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);
271: /* Solve the trust region subproblem */
272: KSPCGSetRadius(tao->ksp,tao->trust);
273: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
274: KSPGetIterationNumber(tao->ksp,&its);
275: tao->ksp_its+=its;
276: tao->ksp_tot_its+=its;
277: KSPCGGetNormD(tao->ksp, &norm_d);
279: if (0.0 == tao->trust) {
280: /* Radius was uninitialized; use the norm of the direction */
281: if (norm_d > 0.0) {
282: tao->trust = norm_d;
284: /* Modify the radius if it is too large or small */
285: tao->trust = PetscMax(tao->trust, tr->min_radius);
286: tao->trust = PetscMin(tao->trust, tr->max_radius);
287: }
288: else {
289: /* The direction was bad; set radius to default value and re-solve
290: the trust-region subproblem to get a direction */
291: tao->trust = tao->trust0;
293: /* Modify the radius if it is too large or small */
294: tao->trust = PetscMax(tao->trust, tr->min_radius);
295: tao->trust = PetscMin(tao->trust, tr->max_radius);
297: KSPCGSetRadius(tao->ksp,tao->trust);
298: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
299: KSPGetIterationNumber(tao->ksp,&its);
300: tao->ksp_its+=its;
301: tao->ksp_tot_its+=its;
302: KSPCGGetNormD(tao->ksp, &norm_d);
304: if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
305: }
306: }
307: VecScale(tao->stepdirection, -1.0);
308: KSPGetConvergedReason(tao->ksp, &ksp_reason);
309: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) {
310: /* Preconditioner is numerically indefinite; reset the
311: approximate if using BFGS preconditioning. */
312: MatLMVMReset(tr->M, PETSC_FALSE);
313: MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
314: bfgsUpdates = 1;
315: }
317: if (NTR_UPDATE_REDUCTION == tr->update_type) {
318: /* Get predicted reduction */
319: KSPCGGetObjFcn(tao->ksp,&prered);
320: if (prered >= 0.0) {
321: /* The predicted reduction has the wrong sign. This cannot
322: happen in infinite precision arithmetic. Step should
323: be rejected! */
324: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
325: }
326: else {
327: /* Compute trial step and function value */
328: VecCopy(tao->solution,tr->W);
329: VecAXPY(tr->W, 1.0, tao->stepdirection);
330: TaoComputeObjective(tao, tr->W, &ftrial);
332: if (PetscIsInfOrNanReal(ftrial)) {
333: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
334: } else {
335: /* Compute and actual reduction */
336: actred = f - ftrial;
337: prered = -prered;
338: if ((PetscAbsScalar(actred) <= tr->epsilon) &&
339: (PetscAbsScalar(prered) <= tr->epsilon)) {
340: kappa = 1.0;
341: }
342: else {
343: kappa = actred / prered;
344: }
346: /* Accept or reject the step and update radius */
347: if (kappa < tr->eta1) {
348: /* Reject the step */
349: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
350: }
351: else {
352: /* Accept the step */
353: if (kappa < tr->eta2) {
354: /* Marginal bad step */
355: tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
356: }
357: else if (kappa < tr->eta3) {
358: /* Reasonable step */
359: tao->trust = tr->alpha3 * tao->trust;
360: }
361: else if (kappa < tr->eta4) {
362: /* Good step */
363: tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
364: }
365: else {
366: /* Very good step */
367: tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
368: }
369: break;
370: }
371: }
372: }
373: }
374: else {
375: /* Get predicted reduction */
376: KSPCGGetObjFcn(tao->ksp,&prered);
377: if (prered >= 0.0) {
378: /* The predicted reduction has the wrong sign. This cannot
379: happen in infinite precision arithmetic. Step should
380: be rejected! */
381: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
382: }
383: else {
384: VecCopy(tao->solution, tr->W);
385: VecAXPY(tr->W, 1.0, tao->stepdirection);
386: TaoComputeObjective(tao, tr->W, &ftrial);
387: if (PetscIsInfOrNanReal(ftrial)) {
388: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
389: }
390: else {
391: VecDot(tao->gradient, tao->stepdirection, &beta);
392: actred = f - ftrial;
393: prered = -prered;
394: if ((PetscAbsScalar(actred) <= tr->epsilon) &&
395: (PetscAbsScalar(prered) <= tr->epsilon)) {
396: kappa = 1.0;
397: }
398: else {
399: kappa = actred / prered;
400: }
402: tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
403: tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
404: tau_min = PetscMin(tau_1, tau_2);
405: tau_max = PetscMax(tau_1, tau_2);
407: if (kappa >= 1.0 - tr->mu1) {
408: /* Great agreement; accept step and update radius */
409: if (tau_max < 1.0) {
410: tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
411: }
412: else if (tau_max > tr->gamma4) {
413: tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
414: }
415: else {
416: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
417: }
418: break;
419: }
420: else if (kappa >= 1.0 - tr->mu2) {
421: /* Good agreement */
423: if (tau_max < tr->gamma2) {
424: tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
425: }
426: else if (tau_max > tr->gamma3) {
427: tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
428: }
429: else if (tau_max < 1.0) {
430: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
431: }
432: else {
433: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
434: }
435: break;
436: }
437: else {
438: /* Not good agreement */
439: if (tau_min > 1.0) {
440: tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
441: }
442: else if (tau_max < tr->gamma1) {
443: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
444: }
445: else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
446: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
447: }
448: else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) &&
449: ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
450: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
451: }
452: else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) &&
453: ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
454: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
455: }
456: else {
457: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
458: }
459: }
460: }
461: }
462: }
464: /* The step computed was not good and the radius was decreased.
465: Monitor the radius to terminate. */
466: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
467: TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);
468: (*tao->ops->convergencetest)(tao,tao->cnvP);
469: }
471: /* The radius may have been increased; modify if it is too large */
472: tao->trust = PetscMin(tao->trust, tr->max_radius);
474: if (tao->reason == TAO_CONTINUE_ITERATING) {
475: VecCopy(tr->W, tao->solution);
476: f = ftrial;
477: TaoComputeGradient(tao, tao->solution, tao->gradient);
478: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
479: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
480: needH = 1;
481: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
482: TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);
483: (*tao->ops->convergencetest)(tao,tao->cnvP);
484: }
485: }
486: return(0);
487: }
489: /*------------------------------------------------------------*/
490: static PetscErrorCode TaoSetUp_NTR(Tao tao)
491: {
492: TAO_NTR *tr = (TAO_NTR *)tao->data;
496: if (!tao->gradient) {VecDuplicate(tao->solution, &tao->gradient);}
497: if (!tao->stepdirection) {VecDuplicate(tao->solution, &tao->stepdirection);}
498: if (!tr->W) {VecDuplicate(tao->solution, &tr->W);}
500: tr->bfgs_pre = 0;
501: tr->M = 0;
502: return(0);
503: }
505: /*------------------------------------------------------------*/
506: static PetscErrorCode TaoDestroy_NTR(Tao tao)
507: {
508: TAO_NTR *tr = (TAO_NTR *)tao->data;
512: if (tao->setupcalled) {
513: VecDestroy(&tr->W);
514: }
515: PetscFree(tao->data);
516: return(0);
517: }
519: /*------------------------------------------------------------*/
520: static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao)
521: {
522: TAO_NTR *tr = (TAO_NTR *)tao->data;
526: PetscOptionsHead(PetscOptionsObject,"Newton trust region method for unconstrained optimization");
527: PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type,NULL);
528: PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type,NULL);
529: PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL);
530: PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL);
531: PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL);
532: PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL);
533: PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL);
534: PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL);
535: PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL);
536: PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL);
537: PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL);
538: PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL);
539: PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL);
540: PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL);
541: PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL);
542: PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL);
543: PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL);
544: PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL);
545: PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL);
546: PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL);
547: PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL);
548: PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL);
549: PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL);
550: PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL);
551: PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL);
552: PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL);
553: PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL);
554: PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL);
555: PetscOptionsTail();
556: KSPSetFromOptions(tao->ksp);
557: return(0);
558: }
560: /*------------------------------------------------------------*/
561: /*MC
562: TAONTR - Newton's method with trust region for unconstrained minimization.
563: At each iteration, the Newton trust region method solves the system.
564: NTR expects a KSP solver with a trust region radius.
565: min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k
567: Options Database Keys:
568: + -tao_ntr_init_type - "constant","direction","interpolation"
569: . -tao_ntr_update_type - "reduction","interpolation"
570: . -tao_ntr_min_radius - lower bound on trust region radius
571: . -tao_ntr_max_radius - upper bound on trust region radius
572: . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
573: . -tao_ntr_mu1_i - mu1 interpolation init factor
574: . -tao_ntr_mu2_i - mu2 interpolation init factor
575: . -tao_ntr_gamma1_i - gamma1 interpolation init factor
576: . -tao_ntr_gamma2_i - gamma2 interpolation init factor
577: . -tao_ntr_gamma3_i - gamma3 interpolation init factor
578: . -tao_ntr_gamma4_i - gamma4 interpolation init factor
579: . -tao_ntr_theta_i - theta1 interpolation init factor
580: . -tao_ntr_eta1 - eta1 reduction update factor
581: . -tao_ntr_eta2 - eta2 reduction update factor
582: . -tao_ntr_eta3 - eta3 reduction update factor
583: . -tao_ntr_eta4 - eta4 reduction update factor
584: . -tao_ntr_alpha1 - alpha1 reduction update factor
585: . -tao_ntr_alpha2 - alpha2 reduction update factor
586: . -tao_ntr_alpha3 - alpha3 reduction update factor
587: . -tao_ntr_alpha4 - alpha4 reduction update factor
588: . -tao_ntr_alpha4 - alpha4 reduction update factor
589: . -tao_ntr_mu1 - mu1 interpolation update
590: . -tao_ntr_mu2 - mu2 interpolation update
591: . -tao_ntr_gamma1 - gamma1 interpolcation update
592: . -tao_ntr_gamma2 - gamma2 interpolcation update
593: . -tao_ntr_gamma3 - gamma3 interpolcation update
594: . -tao_ntr_gamma4 - gamma4 interpolation update
595: - -tao_ntr_theta - theta interpolation update
597: Level: beginner
598: M*/
600: PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
601: {
602: TAO_NTR *tr;
607: PetscNewLog(tao,&tr);
609: tao->ops->setup = TaoSetUp_NTR;
610: tao->ops->solve = TaoSolve_NTR;
611: tao->ops->setfromoptions = TaoSetFromOptions_NTR;
612: tao->ops->destroy = TaoDestroy_NTR;
614: /* Override default settings (unless already changed) */
615: if (!tao->max_it_changed) tao->max_it = 50;
616: if (!tao->trust0_changed) tao->trust0 = 100.0;
617: tao->data = (void*)tr;
619: /* Standard trust region update parameters */
620: tr->eta1 = 1.0e-4;
621: tr->eta2 = 0.25;
622: tr->eta3 = 0.50;
623: tr->eta4 = 0.90;
625: tr->alpha1 = 0.25;
626: tr->alpha2 = 0.50;
627: tr->alpha3 = 1.00;
628: tr->alpha4 = 2.00;
629: tr->alpha5 = 4.00;
631: /* Interpolation trust region update parameters */
632: tr->mu1 = 0.10;
633: tr->mu2 = 0.50;
635: tr->gamma1 = 0.25;
636: tr->gamma2 = 0.50;
637: tr->gamma3 = 2.00;
638: tr->gamma4 = 4.00;
640: tr->theta = 0.05;
642: /* Interpolation parameters for initialization */
643: tr->mu1_i = 0.35;
644: tr->mu2_i = 0.50;
646: tr->gamma1_i = 0.0625;
647: tr->gamma2_i = 0.50;
648: tr->gamma3_i = 2.00;
649: tr->gamma4_i = 5.00;
651: tr->theta_i = 0.25;
653: tr->min_radius = 1.0e-10;
654: tr->max_radius = 1.0e10;
655: tr->epsilon = 1.0e-6;
657: tr->init_type = NTR_INIT_INTERPOLATION;
658: tr->update_type = NTR_UPDATE_REDUCTION;
660: /* Set linear solver to default for trust region */
661: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
662: PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1);
663: KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
664: KSPAppendOptionsPrefix(tao->ksp,"tao_ntr_");
665: KSPSetType(tao->ksp,KSPSTCG);
666: return(0);
667: }