Actual source code: ntr.c
petsc-3.6.4 2016-04-12
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 bfgsUpdates = 0;
87: PetscInt needH;
89: PetscInt i_max = 5;
90: PetscInt j_max = 1;
91: PetscInt i, j, N, n, its;
94: if (tao->XL || tao->XU || tao->ops->computebounds) {
95: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");
96: }
98: tao->trust = tao->trust0;
100: /* Modify the radius if it is too large or small */
101: tao->trust = PetscMax(tao->trust, tr->min_radius);
102: tao->trust = PetscMin(tao->trust, tr->max_radius);
105: if (NTR_PC_BFGS == tr->pc_type && !tr->M) {
106: VecGetLocalSize(tao->solution,&n);
107: VecGetSize(tao->solution,&N);
108: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tr->M);
109: MatLMVMAllocateVectors(tr->M,tao->solution);
110: }
112: /* Check convergence criteria */
113: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
114: VecNorm(tao->gradient,NORM_2,&gnorm);
115: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
116: needH = 1;
118: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0, &reason);
119: if (reason != TAO_CONTINUE_ITERATING) return(0);
121: /* Create vectors for the limited memory preconditioner */
122: if ((NTR_PC_BFGS == tr->pc_type) &&
123: (BFGS_SCALE_BFGS != tr->bfgs_scale_type)) {
124: if (!tr->Diag) {
125: VecDuplicate(tao->solution, &tr->Diag);
126: }
127: }
129: switch(tr->ksp_type) {
130: case NTR_KSP_NASH:
131: KSPSetType(tao->ksp, KSPNASH);
132: KSPSetFromOptions(tao->ksp);
133: break;
135: case NTR_KSP_STCG:
136: KSPSetType(tao->ksp, KSPSTCG);
137: KSPSetFromOptions(tao->ksp);
138: break;
140: default:
141: KSPSetType(tao->ksp, KSPGLTR);
142: KSPSetFromOptions(tao->ksp);
143: break;
144: }
146: /* Modify the preconditioner to use the bfgs approximation */
147: KSPGetPC(tao->ksp, &pc);
148: switch(tr->pc_type) {
149: case NTR_PC_NONE:
150: PCSetType(pc, PCNONE);
151: PCSetFromOptions(pc);
152: break;
154: case NTR_PC_AHESS:
155: PCSetType(pc, PCJACOBI);
156: PCSetFromOptions(pc);
157: PCJacobiSetUseAbs(pc,PETSC_TRUE);
158: break;
160: case NTR_PC_BFGS:
161: PCSetType(pc, PCSHELL);
162: PCSetFromOptions(pc);
163: PCShellSetName(pc, "bfgs");
164: PCShellSetContext(pc, tr->M);
165: PCShellSetApply(pc, MatLMVMSolveShell);
166: break;
168: default:
169: /* Use the pc method set by pc_type */
170: break;
171: }
173: /* Initialize trust-region radius */
174: switch(tr->init_type) {
175: case NTR_INIT_CONSTANT:
176: /* Use the initial radius specified */
177: break;
179: case NTR_INIT_INTERPOLATION:
180: /* Use the initial radius specified */
181: max_radius = 0.0;
183: for (j = 0; j < j_max; ++j) {
184: fmin = f;
185: sigma = 0.0;
187: if (needH) {
188: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
189: needH = 0;
190: }
192: for (i = 0; i < i_max; ++i) {
194: VecCopy(tao->solution, tr->W);
195: VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);
196: TaoComputeObjective(tao, tr->W, &ftrial);
198: if (PetscIsInfOrNanReal(ftrial)) {
199: tau = tr->gamma1_i;
200: }
201: else {
202: if (ftrial < fmin) {
203: fmin = ftrial;
204: sigma = -tao->trust / gnorm;
205: }
207: MatMult(tao->hessian, tao->gradient, tao->stepdirection);
208: VecDot(tao->gradient, tao->stepdirection, &prered);
210: prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
211: actred = f - ftrial;
212: if ((PetscAbsScalar(actred) <= tr->epsilon) &&
213: (PetscAbsScalar(prered) <= tr->epsilon)) {
214: kappa = 1.0;
215: }
216: else {
217: kappa = actred / prered;
218: }
220: tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
221: tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
222: tau_min = PetscMin(tau_1, tau_2);
223: tau_max = PetscMax(tau_1, tau_2);
225: if (PetscAbsScalar(kappa - 1.0) <= tr->mu1_i) {
226: /* Great agreement */
227: max_radius = PetscMax(max_radius, tao->trust);
229: if (tau_max < 1.0) {
230: tau = tr->gamma3_i;
231: }
232: else if (tau_max > tr->gamma4_i) {
233: tau = tr->gamma4_i;
234: }
235: else {
236: tau = tau_max;
237: }
238: }
239: else if (PetscAbsScalar(kappa - 1.0) <= tr->mu2_i) {
240: /* Good agreement */
241: max_radius = PetscMax(max_radius, tao->trust);
243: if (tau_max < tr->gamma2_i) {
244: tau = tr->gamma2_i;
245: }
246: else if (tau_max > tr->gamma3_i) {
247: tau = tr->gamma3_i;
248: }
249: else {
250: tau = tau_max;
251: }
252: }
253: else {
254: /* Not good agreement */
255: if (tau_min > 1.0) {
256: tau = tr->gamma2_i;
257: }
258: else if (tau_max < tr->gamma1_i) {
259: tau = tr->gamma1_i;
260: }
261: else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
262: tau = tr->gamma1_i;
263: }
264: else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) &&
265: ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
266: tau = tau_1;
267: }
268: else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) &&
269: ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
270: tau = tau_2;
271: }
272: else {
273: tau = tau_max;
274: }
275: }
276: }
277: tao->trust = tau * tao->trust;
278: }
280: if (fmin < f) {
281: f = fmin;
282: VecAXPY(tao->solution, sigma, tao->gradient);
283: TaoComputeGradient(tao,tao->solution, tao->gradient);
285: VecNorm(tao->gradient, NORM_2, &gnorm);
287: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
288: needH = 1;
290: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0, &reason);
291: if (reason != TAO_CONTINUE_ITERATING) {
292: return(0);
293: }
294: }
295: }
296: tao->trust = PetscMax(tao->trust, max_radius);
298: /* Modify the radius if it is too large or small */
299: tao->trust = PetscMax(tao->trust, tr->min_radius);
300: tao->trust = PetscMin(tao->trust, tr->max_radius);
301: break;
303: default:
304: /* Norm of the first direction will initialize radius */
305: tao->trust = 0.0;
306: break;
307: }
309: /* Set initial scaling for the BFGS preconditioner
310: This step is done after computing the initial trust-region radius
311: since the function value may have decreased */
312: if (NTR_PC_BFGS == tr->pc_type) {
313: if (f != 0.0) {
314: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
315: }
316: else {
317: delta = 2.0 / (gnorm*gnorm);
318: }
319: MatLMVMSetDelta(tr->M,delta);
320: }
322: /* Have not converged; continue with Newton method */
323: while (reason == TAO_CONTINUE_ITERATING) {
324: ++tao->niter;
325: tao->ksp_its=0;
326: /* Compute the Hessian */
327: if (needH) {
328: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
329: needH = 0;
330: }
332: if (NTR_PC_BFGS == tr->pc_type) {
333: if (BFGS_SCALE_AHESS == tr->bfgs_scale_type) {
334: /* Obtain diagonal for the bfgs preconditioner */
335: MatGetDiagonal(tao->hessian, tr->Diag);
336: VecAbs(tr->Diag);
337: VecReciprocal(tr->Diag);
338: MatLMVMSetScale(tr->M,tr->Diag);
339: }
341: /* Update the limited memory preconditioner */
342: MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
343: ++bfgsUpdates;
344: }
346: while (reason == TAO_CONTINUE_ITERATING) {
347: KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);
349: /* Solve the trust region subproblem */
350: if (NTR_KSP_NASH == tr->ksp_type) {
351: KSPNASHSetRadius(tao->ksp,tao->trust);
352: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
353: KSPGetIterationNumber(tao->ksp,&its);
354: tao->ksp_its+=its;
355: tao->ksp_tot_its+=its;
356: KSPNASHGetNormD(tao->ksp, &norm_d);
357: } else if (NTR_KSP_STCG == tr->ksp_type) {
358: KSPSTCGSetRadius(tao->ksp,tao->trust);
359: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
360: KSPGetIterationNumber(tao->ksp,&its);
361: tao->ksp_its+=its;
362: tao->ksp_tot_its+=its;
363: KSPSTCGGetNormD(tao->ksp, &norm_d);
364: } else { /* NTR_KSP_GLTR */
365: KSPGLTRSetRadius(tao->ksp,tao->trust);
366: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
367: KSPGetIterationNumber(tao->ksp,&its);
368: tao->ksp_its+=its;
369: tao->ksp_tot_its+=its;
370: KSPGLTRGetNormD(tao->ksp, &norm_d);
371: }
373: if (0.0 == tao->trust) {
374: /* Radius was uninitialized; use the norm of the direction */
375: if (norm_d > 0.0) {
376: tao->trust = norm_d;
378: /* Modify the radius if it is too large or small */
379: tao->trust = PetscMax(tao->trust, tr->min_radius);
380: tao->trust = PetscMin(tao->trust, tr->max_radius);
381: }
382: else {
383: /* The direction was bad; set radius to default value and re-solve
384: the trust-region subproblem to get a direction */
385: tao->trust = tao->trust0;
387: /* Modify the radius if it is too large or small */
388: tao->trust = PetscMax(tao->trust, tr->min_radius);
389: tao->trust = PetscMin(tao->trust, tr->max_radius);
391: if (NTR_KSP_NASH == tr->ksp_type) {
392: KSPNASHSetRadius(tao->ksp,tao->trust);
393: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
394: KSPGetIterationNumber(tao->ksp,&its);
395: tao->ksp_its+=its;
396: tao->ksp_tot_its+=its;
397: KSPNASHGetNormD(tao->ksp, &norm_d);
398: } else if (NTR_KSP_STCG == tr->ksp_type) {
399: KSPSTCGSetRadius(tao->ksp,tao->trust);
400: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
401: KSPGetIterationNumber(tao->ksp,&its);
402: tao->ksp_its+=its;
403: tao->ksp_tot_its+=its;
404: KSPSTCGGetNormD(tao->ksp, &norm_d);
405: } else { /* NTR_KSP_GLTR */
406: KSPGLTRSetRadius(tao->ksp,tao->trust);
407: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
408: KSPGetIterationNumber(tao->ksp,&its);
409: tao->ksp_its+=its;
410: tao->ksp_tot_its+=its;
411: KSPGLTRGetNormD(tao->ksp, &norm_d);
412: }
414: if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
415: }
416: }
417: VecScale(tao->stepdirection, -1.0);
418: KSPGetConvergedReason(tao->ksp, &ksp_reason);
419: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
420: (NTR_PC_BFGS == tr->pc_type) && (bfgsUpdates > 1)) {
421: /* Preconditioner is numerically indefinite; reset the
422: approximate if using BFGS preconditioning. */
424: if (f != 0.0) {
425: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
426: }
427: else {
428: delta = 2.0 / (gnorm*gnorm);
429: }
430: MatLMVMSetDelta(tr->M, delta);
431: MatLMVMReset(tr->M);
432: MatLMVMUpdate(tr->M, tao->solution, tao->gradient);
433: bfgsUpdates = 1;
434: }
436: if (NTR_UPDATE_REDUCTION == tr->update_type) {
437: /* Get predicted reduction */
438: if (NTR_KSP_NASH == tr->ksp_type) {
439: KSPNASHGetObjFcn(tao->ksp,&prered);
440: } else if (NTR_KSP_STCG == tr->ksp_type) {
441: KSPSTCGGetObjFcn(tao->ksp,&prered);
442: } else { /* gltr */
443: KSPGLTRGetObjFcn(tao->ksp,&prered);
444: }
446: if (prered >= 0.0) {
447: /* The predicted reduction has the wrong sign. This cannot
448: happen in infinite precision arithmetic. Step should
449: be rejected! */
450: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
451: }
452: else {
453: /* Compute trial step and function value */
454: VecCopy(tao->solution,tr->W);
455: VecAXPY(tr->W, 1.0, tao->stepdirection);
456: TaoComputeObjective(tao, tr->W, &ftrial);
458: if (PetscIsInfOrNanReal(ftrial)) {
459: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
460: } else {
461: /* Compute and actual reduction */
462: actred = f - ftrial;
463: prered = -prered;
464: if ((PetscAbsScalar(actred) <= tr->epsilon) &&
465: (PetscAbsScalar(prered) <= tr->epsilon)) {
466: kappa = 1.0;
467: }
468: else {
469: kappa = actred / prered;
470: }
472: /* Accept or reject the step and update radius */
473: if (kappa < tr->eta1) {
474: /* Reject the step */
475: tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
476: }
477: else {
478: /* Accept the step */
479: if (kappa < tr->eta2) {
480: /* Marginal bad step */
481: tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
482: }
483: else if (kappa < tr->eta3) {
484: /* Reasonable step */
485: tao->trust = tr->alpha3 * tao->trust;
486: }
487: else if (kappa < tr->eta4) {
488: /* Good step */
489: tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
490: }
491: else {
492: /* Very good step */
493: tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
494: }
495: break;
496: }
497: }
498: }
499: }
500: else {
501: /* Get predicted reduction */
502: if (NTR_KSP_NASH == tr->ksp_type) {
503: KSPNASHGetObjFcn(tao->ksp,&prered);
504: } else if (NTR_KSP_STCG == tr->ksp_type) {
505: KSPSTCGGetObjFcn(tao->ksp,&prered);
506: } else { /* gltr */
507: KSPGLTRGetObjFcn(tao->ksp,&prered);
508: }
510: if (prered >= 0.0) {
511: /* The predicted reduction has the wrong sign. This cannot
512: happen in infinite precision arithmetic. Step should
513: be rejected! */
514: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
515: }
516: else {
517: VecCopy(tao->solution, tr->W);
518: VecAXPY(tr->W, 1.0, tao->stepdirection);
519: TaoComputeObjective(tao, tr->W, &ftrial);
520: if (PetscIsInfOrNanReal(ftrial)) {
521: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
522: }
523: else {
524: VecDot(tao->gradient, tao->stepdirection, &beta);
525: actred = f - ftrial;
526: prered = -prered;
527: if ((PetscAbsScalar(actred) <= tr->epsilon) &&
528: (PetscAbsScalar(prered) <= tr->epsilon)) {
529: kappa = 1.0;
530: }
531: else {
532: kappa = actred / prered;
533: }
535: tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
536: tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
537: tau_min = PetscMin(tau_1, tau_2);
538: tau_max = PetscMax(tau_1, tau_2);
540: if (kappa >= 1.0 - tr->mu1) {
541: /* Great agreement; accept step and update radius */
542: if (tau_max < 1.0) {
543: tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
544: }
545: else if (tau_max > tr->gamma4) {
546: tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
547: }
548: else {
549: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
550: }
551: break;
552: }
553: else if (kappa >= 1.0 - tr->mu2) {
554: /* Good agreement */
556: if (tau_max < tr->gamma2) {
557: tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
558: }
559: else if (tau_max > tr->gamma3) {
560: tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
561: }
562: else if (tau_max < 1.0) {
563: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
564: }
565: else {
566: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
567: }
568: break;
569: }
570: else {
571: /* Not good agreement */
572: if (tau_min > 1.0) {
573: tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
574: }
575: else if (tau_max < tr->gamma1) {
576: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
577: }
578: else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
579: tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
580: }
581: else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) &&
582: ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
583: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
584: }
585: else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) &&
586: ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
587: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
588: }
589: else {
590: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
591: }
592: }
593: }
594: }
595: }
597: /* The step computed was not good and the radius was decreased.
598: Monitor the radius to terminate. */
599: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust, &reason);
600: }
602: /* The radius may have been increased; modify if it is too large */
603: tao->trust = PetscMin(tao->trust, tr->max_radius);
605: if (reason == TAO_CONTINUE_ITERATING) {
606: VecCopy(tr->W, tao->solution);
607: f = ftrial;
608: TaoComputeGradient(tao, tao->solution, tao->gradient);
609: VecNorm(tao->gradient, NORM_2, &gnorm);
610: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
611: needH = 1;
612: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust, &reason);
613: }
614: }
615: return(0);
616: }
618: /*------------------------------------------------------------*/
621: static PetscErrorCode TaoSetUp_NTR(Tao tao)
622: {
623: TAO_NTR *tr = (TAO_NTR *)tao->data;
628: if (!tao->gradient) {VecDuplicate(tao->solution, &tao->gradient);}
629: if (!tao->stepdirection) {VecDuplicate(tao->solution, &tao->stepdirection);}
630: if (!tr->W) {VecDuplicate(tao->solution, &tr->W);}
632: tr->Diag = 0;
633: tr->M = 0;
636: return(0);
637: }
639: /*------------------------------------------------------------*/
642: static PetscErrorCode TaoDestroy_NTR(Tao tao)
643: {
644: TAO_NTR *tr = (TAO_NTR *)tao->data;
648: if (tao->setupcalled) {
649: VecDestroy(&tr->W);
650: }
651: MatDestroy(&tr->M);
652: VecDestroy(&tr->Diag);
653: PetscFree(tao->data);
654: return(0);
655: }
657: /*------------------------------------------------------------*/
660: static PetscErrorCode TaoSetFromOptions_NTR(PetscOptions *PetscOptionsObject,Tao tao)
661: {
662: TAO_NTR *tr = (TAO_NTR *)tao->data;
666: PetscOptionsHead(PetscOptionsObject,"Newton trust region method for unconstrained optimization");
667: PetscOptionsEList("-tao_ntr_ksp_type", "ksp type", "", NTR_KSP, NTR_KSP_TYPES, NTR_KSP[tr->ksp_type], &tr->ksp_type,NULL);
668: PetscOptionsEList("-tao_ntr_pc_type", "pc type", "", NTR_PC, NTR_PC_TYPES, NTR_PC[tr->pc_type], &tr->pc_type,NULL);
669: PetscOptionsEList("-tao_ntr_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tr->bfgs_scale_type], &tr->bfgs_scale_type,NULL);
670: PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type,NULL);
671: PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type,NULL);
672: PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL);
673: PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL);
674: PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL);
675: PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL);
676: PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL);
677: PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL);
678: PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL);
679: PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL);
680: PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL);
681: PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL);
682: PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL);
683: PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL);
684: PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL);
685: PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL);
686: PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL);
687: PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL);
688: PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL);
689: PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL);
690: PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL);
691: PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL);
692: PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL);
693: PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL);
694: PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL);
695: PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL);
696: PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL);
697: PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL);
698: PetscOptionsTail();
699: KSPSetFromOptions(tao->ksp);
700: return(0);
701: }
703: /*------------------------------------------------------------*/
706: static PetscErrorCode TaoView_NTR(Tao tao, PetscViewer viewer)
707: {
708: TAO_NTR *tr = (TAO_NTR *)tao->data;
710: PetscInt nrejects;
711: PetscBool isascii;
714: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
715: if (isascii) {
716: PetscViewerASCIIPushTab(viewer);
717: if (NTR_PC_BFGS == tr->pc_type && tr->M) {
718: MatLMVMGetRejects(tr->M, &nrejects);
719: PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n", nrejects);
720: }
721: PetscViewerASCIIPopTab(viewer);
722: }
723: return(0);
724: }
726: /*------------------------------------------------------------*/
727: /*MC
728: TAONTR - Newton's method with trust region for unconstrained minimization.
729: At each iteration, the Newton trust region method solves the system.
730: NTR expects a KSP solver with a trust region radius.
731: min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k
733: Options Database Keys:
734: + -tao_ntr_ksp_type - "nash","stcg","gltr"
735: . -tao_ntr_pc_type - "none","ahess","bfgs","petsc"
736: . -tao_ntr_bfgs_scale_type - type of scaling with bfgs pc, "ahess" or "bfgs"
737: . -tao_ntr_init_type - "constant","direction","interpolation"
738: . -tao_ntr_update_type - "reduction","interpolation"
739: . -tao_ntr_min_radius - lower bound on trust region radius
740: . -tao_ntr_max_radius - upper bound on trust region radius
741: . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
742: . -tao_ntr_mu1_i - mu1 interpolation init factor
743: . -tao_ntr_mu2_i - mu2 interpolation init factor
744: . -tao_ntr_gamma1_i - gamma1 interpolation init factor
745: . -tao_ntr_gamma2_i - gamma2 interpolation init factor
746: . -tao_ntr_gamma3_i - gamma3 interpolation init factor
747: . -tao_ntr_gamma4_i - gamma4 interpolation init factor
748: . -tao_ntr_theta_i - thetha1 interpolation init factor
749: . -tao_ntr_eta1 - eta1 reduction update factor
750: . -tao_ntr_eta2 - eta2 reduction update factor
751: . -tao_ntr_eta3 - eta3 reduction update factor
752: . -tao_ntr_eta4 - eta4 reduction update factor
753: . -tao_ntr_alpha1 - alpha1 reduction update factor
754: . -tao_ntr_alpha2 - alpha2 reduction update factor
755: . -tao_ntr_alpha3 - alpha3 reduction update factor
756: . -tao_ntr_alpha4 - alpha4 reduction update factor
757: . -tao_ntr_alpha4 - alpha4 reduction update factor
758: . -tao_ntr_mu1 - mu1 interpolation update
759: . -tao_ntr_mu2 - mu2 interpolation update
760: . -tao_ntr_gamma1 - gamma1 interpolcation update
761: . -tao_ntr_gamma2 - gamma2 interpolcation update
762: . -tao_ntr_gamma3 - gamma3 interpolcation update
763: . -tao_ntr_gamma4 - gamma4 interpolation update
764: - -tao_ntr_theta - theta interpolation update
766: Level: beginner
767: M*/
771: PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
772: {
773: TAO_NTR *tr;
778: PetscNewLog(tao,&tr);
780: tao->ops->setup = TaoSetUp_NTR;
781: tao->ops->solve = TaoSolve_NTR;
782: tao->ops->view = TaoView_NTR;
783: tao->ops->setfromoptions = TaoSetFromOptions_NTR;
784: tao->ops->destroy = TaoDestroy_NTR;
786: /* Override default settings (unless already changed) */
787: if (!tao->max_it_changed) tao->max_it = 50;
788: if (!tao->trust0_changed) tao->trust0 = 100.0;
789: #if defined(PETSC_USE_REAL_SINGLE)
790: if (!tao->fatol_changed) tao->fatol = 1.0e-5;
791: if (!tao->frtol_changed) tao->frtol = 1.0e-5;
792: #else
793: if (!tao->fatol_changed) tao->fatol = 1.0e-10;
794: if (!tao->frtol_changed) tao->frtol = 1.0e-10;
795: #endif
796: tao->data = (void*)tr;
798: /* Standard trust region update parameters */
799: tr->eta1 = 1.0e-4;
800: tr->eta2 = 0.25;
801: tr->eta3 = 0.50;
802: tr->eta4 = 0.90;
804: tr->alpha1 = 0.25;
805: tr->alpha2 = 0.50;
806: tr->alpha3 = 1.00;
807: tr->alpha4 = 2.00;
808: tr->alpha5 = 4.00;
810: /* Interpolation parameters */
811: tr->mu1_i = 0.35;
812: tr->mu2_i = 0.50;
814: tr->gamma1_i = 0.0625;
815: tr->gamma2_i = 0.50;
816: tr->gamma3_i = 2.00;
817: tr->gamma4_i = 5.00;
819: tr->theta_i = 0.25;
821: /* Interpolation trust region update parameters */
822: tr->mu1 = 0.10;
823: tr->mu2 = 0.50;
825: tr->gamma1 = 0.25;
826: tr->gamma2 = 0.50;
827: tr->gamma3 = 2.00;
828: tr->gamma4 = 4.00;
830: tr->theta = 0.05;
832: tr->min_radius = 1.0e-10;
833: tr->max_radius = 1.0e10;
834: tr->epsilon = 1.0e-6;
836: tr->ksp_type = NTR_KSP_STCG;
837: tr->pc_type = NTR_PC_BFGS;
838: tr->bfgs_scale_type = BFGS_SCALE_AHESS;
839: tr->init_type = NTR_INIT_INTERPOLATION;
840: tr->update_type = NTR_UPDATE_REDUCTION;
843: /* Set linear solver to default for trust region */
844: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
845: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
846: return(0);
847: }
852: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x)
853: {
855: Mat M;
860: PCShellGetContext(pc,(void**)&M);
861: MatLMVMSolve(M, b, x);
862: return(0);
863: }