Actual source code: nls.c
1: /*$Id$*/
3: #include "nls.h" /*I "tao_solver.h" I*/
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 NLS_KSP_CG 0
15: #define NLS_KSP_NASH 1
16: #define NLS_KSP_STCG 2
17: #define NLS_KSP_GLTR 3
18: #define NLS_KSP_PETSC 4
19: #define NLS_KSP_TYPES 5
21: #define NLS_PC_NONE 0
22: #define NLS_PC_AHESS 1
23: #define NLS_PC_BFGS 2
24: #define NLS_PC_PETSC 3
25: #define NLS_PC_TYPES 4
27: #define BFGS_SCALE_AHESS 0
28: #define BFGS_SCALE_PHESS 1
29: #define BFGS_SCALE_BFGS 2
30: #define BFGS_SCALE_TYPES 3
32: #define NLS_INIT_CONSTANT 0
33: #define NLS_INIT_DIRECTION 1
34: #define NLS_INIT_INTERPOLATION 2
35: #define NLS_INIT_TYPES 3
37: #define NLS_UPDATE_STEP 0
38: #define NLS_UPDATE_REDUCTION 1
39: #define NLS_UPDATE_INTERPOLATION 2
40: #define NLS_UPDATE_TYPES 3
42: static const char *NLS_KSP[64] = {
43: "cg", "nash", "stcg", "gltr", "petsc"
44: };
46: static const char *NLS_PC[64] = {
47: "none", "ahess", "bfgs", "petsc"
48: };
50: static const char *BFGS_SCALE[64] = {
51: "ahess", "phess", "bfgs"
52: };
54: static const char *NLS_INIT[64] = {
55: "constant", "direction", "interpolation"
56: };
58: static const char *NLS_UPDATE[64] = {
59: "step", "reduction", "interpolation"
60: };
62: // Routine for BFGS preconditioner
66: static PetscErrorCode bfgs_apply(PC pc, Vec xin, Vec xout)
67: {
68: TaoLMVMMat *M ;
69: TaoVecPetsc Xin(xin);
70: TaoVecPetsc Xout(xout);
71: TaoTruth info2;
72: int info;
75: info = PCShellGetContext(pc,(void**)&M); CHKERRQ(info);
76: info = M->Solve(&Xin, &Xout, &info2); CHKERRQ(info);
77: return(0);
78: }
80: // Implements Newton's Method with a line search approach for solving
81: // unconstrained minimization problems. A More'-Thuente line search
82: // is used to guarantee that the bfgs preconditioner remains positive
83: // definite.
84: //
85: // The method can shift the Hessian matrix. The shifting procedure is
86: // adapted from the PATH algorithm for solving complementarity
87: // problems.
88: //
89: // The linear system solve should be done with a conjugate gradient
90: // method, although any method can be used.
92: #define NLS_NEWTON 0
93: #define NLS_BFGS 1
94: #define NLS_SCALED_GRADIENT 2
95: #define NLS_GRADIENT 3
99: static int TaoSolve_NLS(TAO_SOLVER tao, void *solver)
100: {
101: TAO_NLS *ls = (TAO_NLS *)solver;
102: TaoVec *X, *G = ls->G, *D = ls->D, *W = ls->W;
103: TaoVec *Xold = ls->Xold, *Gold = ls->Gold, *Diag = ls->Diag;
104: TaoMat *H;
105: TaoLMVMMat *M = ls->M;
107: TaoLinearSolver *tls;
108: TaoLinearSolverPetsc *pls;
110: KSP pksp;
111: PC ppc;
113: KSPConvergedReason ksp_reason;
114: TaoTerminateReason reason;
115: TaoTruth success;
116:
117: double fmin, ftrial, f_full, prered, actred, kappa, sigma;
118: double tau, tau_1, tau_2, tau_max, tau_min, max_radius;
119: double f, fold, gdx, gnorm, pert;
120: double step = 1.0;
122: double delta;
123: double radius, norm_d = 0.0, e_min;
125: int info;
126: TaoInt stepType;
127: TaoInt iter = 0, status = 0;
128: TaoInt bfgsUpdates = 0;
129: TaoInt needH;
131: TaoInt i_max = 5;
132: TaoInt j_max = 1;
133: TaoInt i, j;
135: TaoFunctionBegin;
136: // Initialized variables
137: pert = ls->sval;
139: ls->ksp_atol = 0;
140: ls->ksp_rtol = 0;
141: ls->ksp_dtol = 0;
142: ls->ksp_ctol = 0;
143: ls->ksp_negc = 0;
144: ls->ksp_iter = 0;
145: ls->ksp_othr = 0;
147: // Initialize trust-region radius when using nash, stcg, or gltr
148: // Will be reset during the first iteration
149: if (NLS_KSP_NASH == ls->ksp_type ||
150: NLS_KSP_STCG == ls->ksp_type ||
151: NLS_KSP_GLTR == ls->ksp_type) {
152: info = TaoGetInitialTrustRegionRadius(tao, &radius); CHKERRQ(info);
153: if (radius < 0.0) {
154: SETERRQ(1, "Initial radius negative");
155: }
157: // Modify the radius if it is too large or small
158: radius = TaoMax(radius, ls->min_radius);
159: radius = TaoMin(radius, ls->max_radius);
160: }
162: // Get vectors we will need
163: info = TaoGetSolution(tao, &X); CHKERRQ(info);
164: info = TaoGetHessian(tao, &H); CHKERRQ(info);
166: if (NLS_PC_BFGS == ls->pc_type && !M) {
167: ls->M = new TaoLMVMMat(X);
168: M = ls->M;
169: }
171: // Check convergence criteria
172: info = TaoComputeFunctionGradient(tao, X, &f, G); CHKERRQ(info);
173: info = G->Norm2(&gnorm); CHKERRQ(info);
174: if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
175: SETERRQ(1, "User provided compute function generated Inf or NaN");
176: }
177: needH = 1;
179: info = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); CHKERRQ(info);
180: if (reason != TAO_CONTINUE_ITERATING) {
181: TaoFunctionReturn(0);
182: }
184: // Create vectors for the limited memory preconditioner
185: if ((NLS_PC_BFGS == ls->pc_type) &&
186: (BFGS_SCALE_BFGS != ls->bfgs_scale_type)) {
187: if (!Diag) {
188: info = X->Clone(&ls->Diag); CHKERRQ(info);
189: Diag = ls->Diag;
190: }
191: }
193: // Modify the linear solver to a conjugate gradient method
194: info = TaoGetLinearSolver(tao, &tls); CHKERRQ(info);
195: pls = dynamic_cast <TaoLinearSolverPetsc *> (tls);
197: pksp = pls->GetKSP();
198: switch(ls->ksp_type) {
199: case NLS_KSP_CG:
200: info = KSPSetType(pksp, KSPCG); CHKERRQ(info);
201: if (pksp->ops->setfromoptions) {
202: (*pksp->ops->setfromoptions)(pksp);
203: }
204: break;
206: case NLS_KSP_NASH:
207: info = KSPSetType(pksp, KSPNASH); CHKERRQ(info);
208: if (pksp->ops->setfromoptions) {
209: (*pksp->ops->setfromoptions)(pksp);
210: }
211: break;
213: case NLS_KSP_STCG:
214: info = KSPSetType(pksp, KSPSTCG); CHKERRQ(info);
215: if (pksp->ops->setfromoptions) {
216: (*pksp->ops->setfromoptions)(pksp);
217: }
218: break;
220: case NLS_KSP_GLTR:
221: info = KSPSetType(pksp, KSPGLTR); CHKERRQ(info);
222: if (pksp->ops->setfromoptions) {
223: (*pksp->ops->setfromoptions)(pksp);
224: }
225: break;
227: default:
228: // Use the method set by the ksp_type
229: break;
230: }
232: // Modify the preconditioner to use the bfgs approximation
233: info = KSPGetPC(pksp, &ppc); CHKERRQ(info);
234: switch(ls->pc_type) {
235: case NLS_PC_NONE:
236: info = PCSetType(ppc, PCNONE); CHKERRQ(info);
237: if (ppc->ops->setfromoptions) {
238: (*ppc->ops->setfromoptions)(ppc);
239: }
240: break;
242: case NLS_PC_AHESS:
243: info = PCSetType(ppc, PCJACOBI); CHKERRQ(info);
244: if (ppc->ops->setfromoptions) {
245: (*ppc->ops->setfromoptions)(ppc);
246: }
247: info = PCJacobiSetUseAbs(ppc); CHKERRQ(info);
248: break;
250: case NLS_PC_BFGS:
251: info = PCSetType(ppc, PCSHELL); CHKERRQ(info);
252: if (ppc->ops->setfromoptions) {
253: (*ppc->ops->setfromoptions)(ppc);
254: }
255: info = PCShellSetName(ppc, "bfgs"); CHKERRQ(info);
256: info = PCShellSetContext(ppc, M); CHKERRQ(info);
257: info = PCShellSetApply(ppc, bfgs_apply); CHKERRQ(info);
258: break;
260: default:
261: // Use the pc method set by pc_type
262: break;
263: }
265: // Initialize trust-region radius. The initialization is only performed
266: // when we are using Steihaug-Toint or the Generalized Lanczos method.
267: if (NLS_KSP_NASH == ls->ksp_type ||
268: NLS_KSP_STCG == ls->ksp_type ||
269: NLS_KSP_GLTR == ls->ksp_type) {
270: switch(ls->init_type) {
271: case NLS_INIT_CONSTANT:
272: // Use the initial radius specified
273: break;
275: case NLS_INIT_INTERPOLATION:
276: // Use the initial radius specified
277: max_radius = 0.0;
278:
279: for (j = 0; j < j_max; ++j) {
280: fmin = f;
281: sigma = 0.0;
282:
283: if (needH) {
284: info = TaoComputeHessian(tao, X, H); CHKERRQ(info);
285: needH = 0;
286: }
287:
288: for (i = 0; i < i_max; ++i) {
289: info = W->Waxpby(1.0, X, -radius / gnorm, G); CHKERRQ(info);
290:
291: info = TaoComputeFunction(tao, W, &ftrial); CHKERRQ(info);
292: if (TaoInfOrNaN(ftrial)) {
293: tau = ls->gamma1_i;
294: }
295: else {
296: if (ftrial < fmin) {
297: fmin = ftrial;
298: sigma = -radius / gnorm;
299: }
300:
301: info = H->Multiply(G, D); CHKERRQ(info);
302: info = D->Dot(G, &prered); CHKERRQ(info);
303:
304: prered = radius * (gnorm - 0.5 * radius * prered / (gnorm * gnorm));
305: actred = f - ftrial;
306: if ((fabs(actred) <= ls->epsilon) &&
307: (fabs(prered) <= ls->epsilon)) {
308: kappa = 1.0;
309: }
310: else {
311: kappa = actred / prered;
312: }
313:
314: tau_1 = ls->theta_i * gnorm * radius / (ls->theta_i * gnorm * radius + (1.0 - ls->theta_i) * prered - actred);
315: tau_2 = ls->theta_i * gnorm * radius / (ls->theta_i * gnorm * radius - (1.0 + ls->theta_i) * prered + actred);
316: tau_min = TaoMin(tau_1, tau_2);
317: tau_max = TaoMax(tau_1, tau_2);
318:
319: if (fabs(kappa - 1.0) <= ls->mu1_i) {
320: // Great agreement
321: max_radius = TaoMax(max_radius, radius);
322:
323: if (tau_max < 1.0) {
324: tau = ls->gamma3_i;
325: }
326: else if (tau_max > ls->gamma4_i) {
327: tau = ls->gamma4_i;
328: }
329: else if (tau_1 >= 1.0 && tau_1 <= ls->gamma4_i && tau_2 < 1.0) {
330: tau = tau_1;
331: }
332: else if (tau_2 >= 1.0 && tau_2 <= ls->gamma4_i && tau_1 < 1.0) {
333: tau = tau_2;
334: }
335: else {
336: tau = tau_max;
337: }
338: }
339: else if (fabs(kappa - 1.0) <= ls->mu2_i) {
340: // Good agreement
341: max_radius = TaoMax(max_radius, radius);
342:
343: if (tau_max < ls->gamma2_i) {
344: tau = ls->gamma2_i;
345: }
346: else if (tau_max > ls->gamma3_i) {
347: tau = ls->gamma3_i;
348: }
349: else {
350: tau = tau_max;
351: }
352: }
353: else {
354: // Not good agreement
355: if (tau_min > 1.0) {
356: tau = ls->gamma2_i;
357: }
358: else if (tau_max < ls->gamma1_i) {
359: tau = ls->gamma1_i;
360: }
361: else if ((tau_min < ls->gamma1_i) && (tau_max >= 1.0)) {
362: tau = ls->gamma1_i;
363: }
364: else if ((tau_1 >= ls->gamma1_i) && (tau_1 < 1.0) &&
365: ((tau_2 < ls->gamma1_i) || (tau_2 >= 1.0))) {
366: tau = tau_1;
367: }
368: else if ((tau_2 >= ls->gamma1_i) && (tau_2 < 1.0) &&
369: ((tau_1 < ls->gamma1_i) || (tau_2 >= 1.0))) {
370: tau = tau_2;
371: }
372: else {
373: tau = tau_max;
374: }
375: }
376: }
377: radius = tau * radius;
378: }
379:
380: if (fmin < f) {
381: f = fmin;
382: info = X->Axpy(sigma, G); CHKERRQ(info);
383: info = TaoComputeGradient(tao, X, G); CHKERRQ(info);
384:
385: info = G->Norm2(&gnorm); CHKERRQ(info);
386: if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
387: SETERRQ(1, "User provided compute function generated Inf or NaN");
388: }
389: needH = 1;
390:
391: info = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); CHKERRQ(info);
392: if (reason != TAO_CONTINUE_ITERATING) {
393: TaoFunctionReturn(0);
394: }
395: }
396: }
397: radius = TaoMax(radius, max_radius);
399: // Modify the radius if it is too large or small
400: radius = TaoMax(radius, ls->min_radius);
401: radius = TaoMin(radius, ls->max_radius);
402: break;
404: default:
405: // Norm of the first direction will initialize radius
406: radius = 0.0;
407: break;
408: }
409: }
411: // Set initial scaling for the BFGS preconditioner
412: // This step is done after computing the initial trust-region radius
413: // since the function value may have decreased
414: if (NLS_PC_BFGS == ls->pc_type) {
415: if (f != 0.0) {
416: delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
417: }
418: else {
419: delta = 2.0 / (gnorm*gnorm);
420: }
421: info = M->SetDelta(delta); CHKERRQ(info);
422: }
424: // Set counter for gradient/reset steps
425: ls->newt = 0;
426: ls->bfgs = 0;
427: ls->sgrad = 0;
428: ls->grad = 0;
430: // Have not converged; continue with Newton method
431: while (reason == TAO_CONTINUE_ITERATING) {
432: ++iter;
434: // Compute the Hessian
435: if (needH) {
436: info = TaoComputeHessian(tao, X, H); CHKERRQ(info);
437: needH = 0;
438: }
440: if ((NLS_PC_BFGS == ls->pc_type) &&
441: (BFGS_SCALE_AHESS == ls->bfgs_scale_type)) {
442: // Obtain diagonal for the bfgs preconditioner
443: info = H->GetDiagonal(Diag); CHKERRQ(info);
444: info = Diag->AbsoluteValue(); CHKERRQ(info);
445: info = Diag->Reciprocal(); CHKERRQ(info);
446: info = M->SetScale(Diag); CHKERRQ(info);
447: }
448:
449: // Shift the Hessian matrix
450: if (pert > 0) {
451: info = H->ShiftDiagonal(pert); CHKERRQ(info);
452: }
453:
454: if (NLS_PC_BFGS == ls->pc_type) {
455: if (BFGS_SCALE_PHESS == ls->bfgs_scale_type) {
456: // Obtain diagonal for the bfgs preconditioner
457: info = H->GetDiagonal(Diag); CHKERRQ(info);
458: info = Diag->AbsoluteValue(); CHKERRQ(info);
459: info = Diag->Reciprocal(); CHKERRQ(info);
460: info = M->SetScale(Diag); CHKERRQ(info);
461: }
463: // Update the limited memory preconditioner
464: info = M->Update(X, G); CHKERRQ(info);
465: ++bfgsUpdates;
466: }
468: // Solve the Newton system of equations
469: info = TaoPreLinearSolve(tao, H); CHKERRQ(info);
470: if (NLS_KSP_NASH == ls->ksp_type ||
471: NLS_KSP_STCG == ls->ksp_type ||
472: NLS_KSP_GLTR == ls->ksp_type) {
473: info = TaoLinearSolveTrustRegion(tao, H, G, D, radius, &success); CHKERRQ(info);
474: info = pls->GetNormDirection(&norm_d); CHKERRQ(info);
475: if (0.0 == radius) {
476: // Radius was uninitialized; use the norm of the direction
477: if (norm_d > 0.0) {
478: radius = norm_d;
480: // Modify the radius if it is too large or small
481: radius = TaoMax(radius, ls->min_radius);
482: radius = TaoMin(radius, ls->max_radius);
483: }
484: else {
485: // The direction was bad; set radius to default value and re-solve
486: // the trust-region subproblem to get a direction
487: info = TaoGetInitialTrustRegionRadius(tao, &radius); CHKERRQ(info);
489: // Modify the radius if it is too large or small
490: radius = TaoMax(radius, ls->min_radius);
491: radius = TaoMin(radius, ls->max_radius);
493: info = TaoLinearSolveTrustRegion(tao, H, G, D, radius, &success); CHKERRQ(info);
494: info = pls->GetNormDirection(&norm_d); CHKERRQ(info);
495: if (norm_d == 0.0) {
496: SETERRQ(1, "Initial direction zero");
497: }
498: }
499: }
500: }
501: else {
502: info = TaoLinearSolve(tao, H, G, D, &success); CHKERRQ(info);
503: }
504: info = D->Negate(); CHKERRQ(info);
506: info = KSPGetConvergedReason(pksp, &ksp_reason); CHKERRQ(info);
507: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
508: (NLS_PC_BFGS == ls->pc_type) && (bfgsUpdates > 1)) {
509: // Preconditioner is numerically indefinite; reset the
510: // approximate if using BFGS preconditioning.
512: if (f != 0.0) {
513: delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
514: }
515: else {
516: delta = 2.0 / (gnorm*gnorm);
517: }
518: info = M->SetDelta(delta); CHKERRQ(info);
519: info = M->Reset(); CHKERRQ(info);
520: info = M->Update(X, G); CHKERRQ(info);
521: bfgsUpdates = 1;
522: }
524: if (KSP_CONVERGED_ATOL == ksp_reason) {
525: ++ls->ksp_atol;
526: }
527: else if (KSP_CONVERGED_RTOL == ksp_reason) {
528: ++ls->ksp_rtol;
529: }
530: else if (KSP_CONVERGED_CG_CONSTRAINED == ksp_reason) {
531: ++ls->ksp_ctol;
532: }
533: else if (KSP_CONVERGED_CG_NEG_CURVE == ksp_reason) {
534: ++ls->ksp_negc;
535: }
536: else if (KSP_DIVERGED_DTOL == ksp_reason) {
537: ++ls->ksp_dtol;
538: }
539: else if (KSP_DIVERGED_ITS == ksp_reason) {
540: ++ls->ksp_iter;
541: }
542: else {
543: ++ls->ksp_othr;
544: }
546: // Check for success (descent direction)
547: info = D->Dot(G, &gdx); CHKERRQ(info);
548: if ((gdx >= 0.0) || TaoInfOrNaN(gdx)) {
549: // Newton step is not descent or direction produced Inf or NaN
550: // Update the perturbation for next time
551: if (pert <= 0.0) {
552: // Initialize the perturbation
553: pert = TaoMin(ls->imax, TaoMax(ls->imin, ls->imfac * gnorm));
554: if (NLS_KSP_GLTR == ls->ksp_type) {
555: info = pls->GetMinEig(&e_min); CHKERRQ(info);
556: pert = TaoMax(pert, -e_min);
557: }
558: }
559: else {
560: // Increase the perturbation
561: pert = TaoMin(ls->pmax, TaoMax(ls->pgfac * pert, ls->pmgfac * gnorm));
562: }
564: if (NLS_PC_BFGS != ls->pc_type) {
565: // We don't have the bfgs matrix around and updated
566: // Must use gradient direction in this case
567: info = D->CopyFrom(G); CHKERRQ(info);
568: info = D->Negate(); CHKERRQ(info);
570: ++ls->grad;
571: stepType = NLS_GRADIENT;
572: }
573: else {
574: // Attempt to use the BFGS direction
575: info = M->Solve(G, D, &success); CHKERRQ(info);
576: info = D->Negate(); CHKERRQ(info);
578: // Check for success (descent direction)
579: info = D->Dot(G, &gdx); CHKERRQ(info);
580: if ((gdx >= 0) || TaoInfOrNaN(gdx)) {
581: // BFGS direction is not descent or direction produced not a number
582: // We can assert bfgsUpdates > 1 in this case because
583: // the first solve produces the scaled gradient direction,
584: // which is guaranteed to be descent
585: //
586: // Use steepest descent direction (scaled)
588: if (f != 0.0) {
589: delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
590: }
591: else {
592: delta = 2.0 / (gnorm*gnorm);
593: }
594: info = M->SetDelta(delta); CHKERRQ(info);
595: info = M->Reset(); CHKERRQ(info);
596: info = M->Update(X, G); CHKERRQ(info);
597: info = M->Solve(G, D, &success); CHKERRQ(info);
598: info = D->Negate(); CHKERRQ(info);
599:
600: bfgsUpdates = 1;
601: ++ls->sgrad;
602: stepType = NLS_SCALED_GRADIENT;
603: }
604: else {
605: if (1 == bfgsUpdates) {
606: // The first BFGS direction is always the scaled gradient
607: ++ls->sgrad;
608: stepType = NLS_SCALED_GRADIENT;
609: }
610: else {
611: ++ls->bfgs;
612: stepType = NLS_BFGS;
613: }
614: }
615: }
616: }
617: else {
618: // Computed Newton step is descent
619: switch (ksp_reason) {
620: case KSP_DIVERGED_NAN:
621: case KSP_DIVERGED_BREAKDOWN:
622: case KSP_DIVERGED_INDEFINITE_MAT:
623: case KSP_DIVERGED_INDEFINITE_PC:
624: case KSP_CONVERGED_CG_NEG_CURVE:
625: // Matrix or preconditioner is indefinite; increase perturbation
626: if (pert <= 0.0) {
627: // Initialize the perturbation
628: pert = TaoMin(ls->imax, TaoMax(ls->imin, ls->imfac * gnorm));
629: if (NLS_KSP_GLTR == ls->ksp_type) {
630: info = pls->GetMinEig(&e_min); CHKERRQ(info);
631: pert = TaoMax(pert, -e_min);
632: }
633: }
634: else {
635: // Increase the perturbation
636: pert = TaoMin(ls->pmax, TaoMax(ls->pgfac * pert, ls->pmgfac * gnorm));
637: }
638: break;
640: default:
641: // Newton step computation is good; decrease perturbation
642: pert = TaoMin(ls->psfac * pert, ls->pmsfac * gnorm);
643: if (pert < ls->pmin) {
644: pert = 0.0;
645: }
646: break;
647: }
649: ++ls->newt;
650: stepType = NLS_NEWTON;
651: }
653: // Perform the linesearch
654: fold = f;
655: info = Xold->CopyFrom(X); CHKERRQ(info);
656: info = Gold->CopyFrom(G); CHKERRQ(info);
658: step = 1.0;
659: info = TaoLineSearchApply(tao, X, G, D, W, &f, &f_full, &step, &status); CHKERRQ(info);
661: while (status && stepType != NLS_GRADIENT) {
662: // Linesearch failed
663: f = fold;
664: info = X->CopyFrom(Xold); CHKERRQ(info);
665: info = G->CopyFrom(Gold); CHKERRQ(info);
667: switch(stepType) {
668: case NLS_NEWTON:
669: // Failed to obtain acceptable iterate with Newton step
670: // Update the perturbation for next time
671: if (pert <= 0.0) {
672: // Initialize the perturbation
673: pert = TaoMin(ls->imax, TaoMax(ls->imin, ls->imfac * gnorm));
674: if (NLS_KSP_GLTR == ls->ksp_type) {
675: info = pls->GetMinEig(&e_min); CHKERRQ(info);
676: pert = TaoMax(pert, -e_min);
677: }
678: }
679: else {
680: // Increase the perturbation
681: pert = TaoMin(ls->pmax, TaoMax(ls->pgfac * pert, ls->pmgfac * gnorm));
682: }
684: if (NLS_PC_BFGS != ls->pc_type) {
685: // We don't have the bfgs matrix around and being updated
686: // Must use gradient direction in this case
687: info = D->CopyFrom(G); CHKERRQ(info);
689: ++ls->grad;
690: stepType = NLS_GRADIENT;
691: }
692: else {
693: // Attempt to use the BFGS direction
694: info = M->Solve(G, D, &success); CHKERRQ(info);
696: // Check for success (descent direction)
697: info = D->Dot(G, &gdx); CHKERRQ(info);
698: if ((gdx <= 0) || TaoInfOrNaN(gdx)) {
699: // BFGS direction is not descent or direction produced not a number
700: // We can assert bfgsUpdates > 1 in this case
701: // Use steepest descent direction (scaled)
702:
703: if (f != 0.0) {
704: delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
705: }
706: else {
707: delta = 2.0 / (gnorm*gnorm);
708: }
709: info = M->SetDelta(delta); CHKERRQ(info);
710: info = M->Reset(); CHKERRQ(info);
711: info = M->Update(X, G); CHKERRQ(info);
712: info = M->Solve(G, D, &success); CHKERRQ(info);
713:
714: bfgsUpdates = 1;
715: ++ls->sgrad;
716: stepType = NLS_SCALED_GRADIENT;
717: }
718: else {
719: if (1 == bfgsUpdates) {
720: // The first BFGS direction is always the scaled gradient
721: ++ls->sgrad;
722: stepType = NLS_SCALED_GRADIENT;
723: }
724: else {
725: ++ls->bfgs;
726: stepType = NLS_BFGS;
727: }
728: }
729: }
730: break;
732: case NLS_BFGS:
733: // Can only enter if pc_type == NLS_PC_BFGS
734: // Failed to obtain acceptable iterate with BFGS step
735: // Attempt to use the scaled gradient direction
737: if (f != 0.0) {
738: delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
739: }
740: else {
741: delta = 2.0 / (gnorm*gnorm);
742: }
743: info = M->SetDelta(delta); CHKERRQ(info);
744: info = M->Reset(); CHKERRQ(info);
745: info = M->Update(X, G); CHKERRQ(info);
746: info = M->Solve(G, D, &success); CHKERRQ(info);
748: bfgsUpdates = 1;
749: ++ls->sgrad;
750: stepType = NLS_SCALED_GRADIENT;
751: break;
753: case NLS_SCALED_GRADIENT:
754: // Can only enter if pc_type == NLS_PC_BFGS
755: // The scaled gradient step did not produce a new iterate;
756: // attemp to use the gradient direction.
757: // Need to make sure we are not using a different diagonal scaling
758: info = M->SetScale(0); CHKERRQ(info);
759: info = M->SetDelta(1.0); CHKERRQ(info);
760: info = M->Reset(); CHKERRQ(info);
761: info = M->Update(X, G); CHKERRQ(info);
762: info = M->Solve(G, D, &success); CHKERRQ(info);
764: bfgsUpdates = 1;
765: ++ls->grad;
766: stepType = NLS_GRADIENT;
767: break;
768: }
769: info = D->Negate(); CHKERRQ(info);
771: // This may be incorrect; linesearch has values for stepmax and stepmin
772: // that should be reset.
773: step = 1.0;
774: info = TaoLineSearchApply(tao, X, G, D, W, &f, &f_full, &step, &status); CHKERRQ(info);
775: }
777: if (status) {
778: // Failed to find an improving point
779: f = fold;
780: info = X->CopyFrom(Xold); CHKERRQ(info);
781: info = G->CopyFrom(Gold); CHKERRQ(info);
782: step = 0.0;
783: }
785: // Update trust region radius
786: if (NLS_KSP_NASH == ls->ksp_type ||
787: NLS_KSP_STCG == ls->ksp_type ||
788: NLS_KSP_GLTR == ls->ksp_type) {
789: switch(ls->update_type) {
790: case NLS_UPDATE_STEP:
791: if (stepType == NLS_NEWTON) {
792: if (step < ls->nu1) {
793: // Very bad step taken; reduce radius
794: radius = ls->omega1 * TaoMin(norm_d, radius);
795: }
796: else if (step < ls->nu2) {
797: // Reasonably bad step taken; reduce radius
798: radius = ls->omega2 * TaoMin(norm_d, radius);
799: }
800: else if (step < ls->nu3) {
801: // Reasonable step was taken; leave radius alone
802: if (ls->omega3 < 1.0) {
803: radius = ls->omega3 * TaoMin(norm_d, radius);
804: }
805: else if (ls->omega3 > 1.0) {
806: radius = TaoMax(ls->omega3 * norm_d, radius);
807: }
808: }
809: else if (step < ls->nu4) {
810: // Full step taken; increase the radius
811: radius = TaoMax(ls->omega4 * norm_d, radius);
812: }
813: else {
814: // More than full step taken; increase the radius
815: radius = TaoMax(ls->omega5 * norm_d, radius);
816: }
817: }
818: else {
819: // Newton step was not good; reduce the radius
820: radius = ls->omega1 * TaoMin(norm_d, radius);
821: }
822: break;
824: case NLS_UPDATE_REDUCTION:
825: if (stepType == NLS_NEWTON) {
826: // Get predicted reduction
827: info = pls->GetObjFcn(&prered); CHKERRQ(info);
829: if (prered >= 0.0) {
830: // The predicted reduction has the wrong sign. This cannot
831: // happen in infinite precision arithmetic. Step should
832: // be rejected!
833: radius = ls->alpha1 * TaoMin(radius, norm_d);
834: }
835: else {
836: if (TaoInfOrNaN(f_full)) {
837: radius = ls->alpha1 * TaoMin(radius, norm_d);
838: }
839: else {
840: // Compute and actual reduction
841: actred = fold - f_full;
842: prered = -prered;
843: if ((fabs(actred) <= ls->epsilon) &&
844: (fabs(prered) <= ls->epsilon)) {
845: kappa = 1.0;
846: }
847: else {
848: kappa = actred / prered;
849: }
850:
851: // Accept of reject the step and update radius
852: if (kappa < ls->eta1) {
853: // Very bad step
854: radius = ls->alpha1 * TaoMin(radius, norm_d);
855: }
856: else if (kappa < ls->eta2) {
857: // Marginal bad step
858: radius = ls->alpha2 * TaoMin(radius, norm_d);
859: }
860: else if (kappa < ls->eta3) {
861: // Reasonable step
862: if (ls->alpha3 < 1.0) {
863: radius = ls->alpha3 * TaoMin(norm_d, radius);
864: }
865: else if (ls->alpha3 > 1.0) {
866: radius = TaoMax(ls->alpha3 * norm_d, radius);
867: }
868: }
869: else if (kappa < ls->eta4) {
870: // Good step
871: radius = TaoMax(ls->alpha4 * norm_d, radius);
872: }
873: else {
874: // Very good step
875: radius = TaoMax(ls->alpha5 * norm_d, radius);
876: }
877: }
878: }
879: }
880: else {
881: // Newton step was not good; reduce the radius
882: radius = ls->alpha1 * TaoMin(norm_d, radius);
883: }
884: break;
886: default:
887: if (stepType == NLS_NEWTON) {
888: // Get predicted reduction
889: info = pls->GetObjFcn(&prered); CHKERRQ(info);
891: if (prered >= 0.0) {
892: // The predicted reduction has the wrong sign. This cannot
893: // happen in infinite precision arithmetic. Step should
894: // be rejected!
895: radius = ls->gamma1 * TaoMin(radius, norm_d);
896: }
897: else {
898: if (TaoInfOrNaN(f_full)) {
899: radius = ls->gamma1 * TaoMin(radius, norm_d);
900: }
901: else {
902: actred = fold - f_full;
903: prered = -prered;
904: if ((fabs(actred) <= ls->epsilon) &&
905: (fabs(prered) <= ls->epsilon)) {
906: kappa = 1.0;
907: }
908: else {
909: kappa = actred / prered;
910: }
912: tau_1 = ls->theta * gdx / (ls->theta * gdx - (1.0 - ls->theta) * prered + actred);
913: tau_2 = ls->theta * gdx / (ls->theta * gdx + (1.0 + ls->theta) * prered - actred);
914: tau_min = TaoMin(tau_1, tau_2);
915: tau_max = TaoMax(tau_1, tau_2);
917: if (kappa >= 1.0 - ls->mu1) {
918: // Great agreement
919: if (tau_max < 1.0) {
920: radius = TaoMax(radius, ls->gamma3 * norm_d);
921: }
922: else if (tau_max > ls->gamma4) {
923: radius = TaoMax(radius, ls->gamma4 * norm_d);
924: }
925: else {
926: radius = TaoMax(radius, tau_max * norm_d);
927: }
928: }
929: else if (kappa >= 1.0 - ls->mu2) {
930: // Good agreement
932: if (tau_max < ls->gamma2) {
933: radius = ls->gamma2 * TaoMin(radius, norm_d);
934: }
935: else if (tau_max > ls->gamma3) {
936: radius = TaoMax(radius, ls->gamma3 * norm_d);
937: }
938: else if (tau_max < 1.0) {
939: radius = tau_max * TaoMin(radius, norm_d);
940: }
941: else {
942: radius = TaoMax(radius, tau_max * norm_d);
943: }
944: }
945: else {
946: // Not good agreement
947: if (tau_min > 1.0) {
948: radius = ls->gamma2 * TaoMin(radius, norm_d);
949: }
950: else if (tau_max < ls->gamma1) {
951: radius = ls->gamma1 * TaoMin(radius, norm_d);
952: }
953: else if ((tau_min < ls->gamma1) && (tau_max >= 1.0)) {
954: radius = ls->gamma1 * TaoMin(radius, norm_d);
955: }
956: else if ((tau_1 >= ls->gamma1) && (tau_1 < 1.0) &&
957: ((tau_2 < ls->gamma1) || (tau_2 >= 1.0))) {
958: radius = tau_1 * TaoMin(radius, norm_d);
959: }
960: else if ((tau_2 >= ls->gamma1) && (tau_2 < 1.0) &&
961: ((tau_1 < ls->gamma1) || (tau_2 >= 1.0))) {
962: radius = tau_2 * TaoMin(radius, norm_d);
963: }
964: else {
965: radius = tau_max * TaoMin(radius, norm_d);
966: }
967: }
968: }
969: }
970: }
971: else {
972: // Newton step was not good; reduce the radius
973: radius = ls->gamma1 * TaoMin(norm_d, radius);
974: }
975: break;
976: }
978: // The radius may have been increased; modify if it is too large
979: radius = TaoMin(radius, ls->max_radius);
980: }
982: // Check for termination
983: info = G->Norm2(&gnorm); CHKERRQ(info);
984: if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
985: SETERRQ(1,"User provided compute function generated Not-a-Number");
986: }
987: needH = 1;
989: info = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason); CHKERRQ(info);
990: }
991: TaoFunctionReturn(0);
992: }
994: /* ---------------------------------------------------------- */
997: static int TaoSetUp_NLS(TAO_SOLVER tao, void *solver)
998: {
999: TAO_NLS *ls = (TAO_NLS *)solver;
1000: TaoVec *X;
1001: TaoMat *H;
1002: int info;
1004: TaoFunctionBegin;
1006: info = TaoGetSolution(tao, &X); CHKERRQ(info);
1007: info = X->Clone(&ls->G); CHKERRQ(info);
1008: info = X->Clone(&ls->D); CHKERRQ(info);
1009: info = X->Clone(&ls->W); CHKERRQ(info);
1011: info = X->Clone(&ls->Xold); CHKERRQ(info);
1012: info = X->Clone(&ls->Gold); CHKERRQ(info);
1014: ls->Diag = 0;
1015: ls->M = 0;
1017: info = TaoSetLagrangianGradientVector(tao, ls->G); CHKERRQ(info);
1018: info = TaoSetStepDirectionVector(tao, ls->D); CHKERRQ(info);
1020: // Set linear solver to default for symmetric matrices
1021: info = TaoGetHessian(tao, &H); CHKERRQ(info);
1022: info = TaoCreateLinearSolver(tao, H, 200, 0); CHKERRQ(info);
1024: // Check sizes for compatability
1025: info = TaoCheckFGH(tao); CHKERRQ(info);
1026: TaoFunctionReturn(0);
1027: }
1029: /*------------------------------------------------------------*/
1032: static int TaoDestroy_NLS(TAO_SOLVER tao, void *solver)
1033: {
1034: TAO_NLS *ls = (TAO_NLS *)solver;
1035: int info;
1037: TaoFunctionBegin;
1038: info = TaoVecDestroy(ls->G); CHKERRQ(info);
1039: info = TaoVecDestroy(ls->D); CHKERRQ(info);
1040: info = TaoVecDestroy(ls->W); CHKERRQ(info);
1042: info = TaoVecDestroy(ls->Xold); CHKERRQ(info);
1043: info = TaoVecDestroy(ls->Gold); CHKERRQ(info);
1045: info = TaoSetLagrangianGradientVector(tao, 0); CHKERRQ(info);
1046: info = TaoSetStepDirectionVector(tao, 0); CHKERRQ(info);
1048: if (ls->Diag) {
1049: info = TaoVecDestroy(ls->Diag); CHKERRQ(info);
1050: ls->Diag = 0;
1051: }
1053: if (ls->M) {
1054: info = TaoMatDestroy(ls->M); CHKERRQ(info);
1055: ls->M = 0;
1056: }
1057: TaoFunctionReturn(0);
1058: }
1060: /*------------------------------------------------------------*/
1063: static int TaoSetOptions_NLS(TAO_SOLVER tao, void *solver)
1064: {
1065: TAO_NLS *ls = (TAO_NLS *)solver;
1066: int info;
1068: TaoFunctionBegin;
1069: info = TaoOptionsHead("Newton line search method for unconstrained optimization"); CHKERRQ(info);
1070: info = TaoOptionList("-tao_nls_ksp_type", "ksp type", "", NLS_KSP, NLS_KSP_TYPES, NLS_KSP[ls->ksp_type], &ls->ksp_type, 0); CHKERRQ(info);
1071: info = TaoOptionList("-tao_nls_pc_type", "pc type", "", NLS_PC, NLS_PC_TYPES, NLS_PC[ls->pc_type], &ls->pc_type, 0); CHKERRQ(info);
1072: info = TaoOptionList("-tao_nls_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[ls->bfgs_scale_type], &ls->bfgs_scale_type, 0); CHKERRQ(info);
1073: info = TaoOptionList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[ls->init_type], &ls->init_type, 0); CHKERRQ(info);
1074: info = TaoOptionList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[ls->update_type], &ls->update_type, 0); CHKERRQ(info);
1075: info = TaoOptionDouble("-tao_nls_sval", "perturbation starting value", "", ls->sval, &ls->sval, 0); CHKERRQ(info);
1076: info = TaoOptionDouble("-tao_nls_imin", "minimum initial perturbation", "", ls->imin, &ls->imin, 0); CHKERRQ(info);
1077: info = TaoOptionDouble("-tao_nls_imax", "maximum initial perturbation", "", ls->imax, &ls->imax, 0); CHKERRQ(info);
1078: info = TaoOptionDouble("-tao_nls_imfac", "initial merit factor", "", ls->imfac, &ls->imfac, 0); CHKERRQ(info);
1079: info = TaoOptionDouble("-tao_nls_pmin", "minimum perturbation", "", ls->pmin, &ls->pmin, 0); CHKERRQ(info);
1080: info = TaoOptionDouble("-tao_nls_pmax", "maximum perturbation", "", ls->pmax, &ls->pmax, 0); CHKERRQ(info);
1081: info = TaoOptionDouble("-tao_nls_pgfac", "growth factor", "", ls->pgfac, &ls->pgfac, 0); CHKERRQ(info);
1082: info = TaoOptionDouble("-tao_nls_psfac", "shrink factor", "", ls->psfac, &ls->psfac, 0); CHKERRQ(info);
1083: info = TaoOptionDouble("-tao_nls_pmgfac", "merit growth factor", "", ls->pmgfac, &ls->pmgfac, 0); CHKERRQ(info);
1084: info = TaoOptionDouble("-tao_nls_pmsfac", "merit shrink factor", "", ls->pmsfac, &ls->pmsfac, 0); CHKERRQ(info);
1085: info = TaoOptionDouble("-tao_nls_eta1", "poor steplength; reduce radius", "", ls->eta1, &ls->eta1, 0); CHKERRQ(info);
1086: info = TaoOptionDouble("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", ls->eta2, &ls->eta2, 0); CHKERRQ(info);
1087: info = TaoOptionDouble("-tao_nls_eta3", "good steplength; increase radius", "", ls->eta3, &ls->eta3, 0); CHKERRQ(info);
1088: info = TaoOptionDouble("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", ls->eta4, &ls->eta4, 0); CHKERRQ(info);
1089: info = TaoOptionDouble("-tao_nls_alpha1", "", "", ls->alpha1, &ls->alpha1, 0); CHKERRQ(info);
1090: info = TaoOptionDouble("-tao_nls_alpha2", "", "", ls->alpha2, &ls->alpha2, 0); CHKERRQ(info);
1091: info = TaoOptionDouble("-tao_nls_alpha3", "", "", ls->alpha3, &ls->alpha3, 0); CHKERRQ(info);
1092: info = TaoOptionDouble("-tao_nls_alpha4", "", "", ls->alpha4, &ls->alpha4, 0); CHKERRQ(info);
1093: info = TaoOptionDouble("-tao_nls_alpha5", "", "", ls->alpha5, &ls->alpha5, 0); CHKERRQ(info);
1094: info = TaoOptionDouble("-tao_nls_nu1", "poor steplength; reduce radius", "", ls->nu1, &ls->nu1, 0); CHKERRQ(info);
1095: info = TaoOptionDouble("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", ls->nu2, &ls->nu2, 0); CHKERRQ(info);
1096: info = TaoOptionDouble("-tao_nls_nu3", "good steplength; increase radius", "", ls->nu3, &ls->nu3, 0); CHKERRQ(info);
1097: info = TaoOptionDouble("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", ls->nu4, &ls->nu4, 0); CHKERRQ(info);
1098: info = TaoOptionDouble("-tao_nls_omega1", "", "", ls->omega1, &ls->omega1, 0); CHKERRQ(info);
1099: info = TaoOptionDouble("-tao_nls_omega2", "", "", ls->omega2, &ls->omega2, 0); CHKERRQ(info);
1100: info = TaoOptionDouble("-tao_nls_omega3", "", "", ls->omega3, &ls->omega3, 0); CHKERRQ(info);
1101: info = TaoOptionDouble("-tao_nls_omega4", "", "", ls->omega4, &ls->omega4, 0); CHKERRQ(info);
1102: info = TaoOptionDouble("-tao_nls_omega5", "", "", ls->omega5, &ls->omega5, 0); CHKERRQ(info);
1103: info = TaoOptionDouble("-tao_nls_mu1_i", "", "", ls->mu1_i, &ls->mu1_i, 0); CHKERRQ(info);
1104: info = TaoOptionDouble("-tao_nls_mu2_i", "", "", ls->mu2_i, &ls->mu2_i, 0); CHKERRQ(info);
1105: info = TaoOptionDouble("-tao_nls_gamma1_i", "", "", ls->gamma1_i, &ls->gamma1_i, 0); CHKERRQ(info);
1106: info = TaoOptionDouble("-tao_nls_gamma2_i", "", "", ls->gamma2_i, &ls->gamma2_i, 0); CHKERRQ(info);
1107: info = TaoOptionDouble("-tao_nls_gamma3_i", "", "", ls->gamma3_i, &ls->gamma3_i, 0); CHKERRQ(info);
1108: info = TaoOptionDouble("-tao_nls_gamma4_i", "", "", ls->gamma4_i, &ls->gamma4_i, 0); CHKERRQ(info);
1109: info = TaoOptionDouble("-tao_nls_theta_i", "", "", ls->theta_i, &ls->theta_i, 0); CHKERRQ(info);
1110: info = TaoOptionDouble("-tao_nls_mu1", "", "", ls->mu1, &ls->mu1, 0); CHKERRQ(info);
1111: info = TaoOptionDouble("-tao_nls_mu2", "", "", ls->mu2, &ls->mu2, 0); CHKERRQ(info);
1112: info = TaoOptionDouble("-tao_nls_gamma1", "", "", ls->gamma1, &ls->gamma1, 0); CHKERRQ(info);
1113: info = TaoOptionDouble("-tao_nls_gamma2", "", "", ls->gamma2, &ls->gamma2, 0); CHKERRQ(info);
1114: info = TaoOptionDouble("-tao_nls_gamma3", "", "", ls->gamma3, &ls->gamma3, 0); CHKERRQ(info);
1115: info = TaoOptionDouble("-tao_nls_gamma4", "", "", ls->gamma4, &ls->gamma4, 0); CHKERRQ(info);
1116: info = TaoOptionDouble("-tao_nls_theta", "", "", ls->theta, &ls->theta, 0); CHKERRQ(info);
1117: info = TaoOptionDouble("-tao_nls_min_radius", "lower bound on initial radius", "", ls->min_radius, &ls->min_radius, 0); CHKERRQ(info);
1118: info = TaoOptionDouble("-tao_nls_max_radius", "upper bound on radius", "", ls->max_radius, &ls->max_radius, 0); CHKERRQ(info);
1119: info = TaoOptionDouble("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", ls->epsilon, &ls->epsilon, 0); CHKERRQ(info);
1120: info = TaoLineSearchSetFromOptions(tao); CHKERRQ(info);
1121: info = TaoOptionsTail(); CHKERRQ(info);
1122: TaoFunctionReturn(0);
1123: }
1126: /*------------------------------------------------------------*/
1129: static int TaoView_NLS(TAO_SOLVER tao,void* solver)
1130: {
1131: TAO_NLS *ls = (TAO_NLS *)solver;
1132: int info;
1134: TaoFunctionBegin;
1135: if (NLS_PC_BFGS == ls->pc_type && ls->M) {
1136: info = TaoPrintInt(tao, " Rejected matrix updates: %d\n", ls->M->GetRejects()); CHKERRQ(info);
1137: }
1138: info = TaoPrintInt(tao, " Newton steps: %d\n", ls->newt); CHKERRQ(info);
1139: info = TaoPrintInt(tao, " BFGS steps: %d\n", ls->bfgs); CHKERRQ(info);
1140: info = TaoPrintInt(tao, " Scaled gradient steps: %d\n", ls->sgrad); CHKERRQ(info);
1141: info = TaoPrintInt(tao, " Gradient steps: %d\n", ls->grad); CHKERRQ(info);
1143: info = TaoPrintInt(tao, " nls ksp atol: %d\n", ls->ksp_atol); CHKERRQ(info);
1144: info = TaoPrintInt(tao, " nls ksp rtol: %d\n", ls->ksp_rtol); CHKERRQ(info);
1145: info = TaoPrintInt(tao, " nls ksp ctol: %d\n", ls->ksp_ctol); CHKERRQ(info);
1146: info = TaoPrintInt(tao, " nls ksp negc: %d\n", ls->ksp_negc); CHKERRQ(info);
1147: info = TaoPrintInt(tao, " nls ksp dtol: %d\n", ls->ksp_dtol); CHKERRQ(info);
1148: info = TaoPrintInt(tao, " nls ksp iter: %d\n", ls->ksp_iter); CHKERRQ(info);
1149: info = TaoPrintInt(tao, " nls ksp othr: %d\n", ls->ksp_othr); CHKERRQ(info);
1150: info = TaoLineSearchView(tao); CHKERRQ(info);
1151: TaoFunctionReturn(0);
1152: }
1154: /* ---------------------------------------------------------- */
1158: int TaoCreate_NLS(TAO_SOLVER tao)
1159: {
1160: TAO_NLS *ls;
1161: int info;
1163: TaoFunctionBegin;
1165: info = TaoNew(TAO_NLS, &ls); CHKERRQ(info);
1166: info = PetscLogObjectMemory(tao, sizeof(TAO_NLS)); CHKERRQ(info);
1168: info = TaoSetTaoSolveRoutine(tao, TaoSolve_NLS, (void *)ls); CHKERRQ(info);
1169: info = TaoSetTaoSetUpDownRoutines(tao, TaoSetUp_NLS, TaoDestroy_NLS); CHKERRQ(info);
1170: info = TaoSetTaoOptionsRoutine(tao, TaoSetOptions_NLS); CHKERRQ(info);
1171: info = TaoSetTaoViewRoutine(tao, TaoView_NLS); CHKERRQ(info);
1173: info = TaoSetMaximumIterates(tao, 50); CHKERRQ(info);
1174: info = TaoSetTolerances(tao, 1e-10, 1e-10, 0, 0); CHKERRQ(info);
1176: info = TaoSetTrustRegionRadius(tao, 100.0); CHKERRQ(info);
1177: info = TaoSetTrustRegionTolerance(tao, 1.0e-12); CHKERRQ(info);
1179: ls->sval = 0.0;
1180: ls->imin = 1.0e-4;
1181: ls->imax = 1.0e+2;
1182: ls->imfac = 1.0e-1;
1184: ls->pmin = 1.0e-12;
1185: ls->pmax = 1.0e+2;
1186: ls->pgfac = 1.0e+1;
1187: ls->psfac = 4.0e-1;
1188: ls->pmgfac = 1.0e-1;
1189: ls->pmsfac = 1.0e-1;
1191: // Default values for trust-region radius update based on steplength
1192: ls->nu1 = 0.25;
1193: ls->nu2 = 0.50;
1194: ls->nu3 = 1.00;
1195: ls->nu4 = 1.25;
1197: ls->omega1 = 0.25;
1198: ls->omega2 = 0.50;
1199: ls->omega3 = 1.00;
1200: ls->omega4 = 2.00;
1201: ls->omega5 = 4.00;
1203: // Default values for trust-region radius update based on reduction
1204: ls->eta1 = 1.0e-4;
1205: ls->eta2 = 0.25;
1206: ls->eta3 = 0.50;
1207: ls->eta4 = 0.90;
1209: ls->alpha1 = 0.25;
1210: ls->alpha2 = 0.50;
1211: ls->alpha3 = 1.00;
1212: ls->alpha4 = 2.00;
1213: ls->alpha5 = 4.00;
1215: // Default values for trust-region radius update based on interpolation
1216: ls->mu1 = 0.10;
1217: ls->mu2 = 0.50;
1219: ls->gamma1 = 0.25;
1220: ls->gamma2 = 0.50;
1221: ls->gamma3 = 2.00;
1222: ls->gamma4 = 4.00;
1224: ls->theta = 0.05;
1226: // Default values for trust region initialization based on interpolation
1227: ls->mu1_i = 0.35;
1228: ls->mu2_i = 0.50;
1230: ls->gamma1_i = 0.0625;
1231: ls->gamma2_i = 0.5;
1232: ls->gamma3_i = 2.0;
1233: ls->gamma4_i = 5.0;
1234:
1235: ls->theta_i = 0.25;
1237: // Remaining parameters
1238: ls->min_radius = 1.0e-10;
1239: ls->max_radius = 1.0e10;
1240: ls->epsilon = 1.0e-6;
1242: ls->ksp_type = NLS_KSP_STCG;
1243: ls->pc_type = NLS_PC_BFGS;
1244: ls->bfgs_scale_type = BFGS_SCALE_PHESS;
1245: ls->init_type = NLS_INIT_INTERPOLATION;
1246: ls->update_type = NLS_UPDATE_STEP;
1248: info = TaoCreateMoreThuenteLineSearch(tao, 0, 0); CHKERRQ(info);
1249: TaoFunctionReturn(0);
1250: }
1253: #endif