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