Actual source code: ntr.c
1: /*$Id$*/
3: #include "ntr.h"
5: #ifdef TAO_USE_PETSC
6: #include "petscksp.h"
7: #include "petscpc.h"
8: #include "src/petsctao/linearsolver/taolinearsolver_petsc.h"
9: #include "src/petsctao/vector/taovec_petsc.h"
11: #include "private/kspimpl.h"
12: #include "private/pcimpl.h"
14: #define NTR_KSP_NASH 0
15: #define NTR_KSP_STCG 1
16: #define NTR_KSP_GLTR 2
17: #define NTR_KSP_TYPES 3
19: #define NTR_PC_NONE 0
20: #define NTR_PC_AHESS 1
21: #define NTR_PC_BFGS 2
22: #define NTR_PC_PETSC 3
23: #define NTR_PC_TYPES 4
25: #define BFGS_SCALE_AHESS 0
26: #define BFGS_SCALE_BFGS 1
27: #define BFGS_SCALE_TYPES 2
29: #define NTR_INIT_CONSTANT 0
30: #define NTR_INIT_DIRECTION 1
31: #define NTR_INIT_INTERPOLATION 2
32: #define NTR_INIT_TYPES 3
34: #define NTR_UPDATE_REDUCTION 0
35: #define NTR_UPDATE_INTERPOLATION 1
36: #define NTR_UPDATE_TYPES 2
38: static const char *NTR_KSP[64] = {
39: "nash", "stcg", "gltr"
40: };
42: static const char *NTR_PC[64] = {
43: "none", "ahess", "bfgs", "petsc"
44: };
46: static const char *BFGS_SCALE[64] = {
47: "ahess", "bfgs"
48: };
50: static const char *NTR_INIT[64] = {
51: "constant", "direction", "interpolation"
52: };
54: static const char *NTR_UPDATE[64] = {
55: "reduction", "interpolation"
56: };
58: // Routine for BFGS preconditioner
62: static PetscErrorCode bfgs_apply(PC pc, Vec xin, Vec xout)
63: {
64: TaoLMVMMat *M;
65: TaoVecPetsc Xin(xin);
66: TaoVecPetsc Xout(xout);
67: TaoTruth info2;
68: int info;
71: info = PCShellGetContext(pc,(void**)&M); CHKERRQ(info);
72: info = M->Solve(&Xin, &Xout, &info2); CHKERRQ(info);
73: return(0);
74: }
76: // TaoSolve_NTR - Implements Newton's Method with a trust region approach
77: // for solving unconstrained minimization problems.
78: //
79: // The basic algorithm is taken from MINPACK-2 (dstrn).
80: //
81: // TaoSolve_NTR computes a local minimizer of a twice differentiable function
82: // f by applying a trust region variant of Newton's method. At each stage
83: // of the algorithm, we use the prconditioned conjugate gradient method to
84: // determine an approximate minimizer of the quadratic equation
85: //
86: // q(s) = <s, Hs + g>
87: //
88: // subject to the trust region constraint
89: //
90: // || s ||_M <= radius,
91: //
92: // where radius is the trust region radius and M is a symmetric positive
93: // definite matrix (the preconditioner). Here g is the gradient and H
94: // is the Hessian matrix.
95: //
96: // Note: TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG,
97: // or KSPGLTR. Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this
98: // routine regardless of what the user may have previously specified.
102: static int TaoSolve_NTR(TAO_SOLVER tao, void *solver)
103: {
104: TAO_NTR *tr = (TAO_NTR *)solver;
105: TaoVec *X, *G = tr->G, *D = tr->D, *W = tr->W;
106: TaoVec *Diag = tr->Diag;
107: TaoMat *H;
108: TaoLMVMMat *M = tr->M;
110: TaoLinearSolver *tls;
111: TaoLinearSolverPetsc *pls;
113: KSP pksp;
114: PC ppc;
116: KSPConvergedReason ksp_reason;
117: TaoTerminateReason reason;
118: TaoTruth success;
120: double fmin, ftrial, prered, actred, kappa, sigma, beta;
121: double tau, tau_1, tau_2, tau_max, tau_min, max_radius;
122: double f, gnorm;
124: double delta;
125: double radius, norm_d;
126: int info;
128: TaoInt iter = 0;
129: TaoInt bfgsUpdates = 0;
130: TaoInt needH;
132: TaoInt i_max = 5;
133: TaoInt j_max = 1;
134: TaoInt i, j;
136: TaoFunctionBegin;
138: // Get the initial trust-region radius
139: info = TaoGetInitialTrustRegionRadius(tao, &radius); CHKERRQ(info);
140: if (radius < 0.0) {
141: SETERRQ(1, "Initial radius negative");
142: }
144: // Modify the radius if it is too large or small
145: radius = TaoMax(radius, tr->min_radius);
146: radius = TaoMin(radius, tr->max_radius);
148: // Get vectors we will need
149: info = TaoGetSolution(tao, &X); CHKERRQ(info);
150: info = TaoGetHessian(tao, &H); CHKERRQ(info);
152: if (NTR_PC_BFGS == tr->pc_type && !M) {
153: tr->M = new TaoLMVMMat(X);
154: M = tr->M;
155: }
157: // Check convergence criteria
158: info = TaoComputeFunctionGradient(tao, X, &f, G); CHKERRQ(info);
159: info = G->Norm2(&gnorm); CHKERRQ(info);
160: if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
161: SETERRQ(1, "User provided compute function generated Inf or NaN");
162: }
163: needH = 1;
165: info = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); CHKERRQ(info);
166: if (reason != TAO_CONTINUE_ITERATING) {
167: TaoFunctionReturn(0);
168: }
170: // Create vectors for the limited memory preconditioner
171: if ((NTR_PC_BFGS == tr->pc_type) &&
172: (BFGS_SCALE_BFGS != tr->bfgs_scale_type)) {
173: if (!Diag) {
174: info = X->Clone(&tr->Diag); CHKERRQ(info);
175: Diag = tr->Diag;
176: }
177: }
179: // Modify the linear solver to a conjugate gradient method
180: info = TaoGetLinearSolver(tao, &tls); CHKERRQ(info);
181: pls = dynamic_cast <TaoLinearSolverPetsc *> (tls);
183: pksp = pls->GetKSP();
184: switch(tr->ksp_type) {
185: case NTR_KSP_NASH:
186: info = KSPSetType(pksp, KSPNASH); CHKERRQ(info);
187: if (pksp->ops->setfromoptions) {
188: (*pksp->ops->setfromoptions)(pksp);
189: }
190: break;
192: case NTR_KSP_STCG:
193: info = KSPSetType(pksp, KSPSTCG); CHKERRQ(info);
194: if (pksp->ops->setfromoptions) {
195: (*pksp->ops->setfromoptions)(pksp);
196: }
197: break;
199: default:
200: info = KSPSetType(pksp, KSPGLTR); CHKERRQ(info);
201: if (pksp->ops->setfromoptions) {
202: (*pksp->ops->setfromoptions)(pksp);
203: }
204: break;
205: }
207: // Modify the preconditioner to use the bfgs approximation
208: info = KSPGetPC(pksp, &ppc); CHKERRQ(info);
209: switch(tr->pc_type) {
210: case NTR_PC_NONE:
211: info = PCSetType(ppc, PCNONE); CHKERRQ(info);
212: if (ppc->ops->setfromoptions) {
213: (*ppc->ops->setfromoptions)(ppc);
214: }
215: break;
216:
217: case NTR_PC_AHESS:
218: info = PCSetType(ppc, PCJACOBI); CHKERRQ(info);
219: if (ppc->ops->setfromoptions) {
220: (*ppc->ops->setfromoptions)(ppc);
221: }
222: info = PCJacobiSetUseAbs(ppc); CHKERRQ(info);
223: break;
225: case NTR_PC_BFGS:
226: info = PCSetType(ppc, PCSHELL); CHKERRQ(info);
227: if (ppc->ops->setfromoptions) {
228: (*ppc->ops->setfromoptions)(ppc);
229: }
230: info = PCShellSetName(ppc, "bfgs"); CHKERRQ(info);
231: info = PCShellSetContext(ppc, M); CHKERRQ(info);
232: info = PCShellSetApply(ppc, bfgs_apply); CHKERRQ(info);
233: break;
235: default:
236: // Use the pc method set by pc_type
237: break;
238: }
240: // Initialize trust-region radius
241: switch(tr->init_type) {
242: case NTR_INIT_CONSTANT:
243: // Use the initial radius specified
244: break;
246: case NTR_INIT_INTERPOLATION:
247: // Use the initial radius specified
248: max_radius = 0.0;
250: for (j = 0; j < j_max; ++j) {
251: fmin = f;
252: sigma = 0.0;
254: if (needH) {
255: info = TaoComputeHessian(tao, X, H); CHKERRQ(info);
256: needH = 0;
257: }
259: for (i = 0; i < i_max; ++i) {
260: info = W->Waxpby(1.0, X, -radius / gnorm, G); CHKERRQ(info);
262: info = TaoComputeFunction(tao, W, &ftrial); CHKERRQ(info);
263: if (TaoInfOrNaN(ftrial)) {
264: tau = tr->gamma1_i;
265: }
266: else {
267: if (ftrial < fmin) {
268: fmin = ftrial;
269: sigma = -radius / gnorm;
270: }
272: info = H->Multiply(G, D); CHKERRQ(info);
273: info = D->Dot(G, &prered); CHKERRQ(info);
275: prered = radius * (gnorm - 0.5 * radius * prered / (gnorm * gnorm));
276: actred = f - ftrial;
277: if ((fabs(actred) <= tr->epsilon) &&
278: (fabs(prered) <= tr->epsilon)) {
279: kappa = 1.0;
280: }
281: else {
282: kappa = actred / prered;
283: }
285: tau_1 = tr->theta_i * gnorm * radius / (tr->theta_i * gnorm * radius + (1.0 - tr->theta_i) * prered - actred);
286: tau_2 = tr->theta_i * gnorm * radius / (tr->theta_i * gnorm * radius - (1.0 + tr->theta_i) * prered + actred);
287: tau_min = TaoMin(tau_1, tau_2);
288: tau_max = TaoMax(tau_1, tau_2);
290: if (fabs(kappa - 1.0) <= tr->mu1_i) {
291: // Great agreement
292: max_radius = TaoMax(max_radius, radius);
294: if (tau_max < 1.0) {
295: tau = tr->gamma3_i;
296: }
297: else if (tau_max > tr->gamma4_i) {
298: tau = tr->gamma4_i;
299: }
300: else {
301: tau = tau_max;
302: }
303: }
304: else if (fabs(kappa - 1.0) <= tr->mu2_i) {
305: // Good agreement
306: max_radius = TaoMax(max_radius, radius);
308: if (tau_max < tr->gamma2_i) {
309: tau = tr->gamma2_i;
310: }
311: else if (tau_max > tr->gamma3_i) {
312: tau = tr->gamma3_i;
313: }
314: else {
315: tau = tau_max;
316: }
317: }
318: else {
319: // Not good agreement
320: if (tau_min > 1.0) {
321: tau = tr->gamma2_i;
322: }
323: else if (tau_max < tr->gamma1_i) {
324: tau = tr->gamma1_i;
325: }
326: else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
327: tau = tr->gamma1_i;
328: }
329: else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) &&
330: ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
331: tau = tau_1;
332: }
333: else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) &&
334: ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
335: tau = tau_2;
336: }
337: else {
338: tau = tau_max;
339: }
340: }
341: }
342: radius = tau * radius;
343: }
344:
345: if (fmin < f) {
346: f = fmin;
347: info = X->Axpy(sigma, G); CHKERRQ(info);
348: info = TaoComputeGradient(tao, X, G); CHKERRQ(info);
350: info = G->Norm2(&gnorm); CHKERRQ(info);
351: if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
352: SETERRQ(1, "User provided compute function generated Inf or NaN");
353: }
354: needH = 1;
356: info = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); CHKERRQ(info);
357: if (reason != TAO_CONTINUE_ITERATING) {
358: TaoFunctionReturn(0);
359: }
360: }
361: }
362: radius = TaoMax(radius, max_radius);
364: // Modify the radius if it is too large or small
365: radius = TaoMax(radius, tr->min_radius);
366: radius = TaoMin(radius, tr->max_radius);
367: break;
369: default:
370: // Norm of the first direction will initialize radius
371: radius = 0.0;
372: break;
373: }
375: // Set initial scaling for the BFGS preconditioner
376: // This step is done after computing the initial trust-region radius
377: // since the function value may have decreased
378: if (NTR_PC_BFGS == tr->pc_type) {
379: if (f != 0.0) {
380: delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
381: }
382: else {
383: delta = 2.0 / (gnorm*gnorm);
384: }
385: info = M->SetDelta(delta); CHKERRQ(info);
386: }
388: // Have not converged; continue with Newton method
389: while (reason == TAO_CONTINUE_ITERATING) {
390: ++iter;
392: // Compute the Hessian
393: if (needH) {
394: info = TaoComputeHessian(tao, X, H); CHKERRQ(info);
395: needH = 0;
396: }
398: if (NTR_PC_BFGS == tr->pc_type) {
399: if (BFGS_SCALE_AHESS == tr->bfgs_scale_type) {
400: // Obtain diagonal for the bfgs preconditioner
401: info = H->GetDiagonal(Diag); CHKERRQ(info);
402: info = Diag->AbsoluteValue(); CHKERRQ(info);
403: info = Diag->Reciprocal(); CHKERRQ(info);
404: info = M->SetScale(Diag); CHKERRQ(info);
405: }
407: // Update the limited memory preconditioner
408: info = M->Update(X, G); CHKERRQ(info);
409: ++bfgsUpdates;
410: }
412: while (reason == TAO_CONTINUE_ITERATING) {
413: // Solve the trust region subproblem
414: info = TaoPreLinearSolve(tao, H); CHKERRQ(info);
415: info = TaoLinearSolveTrustRegion(tao, H, G, D, radius, &success); CHKERRQ(info);
416: info = pls->GetNormDirection(&norm_d); CHKERRQ(info);
417: if (0.0 == radius) {
418: // Radius was uninitialized; use the norm of the direction
419: if (norm_d > 0.0) {
420: radius = norm_d;
422: // Modify the radius if it is too large or small
423: radius = TaoMax(radius, tr->min_radius);
424: radius = TaoMin(radius, tr->max_radius);
425: }
426: else {
427: // The direction was bad; set radius to default value and re-solve
428: // the trust-region subproblem to get a direction
429: info = TaoGetInitialTrustRegionRadius(tao, &radius); CHKERRQ(info);
431: // Modify the radius if it is too large or small
432: radius = TaoMax(radius, tr->min_radius);
433: radius = TaoMin(radius, tr->max_radius);
435: info = TaoLinearSolveTrustRegion(tao, H, G, D, radius, &success); CHKERRQ(info);
436: info = pls->GetNormDirection(&norm_d); CHKERRQ(info);
437: if (norm_d == 0.0) {
438: SETERRQ(1, "Initial direction zero");
439: }
440: }
441: }
442: info = D->Negate(); CHKERRQ(info);
444: info = KSPGetConvergedReason(pksp, &ksp_reason); CHKERRQ(info);
445: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
446: (NTR_PC_BFGS == tr->pc_type) && (bfgsUpdates > 1)) {
447: // Preconditioner is numerically indefinite; reset the
448: // approximate if using BFGS preconditioning.
449:
450: if (f != 0.0) {
451: delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
452: }
453: else {
454: delta = 2.0 / (gnorm*gnorm);
455: }
456: info = M->SetDelta(delta); CHKERRQ(info);
457: info = M->Reset(); CHKERRQ(info);
458: info = M->Update(X, G); CHKERRQ(info);
459: bfgsUpdates = 1;
460: }
462: if (NTR_UPDATE_REDUCTION == tr->update_type) {
463: // Get predicted reduction
464: info = pls->GetObjFcn(&prered); CHKERRQ(info);
465:
466: if (prered >= 0.0) {
467: // The predicted reduction has the wrong sign. This cannot
468: // happen in infinite precision arithmetic. Step should
469: // be rejected!
470: radius = tr->alpha1 * TaoMin(radius, norm_d);
471: }
472: else {
473: // Compute trial step and function value
474: info = W->Waxpby(1.0, X, 1.0, D); CHKERRQ(info);
475: info = TaoComputeFunction(tao, W, &ftrial); CHKERRQ(info);
476: if (TaoInfOrNaN(ftrial)) {
477: radius = tr->alpha1 * TaoMin(radius, norm_d);
478: }
479: else {
480: // Compute and actual reduction
481: actred = f - ftrial;
482: prered = -prered;
483: if ((fabs(actred) <= tr->epsilon) &&
484: (fabs(prered) <= tr->epsilon)) {
485: kappa = 1.0;
486: }
487: else {
488: kappa = actred / prered;
489: }
490:
491: // Accept or reject the step and update radius
492: if (kappa < tr->eta1) {
493: // Reject the step
494: radius = tr->alpha1 * TaoMin(radius, norm_d);
495: }
496: else {
497: // Accept the step
498: if (kappa < tr->eta2) {
499: // Marginal bad step
500: radius = tr->alpha2 * TaoMin(radius, norm_d);
501: }
502: else if (kappa < tr->eta3) {
503: // Reasonable step
504: radius = tr->alpha3 * radius;
505: }
506: else if (kappa < tr->eta4) {
507: // Good step
508: radius = TaoMax(tr->alpha4 * norm_d, radius);
509: }
510: else {
511: // Very good step
512: radius = TaoMax(tr->alpha5 * norm_d, radius);
513: }
514: break;
515: }
516: }
517: }
518: }
519: else {
520: // Get predicted reduction
521: info = pls->GetObjFcn(&prered); CHKERRQ(info);
523: if (prered >= 0.0) {
524: // The predicted reduction has the wrong sign. This cannot
525: // happen in infinite precision arithmetic. Step should
526: // be rejected!
527: radius = tr->gamma1 * TaoMin(radius, norm_d);
528: }
529: else {
530: info = W->Waxpby(1.0, X, 1.0, D); CHKERRQ(info);
531: info = TaoComputeFunction(tao, W, &ftrial); CHKERRQ(info);
532: if (TaoInfOrNaN(ftrial)) {
533: radius = tr->gamma1 * TaoMin(radius, norm_d);
534: }
535: else {
536: info = D->Dot(G, &beta); CHKERRQ(info);
538: actred = f - ftrial;
539: prered = -prered;
540: if ((fabs(actred) <= tr->epsilon) &&
541: (fabs(prered) <= tr->epsilon)) {
542: kappa = 1.0;
543: }
544: else {
545: kappa = actred / prered;
546: }
548: tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
549: tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
550: tau_min = TaoMin(tau_1, tau_2);
551: tau_max = TaoMax(tau_1, tau_2);
553: if (kappa >= 1.0 - tr->mu1) {
554: // Great agreement; accept step and update radius
555: if (tau_max < 1.0) {
556: radius = TaoMax(radius, tr->gamma3 * norm_d);
557: }
558: else if (tau_max > tr->gamma4) {
559: radius = TaoMax(radius, tr->gamma4 * norm_d);
560: }
561: else {
562: radius = TaoMax(radius, tau_max * norm_d);
563: }
564: break;
565: }
566: else if (kappa >= 1.0 - tr->mu2) {
567: // Good agreement
569: if (tau_max < tr->gamma2) {
570: radius = tr->gamma2 * TaoMin(radius, norm_d);
571: }
572: else if (tau_max > tr->gamma3) {
573: radius = TaoMax(radius, tr->gamma3 * norm_d);
574: }
575: else if (tau_max < 1.0) {
576: radius = tau_max * TaoMin(radius, norm_d);
577: }
578: else {
579: radius = TaoMax(radius, tau_max * norm_d);
580: }
581: break;
582: }
583: else {
584: // Not good agreement
585: if (tau_min > 1.0) {
586: radius = tr->gamma2 * TaoMin(radius, norm_d);
587: }
588: else if (tau_max < tr->gamma1) {
589: radius = tr->gamma1 * TaoMin(radius, norm_d);
590: }
591: else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
592: radius = tr->gamma1 * TaoMin(radius, norm_d);
593: }
594: else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) &&
595: ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
596: radius = tau_1 * TaoMin(radius, norm_d);
597: }
598: else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) &&
599: ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
600: radius = tau_2 * TaoMin(radius, norm_d);
601: }
602: else {
603: radius = tau_max * TaoMin(radius, norm_d);
604: }
605: }
606: }
607: }
608: }
610: // The step computed was not good and the radius was decreased.
611: // Monitor the radius to terminate.
612: info = TaoMonitor(tao, iter, f, gnorm, 0.0, radius, &reason); CHKERRQ(info);
613: }
615: // The radius may have been increased; modify if it is too large
616: radius = TaoMin(radius, tr->max_radius);
618: if (reason == TAO_CONTINUE_ITERATING) {
619: info = X->CopyFrom(W); CHKERRQ(info);
620: f = ftrial;
621: info = TaoComputeGradient(tao, X, G); CHKERRQ(info);
622: info = G->Norm2(&gnorm); CHKERRQ(info);
623: if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
624: SETERRQ(1, "User provided compute function generated Inf or NaN");
625: }
626: needH = 1;
628: info = TaoMonitor(tao, iter, f, gnorm, 0.0, radius, &reason); CHKERRQ(info);
629: }
630: }
631: TaoFunctionReturn(0);
632: }
634: /*------------------------------------------------------------*/
637: static int TaoSetUp_NTR(TAO_SOLVER tao, void *solver)
638: {
639: TAO_NTR *tr = (TAO_NTR *)solver;
640: TaoVec *X;
641: TaoMat *H;
642: int info;
644: TaoFunctionBegin;
646: info = TaoGetSolution(tao, &X); CHKERRQ(info);
647: info = X->Clone(&tr->G); CHKERRQ(info);
648: info = X->Clone(&tr->D); CHKERRQ(info);
649: info = X->Clone(&tr->W); CHKERRQ(info);
651: tr->Diag = 0;
652: tr->M = 0;
654: info = TaoSetLagrangianGradientVector(tao, tr->G); CHKERRQ(info);
655: info = TaoSetStepDirectionVector(tao, tr->D); CHKERRQ(info);
657: // Set linear solver to default for trust region
658: info = TaoGetHessian(tao, &H); CHKERRQ(info);
659: info = TaoCreateLinearSolver(tao, H, 200, 0); CHKERRQ(info);
661: // Check sizes for compatability
662: info = TaoCheckFGH(tao); CHKERRQ(info);
663: TaoFunctionReturn(0);
664: }
666: /*------------------------------------------------------------*/
669: static int TaoDestroy_NTR(TAO_SOLVER tao, void *solver)
670: {
671: TAO_NTR *tr = (TAO_NTR *)solver;
672: int info;
674: TaoFunctionBegin;
675: info = TaoVecDestroy(tr->G); CHKERRQ(info);
676: info = TaoVecDestroy(tr->D); CHKERRQ(info);
677: info = TaoVecDestroy(tr->W); CHKERRQ(info);
679: info = TaoSetLagrangianGradientVector(tao, 0); CHKERRQ(info);
680: info = TaoSetStepDirectionVector(tao, 0); CHKERRQ(info);
682: if (tr->Diag) {
683: info = TaoVecDestroy(tr->Diag); CHKERRQ(info);
684: tr->Diag = 0;
685: }
687: if (tr->M) {
688: info = TaoMatDestroy(tr->M); CHKERRQ(info);
689: tr->M = 0;
690: }
691: TaoFunctionReturn(0);
692: }
694: /*------------------------------------------------------------*/
697: static int TaoSetOptions_NTR(TAO_SOLVER tao, void*solver)
698: {
699: TAO_NTR *tr = (TAO_NTR *)solver;
700: int info;
702: TaoFunctionBegin;
703: info = TaoOptionsHead("Newton trust region method for unconstrained optimization"); CHKERRQ(info);
704: info = TaoOptionList("-tao_ntr_ksp_type", "ksp type", "", NTR_KSP, NTR_KSP_TYPES, NTR_KSP[tr->ksp_type], &tr->ksp_type, 0); CHKERRQ(info);
705: info = TaoOptionList("-tao_ntr_pc_type", "pc type", "", NTR_PC, NTR_PC_TYPES, NTR_PC[tr->pc_type], &tr->pc_type, 0); CHKERRQ(info);
706: info = TaoOptionList("-tao_ntr_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tr->bfgs_scale_type], &tr->bfgs_scale_type, 0); CHKERRQ(info);
707: info = TaoOptionList("-tao_ntr_init_type", "radius initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type, 0); CHKERRQ(info);
708: info = TaoOptionList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type, 0); CHKERRQ(info);
709: info = TaoOptionDouble("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1, 0); CHKERRQ(info);
710: info = TaoOptionDouble("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2, 0); CHKERRQ(info);
711: info = TaoOptionDouble("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3, 0); CHKERRQ(info);
712: info = TaoOptionDouble("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4, 0); CHKERRQ(info);
713: info = TaoOptionDouble("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1, 0); CHKERRQ(info);
714: info = TaoOptionDouble("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2, 0); CHKERRQ(info);
715: info = TaoOptionDouble("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3, 0); CHKERRQ(info);
716: info = TaoOptionDouble("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4, 0); CHKERRQ(info);
717: info = TaoOptionDouble("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5, 0); CHKERRQ(info);
718: info = TaoOptionDouble("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1, 0); CHKERRQ(info);
719: info = TaoOptionDouble("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2, 0); CHKERRQ(info);
720: info = TaoOptionDouble("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1, 0); CHKERRQ(info);
721: info = TaoOptionDouble("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2, 0); CHKERRQ(info);
722: info = TaoOptionDouble("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3, 0); CHKERRQ(info);
723: info = TaoOptionDouble("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4, 0); CHKERRQ(info);
724: info = TaoOptionDouble("-tao_ntr_theta", "", "", tr->theta, &tr->theta, 0); CHKERRQ(info);
725: info = TaoOptionDouble("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i, 0); CHKERRQ(info);
726: info = TaoOptionDouble("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i, 0); CHKERRQ(info);
727: info = TaoOptionDouble("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i, 0); CHKERRQ(info);
728: info = TaoOptionDouble("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i, 0); CHKERRQ(info);
729: info = TaoOptionDouble("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i, 0); CHKERRQ(info);
730: info = TaoOptionDouble("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i, 0); CHKERRQ(info);
731: info = TaoOptionDouble("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i, 0); CHKERRQ(info);
732: info = TaoOptionDouble("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius, 0); CHKERRQ(info);
733: info = TaoOptionDouble("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius, 0); CHKERRQ(info);
734: info = TaoOptionDouble("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon, 0); CHKERRQ(info);
735: info = TaoOptionsTail(); CHKERRQ(info);
736: TaoFunctionReturn(0);
737: }
739: /*------------------------------------------------------------*/
742: static int TaoView_NTR(TAO_SOLVER tao,void*solver)
743: {
744: TAO_NTR *tr = (TAO_NTR *)solver;
745: int info;
747: TaoFunctionBegin;
748: if (NTR_PC_BFGS == tr->pc_type && tr->M) {
749: info = TaoPrintInt(tao, " Rejected matrix updates: %d\n", tr->M->GetRejects()); CHKERRQ(info);
750: }
751: TaoFunctionReturn(0);
752: }
754: /*------------------------------------------------------------*/
758: int TaoCreate_NTR(TAO_SOLVER tao)
759: {
760: TAO_NTR *tr;
761: int info;
763: TaoFunctionBegin;
765: info = TaoNew(TAO_NTR, &tr); CHKERRQ(info);
766: info = PetscLogObjectMemory(tao, sizeof(TAO_NTR)); CHKERRQ(info);
768: info = TaoSetTaoSolveRoutine(tao, TaoSolve_NTR, (void *)tr); CHKERRQ(info);
769: info = TaoSetTaoSetUpDownRoutines(tao, TaoSetUp_NTR, TaoDestroy_NTR); CHKERRQ(info);
770: info = TaoSetTaoOptionsRoutine(tao, TaoSetOptions_NTR); CHKERRQ(info);
771: info = TaoSetTaoViewRoutine(tao, TaoView_NTR); CHKERRQ(info);
773: info = TaoSetMaximumIterates(tao, 50); CHKERRQ(info);
774: info = TaoSetTolerances(tao, 1e-10, 1e-10, 0, 0); CHKERRQ(info);
776: info = TaoSetTrustRegionRadius(tao, 100.0); CHKERRQ(info);
777: info = TaoSetTrustRegionTolerance(tao, 1.0e-12); CHKERRQ(info);
779: // Standard trust region update parameters
780: tr->eta1 = 1.0e-4;
781: tr->eta2 = 0.25;
782: tr->eta3 = 0.50;
783: tr->eta4 = 0.90;
785: tr->alpha1 = 0.25;
786: tr->alpha2 = 0.50;
787: tr->alpha3 = 1.00;
788: tr->alpha4 = 2.00;
789: tr->alpha5 = 4.00;
791: // Interpolation parameters
792: tr->mu1_i = 0.35;
793: tr->mu2_i = 0.50;
795: tr->gamma1_i = 0.0625;
796: tr->gamma2_i = 0.50;
797: tr->gamma3_i = 2.00;
798: tr->gamma4_i = 5.00;
800: tr->theta_i = 0.25;
802: // Interpolation trust region update parameters
803: tr->mu1 = 0.10;
804: tr->mu2 = 0.50;
806: tr->gamma1 = 0.25;
807: tr->gamma2 = 0.50;
808: tr->gamma3 = 2.00;
809: tr->gamma4 = 4.00;
811: tr->theta = 0.05;
813: tr->min_radius = 1.0e-10;
814: tr->max_radius = 1.0e10;
815: tr->epsilon = 1.0e-6;
817: tr->ksp_type = NTR_KSP_STCG;
818: tr->pc_type = NTR_PC_BFGS;
819: tr->bfgs_scale_type = BFGS_SCALE_AHESS;
820: tr->init_type = NTR_INIT_INTERPOLATION;
821: tr->update_type = NTR_UPDATE_REDUCTION;
822: TaoFunctionReturn(0);
823: }
826: #endif