Actual source code: nls.c
1: #include "taolinesearch.h"
2: #include "src/matrix/lmvmmat.h"
3: #include "nls.h"
5: #include "petscksp.h"
6: #include "petscpc.h"
7: #include "private/kspimpl.h"
8: #include "private/pcimpl.h"
10: #define NLS_KSP_CG 0
11: #define NLS_KSP_NASH 1
12: #define NLS_KSP_STCG 2
13: #define NLS_KSP_GLTR 3
14: #define NLS_KSP_PETSC 4
15: #define NLS_KSP_TYPES 5
17: #define NLS_PC_NONE 0
18: #define NLS_PC_AHESS 1
19: #define NLS_PC_BFGS 2
20: #define NLS_PC_PETSC 3
21: #define NLS_PC_TYPES 4
23: #define BFGS_SCALE_AHESS 0
24: #define BFGS_SCALE_PHESS 1
25: #define BFGS_SCALE_BFGS 2
26: #define BFGS_SCALE_TYPES 3
28: #define NLS_INIT_CONSTANT 0
29: #define NLS_INIT_DIRECTION 1
30: #define NLS_INIT_INTERPOLATION 2
31: #define NLS_INIT_TYPES 3
33: #define NLS_UPDATE_STEP 0
34: #define NLS_UPDATE_REDUCTION 1
35: #define NLS_UPDATE_INTERPOLATION 2
36: #define NLS_UPDATE_TYPES 3
38: static const char *NLS_KSP[64] = {
39: "cg", "nash", "stcg", "gltr", "petsc"
40: };
42: static const char *NLS_PC[64] = {
43: "none", "ahess", "bfgs", "petsc"
44: };
46: static const char *BFGS_SCALE[64] = {
47: "ahess", "phess", "bfgs"
48: };
50: static const char *NLS_INIT[64] = {
51: "constant", "direction", "interpolation"
52: };
54: static const char *NLS_UPDATE[64] = {
55: "step", "reduction", "interpolation"
56: };
58: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x) ;
59: /* Routine for BFGS preconditioner
62: Implements Newton's Method with a line search approach for solving
63: unconstrained minimization problems. A More'-Thuente line search
64: is used to guarantee that the bfgs preconditioner remains positive
65: definite.
67: The method can shift the Hessian matrix. The shifting procedure is
68: adapted from the PATH algorithm for solving complementarity
69: problems.
71: The linear system solve should be done with a conjugate gradient
72: method, although any method can be used. */
74: #define NLS_NEWTON 0
75: #define NLS_BFGS 1
76: #define NLS_SCALED_GRADIENT 2
77: #define NLS_GRADIENT 3
81: static PetscErrorCode TaoSolve_NLS(TaoSolver tao)
82: {
84: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
86: PC pc;
88: KSPConvergedReason ksp_reason;
89: TaoLineSearchTerminationReason ls_reason;
90: TaoSolverTerminationReason reason;
91:
92: PetscReal fmin, ftrial, f_full, prered, actred, kappa, sigma;
93: PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
94: PetscReal f, fold, gdx, gnorm, pert;
95: PetscReal step = 1.0;
97: PetscReal delta;
98: PetscReal norm_d = 0.0, e_min;
100: MatStructure matflag;
102: PetscInt stepType;
103: PetscInt iter = 0;
104: PetscInt bfgsUpdates = 0;
105: PetscInt n,N,kspits;
106: PetscInt needH;
107:
108: PetscInt i_max = 5;
109: PetscInt j_max = 1;
110: PetscInt i, j;
114: if (tao->XL || tao->XU || tao->ops->computebounds) {
115: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by nls algorithm\n");
116: }
118: /* Initialized variables */
119: pert = nlsP->sval;
121: nlsP->ksp_atol = 0;
122: nlsP->ksp_rtol = 0;
123: nlsP->ksp_dtol = 0;
124: nlsP->ksp_ctol = 0;
125: nlsP->ksp_negc = 0;
126: nlsP->ksp_iter = 0;
127: nlsP->ksp_othr = 0;
131: /* Modify the linear solver to a trust region method if desired */
133: switch(nlsP->ksp_type) {
134: case NLS_KSP_CG:
135: KSPSetType(tao->ksp, KSPCG);
136: if (tao->ksp->ops->setfromoptions) {
137: (*tao->ksp->ops->setfromoptions)(tao->ksp);
138: }
139: break;
141: case NLS_KSP_NASH:
142: KSPSetType(tao->ksp, KSPNASH);
143: if (tao->ksp->ops->setfromoptions) {
144: (*tao->ksp->ops->setfromoptions)(tao->ksp);
145: }
146: break;
148: case NLS_KSP_STCG:
149: KSPSetType(tao->ksp, KSPSTCG);
150: if (tao->ksp->ops->setfromoptions) {
151: (*tao->ksp->ops->setfromoptions)(tao->ksp);
152: }
153: break;
155: case NLS_KSP_GLTR:
156: KSPSetType(tao->ksp, KSPGLTR);
157: if (tao->ksp->ops->setfromoptions) {
158: (*tao->ksp->ops->setfromoptions)(tao->ksp);
159: }
160: break;
162: default:
163: /* Use the method set by the ksp_type */
164: break;
165: }
167: /* Initialize trust-region radius when using nash, stcg, or gltr
168: Will be reset during the first iteration */
169: if (NLS_KSP_NASH == nlsP->ksp_type) {
170: KSPNASHSetRadius(tao->ksp,nlsP->max_radius);
171: } else if (NLS_KSP_STCG == nlsP->ksp_type) {
172: KSPSTCGSetRadius(tao->ksp,nlsP->max_radius);
173: } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
174: KSPGLTRSetRadius(tao->ksp,nlsP->max_radius);
175: }
176:
177:
178: if (NLS_KSP_NASH == nlsP->ksp_type ||
179: NLS_KSP_STCG == nlsP->ksp_type ||
180: NLS_KSP_GLTR == nlsP->ksp_type) {
181: tao->trust = tao->trust0;
183: if (tao->trust < 0.0) {
184: SETERRQ(PETSC_COMM_SELF,1, "Initial radius negative");
185: }
187: /* Modify the radius if it is too large or small */
188: tao->trust = PetscMax(tao->trust, nlsP->min_radius);
189: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
190: }
192: /* Get vectors we will need */
194: if (NLS_PC_BFGS == nlsP->pc_type && !nlsP->M) {
195: VecGetLocalSize(tao->solution,&n);
196: VecGetSize(tao->solution,&N);
197: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&nlsP->M);
198: MatLMVMAllocateVectors(nlsP->M,tao->solution);
199: }
201: /* Check convergence criteria */
202: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
203: VecNorm(tao->gradient,NORM_2,&gnorm);
204: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
205: SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
206: }
207: needH = 1;
209: TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);
210: if (reason != TAO_CONTINUE_ITERATING) {
211: return(0);
212: }
214: /* create vectors for the limited memory preconditioner */
215: if ((NLS_PC_BFGS == nlsP->pc_type) &&
216: (BFGS_SCALE_BFGS != nlsP->bfgs_scale_type)) {
217: if (!nlsP->Diag) {
218: VecDuplicate(tao->solution,&nlsP->Diag);
219: }
220: }
223: /* Modify the preconditioner to use the bfgs approximation */
224: KSPGetPC(tao->ksp, &pc);
225: switch(nlsP->pc_type) {
226: case NLS_PC_NONE:
227: PCSetType(pc, PCNONE);
228: if (pc->ops->setfromoptions) {
229: (*pc->ops->setfromoptions)(pc);
230: }
231: break;
233: case NLS_PC_AHESS:
234: PCSetType(pc, PCJACOBI);
235: if (pc->ops->setfromoptions) {
236: (*pc->ops->setfromoptions)(pc);
237: }
238: PCJacobiSetUseAbs(pc);
239: break;
241: case NLS_PC_BFGS:
242: PCSetType(pc, PCSHELL);
243: if (pc->ops->setfromoptions) {
244: (*pc->ops->setfromoptions)(pc);
245: }
246: PCShellSetName(pc, "bfgs");
247: PCShellSetContext(pc, nlsP->M);
248: PCShellSetApply(pc, MatLMVMSolveShell);
249: break;
251: default:
252: /* Use the pc method set by pc_type */
253: break;
254: }
256: /* Initialize trust-region radius. The initialization is only performed
257: when we are using Nash, Steihaug-Toint or the Generalized Lanczos method. */
258: if (NLS_KSP_NASH == nlsP->ksp_type ||
259: NLS_KSP_STCG == nlsP->ksp_type ||
260: NLS_KSP_GLTR == nlsP->ksp_type) {
261: switch(nlsP->init_type) {
262: case NLS_INIT_CONSTANT:
263: /* Use the initial radius specified */
264: break;
266: case NLS_INIT_INTERPOLATION:
267: /* Use the initial radius specified */
268: max_radius = 0.0;
269:
270: for (j = 0; j < j_max; ++j) {
271: fmin = f;
272: sigma = 0.0;
273:
274: if (needH) {
275: TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag);
276: needH = 0;
277: }
278:
279: for (i = 0; i < i_max; ++i) {
280: VecCopy(tao->solution,nlsP->W);
281: VecAXPY(nlsP->W,-tao->trust/gnorm,tao->gradient);
282: TaoComputeObjective(tao, nlsP->W, &ftrial);
283: if (PetscIsInfOrNanReal(ftrial)) {
284: tau = nlsP->gamma1_i;
285: }
286: else {
287: if (ftrial < fmin) {
288: fmin = ftrial;
289: sigma = -tao->trust / gnorm;
290: }
291:
292: MatMult(tao->hessian, tao->gradient, nlsP->D);
293: VecDot(tao->gradient, nlsP->D, &prered);
294:
295: prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
296: actred = f - ftrial;
297: if ((PetscAbsScalar(actred) <= nlsP->epsilon) &&
298: (PetscAbsScalar(prered) <= nlsP->epsilon)) {
299: kappa = 1.0;
300: }
301: else {
302: kappa = actred / prered;
303: }
304:
305: tau_1 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust + (1.0 - nlsP->theta_i) * prered - actred);
306: tau_2 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust - (1.0 + nlsP->theta_i) * prered + actred);
307: tau_min = PetscMin(tau_1, tau_2);
308: tau_max = PetscMax(tau_1, tau_2);
309:
310: if (PetscAbsScalar(kappa - 1.0) <= nlsP->mu1_i) {
311: /* Great agreement */
312: max_radius = PetscMax(max_radius, tao->trust);
313:
314: if (tau_max < 1.0) {
315: tau = nlsP->gamma3_i;
316: }
317: else if (tau_max > nlsP->gamma4_i) {
318: tau = nlsP->gamma4_i;
319: }
320: else if (tau_1 >= 1.0 && tau_1 <= nlsP->gamma4_i && tau_2 < 1.0) {
321: tau = tau_1;
322: }
323: else if (tau_2 >= 1.0 && tau_2 <= nlsP->gamma4_i && tau_1 < 1.0) {
324: tau = tau_2;
325: }
326: else {
327: tau = tau_max;
328: }
329: }
330: else if (PetscAbsScalar(kappa - 1.0) <= nlsP->mu2_i) {
331: /* Good agreement */
332: max_radius = PetscMax(max_radius, tao->trust);
333:
334: if (tau_max < nlsP->gamma2_i) {
335: tau = nlsP->gamma2_i;
336: }
337: else if (tau_max > nlsP->gamma3_i) {
338: tau = nlsP->gamma3_i;
339: }
340: else {
341: tau = tau_max;
342: }
343: }
344: else {
345: /* Not good agreement */
346: if (tau_min > 1.0) {
347: tau = nlsP->gamma2_i;
348: }
349: else if (tau_max < nlsP->gamma1_i) {
350: tau = nlsP->gamma1_i;
351: }
352: else if ((tau_min < nlsP->gamma1_i) && (tau_max >= 1.0)) {
353: tau = nlsP->gamma1_i;
354: }
355: else if ((tau_1 >= nlsP->gamma1_i) && (tau_1 < 1.0) &&
356: ((tau_2 < nlsP->gamma1_i) || (tau_2 >= 1.0))) {
357: tau = tau_1;
358: }
359: else if ((tau_2 >= nlsP->gamma1_i) && (tau_2 < 1.0) &&
360: ((tau_1 < nlsP->gamma1_i) || (tau_2 >= 1.0))) {
361: tau = tau_2;
362: }
363: else {
364: tau = tau_max;
365: }
366: }
367: }
368: tao->trust = tau * tao->trust;
369: }
370:
371: if (fmin < f) {
372: f = fmin;
373: VecAXPY(tao->solution,sigma,tao->gradient);
374: TaoComputeGradient(tao,tao->solution,tao->gradient);
376: VecNorm(tao->gradient,NORM_2,&gnorm);
377: if (PetscIsInfOrNanReal(gnorm)) {
378: SETERRQ(PETSC_COMM_SELF,1, "User provided compute gradient generated Inf or NaN");
379: }
380: needH = 1;
381:
382: TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);
383: if (reason != TAO_CONTINUE_ITERATING) {
384: return(0);
385: }
386: }
387: }
388: tao->trust = PetscMax(tao->trust, max_radius);
390: /* Modify the radius if it is too large or small */
391: tao->trust = PetscMax(tao->trust, nlsP->min_radius);
392: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
393: break;
395: default:
396: /* Norm of the first direction will initialize radius */
397: tao->trust = 0.0;
398: break;
399: }
400: }
402: /* Set initial scaling for the BFGS preconditioner
403: This step is done after computing the initial trust-region radius
404: since the function value may have decreased */
405: if (NLS_PC_BFGS == nlsP->pc_type) {
406: if (f != 0.0) {
407: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
408: }
409: else {
410: delta = 2.0 / (gnorm*gnorm);
411: }
412: MatLMVMSetDelta(nlsP->M,delta);
413: }
415: /* Set counter for gradient/reset steps*/
416: nlsP->newt = 0;
417: nlsP->bfgs = 0;
418: nlsP->sgrad = 0;
419: nlsP->grad = 0;
421: /* Have not converged; continue with Newton method */
422: while (reason == TAO_CONTINUE_ITERATING) {
423: ++iter;
425: /* Compute the Hessian */
426: if (needH) {
427: TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag);
428: needH = 0;
429: }
431: if ((NLS_PC_BFGS == nlsP->pc_type) &&
432: (BFGS_SCALE_AHESS == nlsP->bfgs_scale_type)) {
433: /* Obtain diagonal for the bfgs preconditioner */
434: MatGetDiagonal(tao->hessian, nlsP->Diag);
435: VecAbs(nlsP->Diag);
436: VecReciprocal(nlsP->Diag);
437: MatLMVMSetScale(nlsP->M,nlsP->Diag);
438: }
439:
440: /* Shift the Hessian matrix */
441: if (pert > 0) {
442: MatShift(tao->hessian, pert);
443: if (tao->hessian != tao->hessian_pre) {
444: MatShift(tao->hessian_pre, pert);
445: }
446: }
448:
449: if (NLS_PC_BFGS == nlsP->pc_type) {
450: if (BFGS_SCALE_PHESS == nlsP->bfgs_scale_type) {
451: /* Obtain diagonal for the bfgs preconditioner */
452: MatGetDiagonal(tao->hessian, nlsP->Diag);
453: VecAbs(nlsP->Diag);
454: VecReciprocal(nlsP->Diag);
455: MatLMVMSetScale(nlsP->M,nlsP->Diag);
456: }
457: /* Update the limited memory preconditioner */
458: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
459: ++bfgsUpdates;
460: }
462: /* Solve the Newton system of equations */
463: KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre,matflag);
464: if (NLS_KSP_NASH == nlsP->ksp_type ||
465: NLS_KSP_STCG == nlsP->ksp_type ||
466: NLS_KSP_GLTR == nlsP->ksp_type) {
468: if (NLS_KSP_NASH == nlsP->ksp_type) {
469: KSPNASHSetRadius(tao->ksp,nlsP->max_radius);
470: } else if (NLS_KSP_STCG == nlsP->ksp_type) {
471: KSPSTCGSetRadius(tao->ksp,nlsP->max_radius);
472: } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
473: KSPGLTRSetRadius(tao->ksp,nlsP->max_radius);
474: }
475:
476: KSPSolve(tao->ksp, tao->gradient, nlsP->D);
477: KSPGetIterationNumber(tao->ksp,&kspits);
478: tao->ksp_its+=kspits;
480: if (NLS_KSP_NASH == nlsP->ksp_type) {
481: KSPNASHGetNormD(tao->ksp,&norm_d);
482: } else if (NLS_KSP_STCG == nlsP->ksp_type) {
483: KSPSTCGGetNormD(tao->ksp,&norm_d);
484: } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
485: KSPGLTRGetNormD(tao->ksp,&norm_d);
486: }
488: if (0.0 == tao->trust) {
489: /* Radius was uninitialized; use the norm of the direction */
490: if (norm_d > 0.0) {
491: tao->trust = norm_d;
493: /* Modify the radius if it is too large or small */
494: tao->trust = PetscMax(tao->trust, nlsP->min_radius);
495: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
496: }
497: else {
498: /* The direction was bad; set radius to default value and re-solve
499: the trust-region subproblem to get a direction */
500: tao->trust = tao->trust0;
502: /* Modify the radius if it is too large or small */
503: tao->trust = PetscMax(tao->trust, nlsP->min_radius);
504: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
506: if (NLS_KSP_NASH == nlsP->ksp_type) {
507: KSPNASHSetRadius(tao->ksp,nlsP->max_radius);
508: } else if (NLS_KSP_STCG == nlsP->ksp_type) {
509: KSPSTCGSetRadius(tao->ksp,nlsP->max_radius);
510: } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
511: KSPGLTRSetRadius(tao->ksp,nlsP->max_radius);
512: }
513:
514: KSPSolve(tao->ksp, tao->gradient, nlsP->D);
515: KSPGetIterationNumber(tao->ksp,&kspits);
516: tao->ksp_its+=kspits;
517: if (NLS_KSP_NASH == nlsP->ksp_type) {
518: KSPNASHGetNormD(tao->ksp,&norm_d);
519: } else if (NLS_KSP_STCG == nlsP->ksp_type) {
520: KSPSTCGGetNormD(tao->ksp,&norm_d);
521: } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
522: KSPGLTRGetNormD(tao->ksp,&norm_d);
523: }
525: if (norm_d == 0.0) {
526: SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
527: }
528: }
529: }
530: }
531: else {
532: KSPSolve(tao->ksp, tao->gradient, nlsP->D);
533: KSPGetIterationNumber(tao->ksp, &kspits);
534: tao->ksp_its += kspits;
535: }
536: VecScale(nlsP->D, -1.0);
537: KSPGetConvergedReason(tao->ksp, &ksp_reason);
538: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
539: (NLS_PC_BFGS == nlsP->pc_type) && (bfgsUpdates > 1)) {
540: /* Preconditioner is numerically indefinite; reset the
541: approximate if using BFGS preconditioning. */
543: if (f != 0.0) {
544: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
545: }
546: else {
547: delta = 2.0 / (gnorm*gnorm);
548: }
549: MatLMVMSetDelta(nlsP->M,delta);
550: MatLMVMReset(nlsP->M);
551: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
552: bfgsUpdates = 1;
553: }
555: if (KSP_CONVERGED_ATOL == ksp_reason) {
556: ++nlsP->ksp_atol;
557: }
558: else if (KSP_CONVERGED_RTOL == ksp_reason) {
559: ++nlsP->ksp_rtol;
560: }
561: else if (KSP_CONVERGED_CG_CONSTRAINED == ksp_reason) {
562: ++nlsP->ksp_ctol;
563: }
564: else if (KSP_CONVERGED_CG_NEG_CURVE == ksp_reason) {
565: ++nlsP->ksp_negc;
566: }
567: else if (KSP_DIVERGED_DTOL == ksp_reason) {
568: ++nlsP->ksp_dtol;
569: }
570: else if (KSP_DIVERGED_ITS == ksp_reason) {
571: ++nlsP->ksp_iter;
572: }
573: else {
574: ++nlsP->ksp_othr;
575: }
577: /* Check for success (descent direction) */
578: VecDot(nlsP->D, tao->gradient, &gdx);
579: if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
580: /* Newton step is not descent or direction produced Inf or NaN
581: Update the perturbation for next time */
582: if (pert <= 0.0) {
583: /* Initialize the perturbation */
584: pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
585: if (NLS_KSP_GLTR == nlsP->ksp_type) {
586: KSPGLTRGetMinEig(tao->ksp,&e_min);
587: pert = PetscMax(pert, -e_min);
588: }
589: }
590: else {
591: /* Increase the perturbation */
592: pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
593: }
595: if (NLS_PC_BFGS != nlsP->pc_type) {
596: /* We don't have the bfgs matrix around and updated
597: Must use gradient direction in this case */
598: VecCopy(tao->gradient, nlsP->D);
599: VecScale(nlsP->D, -1.0);
600: ++nlsP->grad;
601: stepType = NLS_GRADIENT;
602: }
603: else {
604: /* Attempt to use the BFGS direction */
605: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
606: VecScale(nlsP->D, -1.0);
608: /* Check for success (descent direction) */
609: VecDot(tao->gradient, nlsP->D, &gdx);
610: if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
611: /* BFGS direction is not descent or direction produced not a number
612: We can assert bfgsUpdates > 1 in this case because
613: the first solve produces the scaled gradient direction,
614: which is guaranteed to be descent */
616: /* Use steepest descent direction (scaled) */
618: if (f != 0.0) {
619: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
620: }
621: else {
622: delta = 2.0 / (gnorm*gnorm);
623: }
624: MatLMVMSetDelta(nlsP->M, delta);
625: MatLMVMReset(nlsP->M);
626: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
627: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
628: VecScale(nlsP->D, -1.0);
629:
630: bfgsUpdates = 1;
631: ++nlsP->sgrad;
632: stepType = NLS_SCALED_GRADIENT;
633: }
634: else {
635: if (1 == bfgsUpdates) {
636: /* The first BFGS direction is always the scaled gradient */
637: ++nlsP->sgrad;
638: stepType = NLS_SCALED_GRADIENT;
639: }
640: else {
641: ++nlsP->bfgs;
642: stepType = NLS_BFGS;
643: }
644: }
645: }
646: }
647: else {
648: /* Computed Newton step is descent */
649: switch (ksp_reason) {
650: case KSP_DIVERGED_NAN:
651: case KSP_DIVERGED_BREAKDOWN:
652: case KSP_DIVERGED_INDEFINITE_MAT:
653: case KSP_DIVERGED_INDEFINITE_PC:
654: case KSP_CONVERGED_CG_NEG_CURVE:
655: /* Matrix or preconditioner is indefinite; increase perturbation */
656: if (pert <= 0.0) {
657: /* Initialize the perturbation */
658: pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
659: if (NLS_KSP_GLTR == nlsP->ksp_type) {
660: KSPGLTRGetMinEig(tao->ksp, &e_min);
661: pert = PetscMax(pert, -e_min);
662: }
663: }
664: else {
665: /* Increase the perturbation */
666: pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
667: }
668: break;
670: default:
671: /* Newton step computation is good; decrease perturbation */
672: pert = PetscMin(nlsP->psfac * pert, nlsP->pmsfac * gnorm);
673: if (pert < nlsP->pmin) {
674: pert = 0.0;
675: }
676: break;
677: }
679: ++nlsP->newt;
680: stepType = NLS_NEWTON;
681: }
683: /* Perform the linesearch */
684: fold = f;
685: VecCopy(tao->solution, nlsP->Xold);
686: VecCopy(tao->gradient, nlsP->Gold);
688: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason);
689: TaoAddLineSearchCounts(tao);
691: while (ls_reason < 0 && stepType != NLS_GRADIENT) {
692: /* Linesearch failed */
693: f = fold;
694: VecCopy(nlsP->Xold, tao->solution);
695: VecCopy(nlsP->Gold, tao->gradient);
697: switch(stepType) {
698: case NLS_NEWTON:
699: /* Failed to obtain acceptable iterate with Newton 1step
700: Update the perturbation for next time */
701: if (pert <= 0.0) {
702: /* Initialize the perturbation */
703: pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
704: if (NLS_KSP_GLTR == nlsP->ksp_type) {
705: KSPGLTRGetMinEig(tao->ksp,&e_min);
706: pert = PetscMax(pert, -e_min);
707: }
708: }
709: else {
710: /* Increase the perturbation */
711: pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
712: }
714: if (NLS_PC_BFGS != nlsP->pc_type) {
715: /* We don't have the bfgs matrix around and being updated
716: Must use gradient direction in this case */
717: VecCopy(tao->gradient, nlsP->D);
718: ++nlsP->grad;
719: stepType = NLS_GRADIENT;
720: }
721: else {
722: /* Attempt to use the BFGS direction */
723: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
724: /* Check for success (descent direction) */
725: VecDot(tao->solution, nlsP->D, &gdx);
726: if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
727: /* BFGS direction is not descent or direction produced not a number
728: We can assert bfgsUpdates > 1 in this case
729: Use steepest descent direction (scaled) */
730:
731: if (f != 0.0) {
732: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
733: }
734: else {
735: delta = 2.0 / (gnorm*gnorm);
736: }
737: MatLMVMSetDelta(nlsP->M, delta);
738: MatLMVMReset(nlsP->M);
739: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
740: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
741:
742: bfgsUpdates = 1;
743: ++nlsP->sgrad;
744: stepType = NLS_SCALED_GRADIENT;
745: }
746: else {
747: if (1 == bfgsUpdates) {
748: /* The first BFGS direction is always the scaled gradient */
749: ++nlsP->sgrad;
750: stepType = NLS_SCALED_GRADIENT;
751: }
752: else {
753: ++nlsP->bfgs;
754: stepType = NLS_BFGS;
755: }
756: }
757: }
758: break;
760: case NLS_BFGS:
761: /* Can only enter if pc_type == NLS_PC_BFGS
762: Failed to obtain acceptable iterate with BFGS step
763: Attempt to use the scaled gradient direction */
765: if (f != 0.0) {
766: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
767: }
768: else {
769: delta = 2.0 / (gnorm*gnorm);
770: }
771: MatLMVMSetDelta(nlsP->M, delta);
772: MatLMVMReset(nlsP->M);
773: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
774: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
776: bfgsUpdates = 1;
777: ++nlsP->sgrad;
778: stepType = NLS_SCALED_GRADIENT;
779: break;
781: case NLS_SCALED_GRADIENT:
782: /* Can only enter if pc_type == NLS_PC_BFGS
783: The scaled gradient step did not produce a new iterate;
784: attemp to use the gradient direction.
785: Need to make sure we are not using a different diagonal scaling */
786:
787: MatLMVMSetScale(nlsP->M,0);
788: MatLMVMSetDelta(nlsP->M,1.0);
789: MatLMVMReset(nlsP->M);
790: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
791: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
792:
793: bfgsUpdates = 1;
794: ++nlsP->grad;
795: stepType = NLS_GRADIENT;
796: break;
797: }
798: VecScale(nlsP->D, -1.0);
800: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason);
801: TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full);
802: TaoAddLineSearchCounts(tao);
803: }
805: if (ls_reason < 0) {
806: /* Failed to find an improving point */
807: f = fold;
808: VecCopy(nlsP->Xold, tao->solution);
809: VecCopy(nlsP->Gold, tao->gradient);
810: step = 0.0;
811: reason = TAO_DIVERGED_LS_FAILURE;
812: tao->reason = TAO_DIVERGED_LS_FAILURE;
813: break;
814: }
816: /* Update trust region radius */
817: if (NLS_KSP_NASH == nlsP->ksp_type ||
818: NLS_KSP_STCG == nlsP->ksp_type ||
819: NLS_KSP_GLTR == nlsP->ksp_type) {
820: switch(nlsP->update_type) {
821: case NLS_UPDATE_STEP:
822: if (stepType == NLS_NEWTON) {
823: if (step < nlsP->nu1) {
824: /* Very bad step taken; reduce radius */
825: tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
826: }
827: else if (step < nlsP->nu2) {
828: /* Reasonably bad step taken; reduce radius */
829: tao->trust = nlsP->omega2 * PetscMin(norm_d, tao->trust);
830: }
831: else if (step < nlsP->nu3) {
832: /* Reasonable step was taken; leave radius alone */
833: if (nlsP->omega3 < 1.0) {
834: tao->trust = nlsP->omega3 * PetscMin(norm_d, tao->trust);
835: }
836: else if (nlsP->omega3 > 1.0) {
837: tao->trust = PetscMax(nlsP->omega3 * norm_d, tao->trust);
838: }
839: }
840: else if (step < nlsP->nu4) {
841: /* Full step taken; increase the radius */
842: tao->trust = PetscMax(nlsP->omega4 * norm_d, tao->trust);
843: }
844: else {
845: /* More than full step taken; increase the radius */
846: tao->trust = PetscMax(nlsP->omega5 * norm_d, tao->trust);
847: }
848: }
849: else {
850: /* Newton step was not good; reduce the radius */
851: tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
852: }
853: break;
855: case NLS_UPDATE_REDUCTION:
856: if (stepType == NLS_NEWTON) {
857: /* Get predicted reduction */
859: if (NLS_KSP_STCG == nlsP->ksp_type) {
860: KSPSTCGGetObjFcn(tao->ksp,&prered);
861: } else if (NLS_KSP_NASH == nlsP->ksp_type) {
862: KSPNASHGetObjFcn(tao->ksp,&prered);
863: } else {
864: KSPGLTRGetObjFcn(tao->ksp,&prered);
865: }
866:
867:
868:
870: if (prered >= 0.0) {
871: /* The predicted reduction has the wrong sign. This cannot */
872: /* happen in infinite precision arithmetic. Step should */
873: /* be rejected! */
874: tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
875: }
876: else {
877: if (PetscIsInfOrNanReal(f_full)) {
878: tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
879: }
880: else {
881: /* Compute and actual reduction */
882: actred = fold - f_full;
883: prered = -prered;
884: if ((PetscAbsScalar(actred) <= nlsP->epsilon) &&
885: (PetscAbsScalar(prered) <= nlsP->epsilon)) {
886: kappa = 1.0;
887: }
888: else {
889: kappa = actred / prered;
890: }
891:
892: /* Accept of reject the step and update radius */
893: if (kappa < nlsP->eta1) {
894: /* Very bad step */
895: tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
896: }
897: else if (kappa < nlsP->eta2) {
898: /* Marginal bad step */
899: tao->trust = nlsP->alpha2 * PetscMin(tao->trust, norm_d);
900: }
901: else if (kappa < nlsP->eta3) {
902: /* Reasonable step */
903: if (nlsP->alpha3 < 1.0) {
904: tao->trust = nlsP->alpha3 * PetscMin(norm_d, tao->trust);
905: }
906: else if (nlsP->alpha3 > 1.0) {
907: tao->trust = PetscMax(nlsP->alpha3 * norm_d, tao->trust);
908: }
909: }
910: else if (kappa < nlsP->eta4) {
911: /* Good step */
912: tao->trust = PetscMax(nlsP->alpha4 * norm_d, tao->trust);
913: }
914: else {
915: /* Very good step */
916: tao->trust = PetscMax(nlsP->alpha5 * norm_d, tao->trust);
917: }
918: }
919: }
920: }
921: else {
922: /* Newton step was not good; reduce the radius */
923: tao->trust = nlsP->alpha1 * PetscMin(norm_d, tao->trust);
924: }
925: break;
927: default:
928: if (stepType == NLS_NEWTON) {
930: if (NLS_KSP_STCG == nlsP->ksp_type) {
931: KSPSTCGGetObjFcn(tao->ksp,&prered);
932: } else if (NLS_KSP_NASH == nlsP->ksp_type) {
933: KSPNASHGetObjFcn(tao->ksp,&prered);
934: } else {
935: KSPGLTRGetObjFcn(tao->ksp,&prered);
936: }
937: if (prered >= 0.0) {
938: /* The predicted reduction has the wrong sign. This cannot */
939: /* happen in infinite precision arithmetic. Step should */
940: /* be rejected! */
941: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
942: }
943: else {
944: if (PetscIsInfOrNanReal(f_full)) {
945: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
946: }
947: else {
948: actred = fold - f_full;
949: prered = -prered;
950: if ((PetscAbsScalar(actred) <= nlsP->epsilon) &&
951: (PetscAbsScalar(prered) <= nlsP->epsilon)) {
952: kappa = 1.0;
953: }
954: else {
955: kappa = actred / prered;
956: }
958: tau_1 = nlsP->theta * gdx / (nlsP->theta * gdx - (1.0 - nlsP->theta) * prered + actred);
959: tau_2 = nlsP->theta * gdx / (nlsP->theta * gdx + (1.0 + nlsP->theta) * prered - actred);
960: tau_min = PetscMin(tau_1, tau_2);
961: tau_max = PetscMax(tau_1, tau_2);
963: if (kappa >= 1.0 - nlsP->mu1) {
964: /* Great agreement */
965: if (tau_max < 1.0) {
966: tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
967: }
968: else if (tau_max > nlsP->gamma4) {
969: tao->trust = PetscMax(tao->trust, nlsP->gamma4 * norm_d);
970: }
971: else {
972: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
973: }
974: }
975: else if (kappa >= 1.0 - nlsP->mu2) {
976: /* Good agreement */
978: if (tau_max < nlsP->gamma2) {
979: tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
980: }
981: else if (tau_max > nlsP->gamma3) {
982: tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
983: }
984: else if (tau_max < 1.0) {
985: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
986: }
987: else {
988: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
989: }
990: }
991: else {
992: /* Not good agreement */
993: if (tau_min > 1.0) {
994: tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
995: }
996: else if (tau_max < nlsP->gamma1) {
997: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
998: }
999: else if ((tau_min < nlsP->gamma1) && (tau_max >= 1.0)) {
1000: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
1001: }
1002: else if ((tau_1 >= nlsP->gamma1) && (tau_1 < 1.0) &&
1003: ((tau_2 < nlsP->gamma1) || (tau_2 >= 1.0))) {
1004: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
1005: }
1006: else if ((tau_2 >= nlsP->gamma1) && (tau_2 < 1.0) &&
1007: ((tau_1 < nlsP->gamma1) || (tau_2 >= 1.0))) {
1008: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
1009: }
1010: else {
1011: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
1012: }
1013: }
1014: }
1015: }
1016: }
1017: else {
1018: /* Newton step was not good; reduce the radius */
1019: tao->trust = nlsP->gamma1 * PetscMin(norm_d, tao->trust);
1020: }
1021: break;
1022: }
1024: /* The radius may have been increased; modify if it is too large */
1025: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
1026: }
1028: /* Check for termination */
1029: VecNorm(tao->gradient, NORM_2, &gnorm);
1030: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
1031: SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number");
1032: }
1033: needH = 1;
1035: TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);
1036: }
1037: return(0);
1038: }
1040: /* ---------------------------------------------------------- */
1043: static PetscErrorCode TaoSetUp_NLS(TaoSolver tao)
1044: {
1045: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
1050: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
1051: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
1052: if (!nlsP->W) {VecDuplicate(tao->solution,&nlsP->W); }
1053: if (!nlsP->D) {VecDuplicate(tao->solution,&nlsP->D); }
1054: if (!nlsP->Xold) {VecDuplicate(tao->solution,&nlsP->Xold); }
1055: if (!nlsP->Gold) {VecDuplicate(tao->solution,&nlsP->Gold); }
1057: nlsP->Diag = 0;
1058: nlsP->M = 0;
1060: return(0);
1061: }
1063: /*------------------------------------------------------------*/
1066: static PetscErrorCode TaoDestroy_NLS(TaoSolver tao)
1067: {
1068: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
1072: if (tao->setupcalled) {
1073: VecDestroy(&nlsP->D);
1074: VecDestroy(&nlsP->W);
1075: VecDestroy(&nlsP->Xold);
1076: VecDestroy(&nlsP->Gold);
1077: }
1078: if (nlsP->Diag) {
1079: VecDestroy(&nlsP->Diag);
1080: }
1081: if (nlsP->M) {
1082: MatDestroy(&nlsP->M);
1083: }
1085: PetscFree(tao->data);
1086: tao->data = PETSC_NULL;
1088: return(0);
1089: }
1091: /*------------------------------------------------------------*/
1094: static PetscErrorCode TaoSetFromOptions_NLS(TaoSolver tao)
1095: {
1096: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
1100: PetscOptionsHead("Newton line search method for unconstrained optimization");
1101: PetscOptionsEList("-tao_nls_ksp_type", "ksp type", "", NLS_KSP, NLS_KSP_TYPES, NLS_KSP[nlsP->ksp_type], &nlsP->ksp_type, 0);
1102: PetscOptionsEList("-tao_nls_pc_type", "pc type", "", NLS_PC, NLS_PC_TYPES, NLS_PC[nlsP->pc_type], &nlsP->pc_type, 0);
1103: PetscOptionsEList("-tao_nls_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[nlsP->bfgs_scale_type], &nlsP->bfgs_scale_type, 0);
1104: PetscOptionsEList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[nlsP->init_type], &nlsP->init_type, 0);
1105: PetscOptionsEList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[nlsP->update_type], &nlsP->update_type, 0);
1106: PetscOptionsReal("-tao_nls_sval", "perturbation starting value", "", nlsP->sval, &nlsP->sval, 0);
1107: PetscOptionsReal("-tao_nls_imin", "minimum initial perturbation", "", nlsP->imin, &nlsP->imin, 0);
1108: PetscOptionsReal("-tao_nls_imax", "maximum initial perturbation", "", nlsP->imax, &nlsP->imax, 0);
1109: PetscOptionsReal("-tao_nls_imfac", "initial merit factor", "", nlsP->imfac, &nlsP->imfac, 0);
1110: PetscOptionsReal("-tao_nls_pmin", "minimum perturbation", "", nlsP->pmin, &nlsP->pmin, 0);
1111: PetscOptionsReal("-tao_nls_pmax", "maximum perturbation", "", nlsP->pmax, &nlsP->pmax, 0);
1112: PetscOptionsReal("-tao_nls_pgfac", "growth factor", "", nlsP->pgfac, &nlsP->pgfac, 0);
1113: PetscOptionsReal("-tao_nls_psfac", "shrink factor", "", nlsP->psfac, &nlsP->psfac, 0);
1114: PetscOptionsReal("-tao_nls_pmgfac", "merit growth factor", "", nlsP->pmgfac, &nlsP->pmgfac, 0);
1115: PetscOptionsReal("-tao_nls_pmsfac", "merit shrink factor", "", nlsP->pmsfac, &nlsP->pmsfac, 0);
1116: PetscOptionsReal("-tao_nls_eta1", "poor steplength; reduce radius", "", nlsP->eta1, &nlsP->eta1, 0);
1117: PetscOptionsReal("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", nlsP->eta2, &nlsP->eta2, 0);
1118: PetscOptionsReal("-tao_nls_eta3", "good steplength; increase radius", "", nlsP->eta3, &nlsP->eta3, 0);
1119: PetscOptionsReal("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", nlsP->eta4, &nlsP->eta4, 0);
1120: PetscOptionsReal("-tao_nls_alpha1", "", "", nlsP->alpha1, &nlsP->alpha1, 0);
1121: PetscOptionsReal("-tao_nls_alpha2", "", "", nlsP->alpha2, &nlsP->alpha2, 0);
1122: PetscOptionsReal("-tao_nls_alpha3", "", "", nlsP->alpha3, &nlsP->alpha3, 0);
1123: PetscOptionsReal("-tao_nls_alpha4", "", "", nlsP->alpha4, &nlsP->alpha4, 0);
1124: PetscOptionsReal("-tao_nls_alpha5", "", "", nlsP->alpha5, &nlsP->alpha5, 0);
1125: PetscOptionsReal("-tao_nls_nu1", "poor steplength; reduce radius", "", nlsP->nu1, &nlsP->nu1, 0);
1126: PetscOptionsReal("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", nlsP->nu2, &nlsP->nu2, 0);
1127: PetscOptionsReal("-tao_nls_nu3", "good steplength; increase radius", "", nlsP->nu3, &nlsP->nu3, 0);
1128: PetscOptionsReal("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", nlsP->nu4, &nlsP->nu4, 0);
1129: PetscOptionsReal("-tao_nls_omega1", "", "", nlsP->omega1, &nlsP->omega1, 0);
1130: PetscOptionsReal("-tao_nls_omega2", "", "", nlsP->omega2, &nlsP->omega2, 0);
1131: PetscOptionsReal("-tao_nls_omega3", "", "", nlsP->omega3, &nlsP->omega3, 0);
1132: PetscOptionsReal("-tao_nls_omega4", "", "", nlsP->omega4, &nlsP->omega4, 0);
1133: PetscOptionsReal("-tao_nls_omega5", "", "", nlsP->omega5, &nlsP->omega5, 0);
1134: PetscOptionsReal("-tao_nls_mu1_i", "", "", nlsP->mu1_i, &nlsP->mu1_i, 0);
1135: PetscOptionsReal("-tao_nls_mu2_i", "", "", nlsP->mu2_i, &nlsP->mu2_i, 0);
1136: PetscOptionsReal("-tao_nls_gamma1_i", "", "", nlsP->gamma1_i, &nlsP->gamma1_i, 0);
1137: PetscOptionsReal("-tao_nls_gamma2_i", "", "", nlsP->gamma2_i, &nlsP->gamma2_i, 0);
1138: PetscOptionsReal("-tao_nls_gamma3_i", "", "", nlsP->gamma3_i, &nlsP->gamma3_i, 0);
1139: PetscOptionsReal("-tao_nls_gamma4_i", "", "", nlsP->gamma4_i, &nlsP->gamma4_i, 0);
1140: PetscOptionsReal("-tao_nls_theta_i", "", "", nlsP->theta_i, &nlsP->theta_i, 0);
1141: PetscOptionsReal("-tao_nls_mu1", "", "", nlsP->mu1, &nlsP->mu1, 0);
1142: PetscOptionsReal("-tao_nls_mu2", "", "", nlsP->mu2, &nlsP->mu2, 0);
1143: PetscOptionsReal("-tao_nls_gamma1", "", "", nlsP->gamma1, &nlsP->gamma1, 0);
1144: PetscOptionsReal("-tao_nls_gamma2", "", "", nlsP->gamma2, &nlsP->gamma2, 0);
1145: PetscOptionsReal("-tao_nls_gamma3", "", "", nlsP->gamma3, &nlsP->gamma3, 0);
1146: PetscOptionsReal("-tao_nls_gamma4", "", "", nlsP->gamma4, &nlsP->gamma4, 0);
1147: PetscOptionsReal("-tao_nls_theta", "", "", nlsP->theta, &nlsP->theta, 0);
1148: PetscOptionsReal("-tao_nls_min_radius", "lower bound on initial radius", "", nlsP->min_radius, &nlsP->min_radius, 0);
1149: PetscOptionsReal("-tao_nls_max_radius", "upper bound on radius", "", nlsP->max_radius, &nlsP->max_radius, 0);
1150: PetscOptionsReal("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", nlsP->epsilon, &nlsP->epsilon, 0);
1151: PetscOptionsTail();
1152: TaoLineSearchSetFromOptions(tao->linesearch);
1153: KSPSetFromOptions(tao->ksp);
1154: return(0);
1155: }
1158: /*------------------------------------------------------------*/
1161: static PetscErrorCode TaoView_NLS(TaoSolver tao, PetscViewer viewer)
1162: {
1163: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
1164: PetscInt nrejects;
1165: PetscBool isascii;
1170: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1171: if (isascii) {
1172: PetscViewerASCIIPushTab(viewer);
1173: if (NLS_PC_BFGS == nlsP->pc_type && nlsP->M) {
1174: MatLMVMGetRejects(nlsP->M,&nrejects);
1175: PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n",nrejects);
1176: }
1177: PetscViewerASCIIPrintf(viewer, "Newton steps: %D\n", nlsP->newt);
1178: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", nlsP->bfgs);
1179: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", nlsP->sgrad);
1180: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", nlsP->grad);
1182: PetscViewerASCIIPrintf(viewer, "nls ksp atol: %D\n", nlsP->ksp_atol);
1183: PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %D\n", nlsP->ksp_rtol);
1184: PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %D\n", nlsP->ksp_ctol);
1185: PetscViewerASCIIPrintf(viewer, "nls ksp negc: %D\n", nlsP->ksp_negc);
1186: PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %D\n", nlsP->ksp_dtol);
1187: PetscViewerASCIIPrintf(viewer, "nls ksp iter: %D\n", nlsP->ksp_iter);
1188: PetscViewerASCIIPrintf(viewer, "nls ksp othr: %D\n", nlsP->ksp_othr);
1189: PetscViewerASCIIPopTab(viewer);
1190: } else {
1191: SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO NLS",((PetscObject)viewer)->type_name);
1192: }
1193: return(0);
1194: }
1196: /* ---------------------------------------------------------- */
1200: PetscErrorCode TaoCreate_NLS(TaoSolver tao)
1201: {
1202: TAO_NLS *nlsP;
1203: const char *morethuente_type = TAOLINESEARCH_MT;
1207: PetscNewLog(tao,TAO_NLS,&nlsP);
1209: tao->ops->setup = TaoSetUp_NLS;
1210: tao->ops->solve = TaoSolve_NLS;
1211: tao->ops->view = TaoView_NLS;
1212: tao->ops->setfromoptions = TaoSetFromOptions_NLS;
1213: tao->ops->destroy = TaoDestroy_NLS;
1215: tao->max_it = 50;
1216: tao->fatol = 1e-10;
1217: tao->frtol = 1e-10;
1218: tao->data = (void*)nlsP;
1219: tao->trust0 = 100.0;
1221: nlsP->sval = 0.0;
1222: nlsP->imin = 1.0e-4;
1223: nlsP->imax = 1.0e+2;
1224: nlsP->imfac = 1.0e-1;
1226: nlsP->pmin = 1.0e-12;
1227: nlsP->pmax = 1.0e+2;
1228: nlsP->pgfac = 1.0e+1;
1229: nlsP->psfac = 4.0e-1;
1230: nlsP->pmgfac = 1.0e-1;
1231: nlsP->pmsfac = 1.0e-1;
1233: /* Default values for trust-region radius update based on steplength */
1234: nlsP->nu1 = 0.25;
1235: nlsP->nu2 = 0.50;
1236: nlsP->nu3 = 1.00;
1237: nlsP->nu4 = 1.25;
1239: nlsP->omega1 = 0.25;
1240: nlsP->omega2 = 0.50;
1241: nlsP->omega3 = 1.00;
1242: nlsP->omega4 = 2.00;
1243: nlsP->omega5 = 4.00;
1245: /* Default values for trust-region radius update based on reduction */
1246: nlsP->eta1 = 1.0e-4;
1247: nlsP->eta2 = 0.25;
1248: nlsP->eta3 = 0.50;
1249: nlsP->eta4 = 0.90;
1251: nlsP->alpha1 = 0.25;
1252: nlsP->alpha2 = 0.50;
1253: nlsP->alpha3 = 1.00;
1254: nlsP->alpha4 = 2.00;
1255: nlsP->alpha5 = 4.00;
1257: /* Default values for trust-region radius update based on interpolation */
1258: nlsP->mu1 = 0.10;
1259: nlsP->mu2 = 0.50;
1261: nlsP->gamma1 = 0.25;
1262: nlsP->gamma2 = 0.50;
1263: nlsP->gamma3 = 2.00;
1264: nlsP->gamma4 = 4.00;
1266: nlsP->theta = 0.05;
1268: /* Default values for trust region initialization based on interpolation */
1269: nlsP->mu1_i = 0.35;
1270: nlsP->mu2_i = 0.50;
1272: nlsP->gamma1_i = 0.0625;
1273: nlsP->gamma2_i = 0.5;
1274: nlsP->gamma3_i = 2.0;
1275: nlsP->gamma4_i = 5.0;
1276:
1277: nlsP->theta_i = 0.25;
1279: /* Remaining parameters */
1280: nlsP->min_radius = 1.0e-10;
1281: nlsP->max_radius = 1.0e10;
1282: nlsP->epsilon = 1.0e-6;
1284: nlsP->ksp_type = NLS_KSP_STCG;
1285: nlsP->pc_type = NLS_PC_BFGS;
1286: nlsP->bfgs_scale_type = BFGS_SCALE_PHESS;
1287: nlsP->init_type = NLS_INIT_INTERPOLATION;
1288: nlsP->update_type = NLS_UPDATE_STEP;
1290: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
1291: TaoLineSearchSetType(tao->linesearch,morethuente_type);
1292: TaoLineSearchUseTaoSolverRoutines(tao->linesearch,tao);
1294: /* Set linear solver to default for symmetric matrices */
1295: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
1297: return(0);
1298: }
1307: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x)
1308: {
1310: Mat M;
1315: PCShellGetContext(pc,(void**)&M);
1316: MatLMVMSolve(M, b, x);
1317: return(0);
1318: }