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 != TAOLINESEARCH_SUCCESS &&
692: ls_reason != TAOLINESEARCH_SUCCESS_USER &&
693: stepType != NLS_GRADIENT) { /* Linesearch failed */
694: f = fold;
695: VecCopy(nlsP->Xold, tao->solution);
696: VecCopy(nlsP->Gold, tao->gradient);
698: switch(stepType) {
699: case NLS_NEWTON:
700: /* Failed to obtain acceptable iterate with Newton 1step
701: Update the perturbation for next time */
702: if (pert <= 0.0) {
703: /* Initialize the perturbation */
704: pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
705: if (NLS_KSP_GLTR == nlsP->ksp_type) {
706: KSPGLTRGetMinEig(tao->ksp,&e_min);
707: pert = PetscMax(pert, -e_min);
708: }
709: }
710: else {
711: /* Increase the perturbation */
712: pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
713: }
715: if (NLS_PC_BFGS != nlsP->pc_type) {
716: /* We don't have the bfgs matrix around and being updated
717: Must use gradient direction in this case */
718: VecCopy(tao->gradient, nlsP->D);
719: ++nlsP->grad;
720: stepType = NLS_GRADIENT;
721: }
722: else {
723: /* Attempt to use the BFGS direction */
724: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
725: /* Check for success (descent direction) */
726: VecDot(tao->solution, nlsP->D, &gdx);
727: if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
728: /* BFGS direction is not descent or direction produced not a number
729: We can assert bfgsUpdates > 1 in this case
730: Use steepest descent direction (scaled) */
731:
732: if (f != 0.0) {
733: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
734: }
735: else {
736: delta = 2.0 / (gnorm*gnorm);
737: }
738: MatLMVMSetDelta(nlsP->M, delta);
739: MatLMVMReset(nlsP->M);
740: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
741: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
742:
743: bfgsUpdates = 1;
744: ++nlsP->sgrad;
745: stepType = NLS_SCALED_GRADIENT;
746: }
747: else {
748: if (1 == bfgsUpdates) {
749: /* The first BFGS direction is always the scaled gradient */
750: ++nlsP->sgrad;
751: stepType = NLS_SCALED_GRADIENT;
752: }
753: else {
754: ++nlsP->bfgs;
755: stepType = NLS_BFGS;
756: }
757: }
758: }
759: break;
761: case NLS_BFGS:
762: /* Can only enter if pc_type == NLS_PC_BFGS
763: Failed to obtain acceptable iterate with BFGS step
764: Attempt to use the scaled gradient direction */
766: if (f != 0.0) {
767: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
768: }
769: else {
770: delta = 2.0 / (gnorm*gnorm);
771: }
772: MatLMVMSetDelta(nlsP->M, delta);
773: MatLMVMReset(nlsP->M);
774: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
775: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
777: bfgsUpdates = 1;
778: ++nlsP->sgrad;
779: stepType = NLS_SCALED_GRADIENT;
780: break;
782: case NLS_SCALED_GRADIENT:
783: /* Can only enter if pc_type == NLS_PC_BFGS
784: The scaled gradient step did not produce a new iterate;
785: attemp to use the gradient direction.
786: Need to make sure we are not using a different diagonal scaling */
787:
788: MatLMVMSetScale(nlsP->M,0);
789: MatLMVMSetDelta(nlsP->M,1.0);
790: MatLMVMReset(nlsP->M);
791: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
792: MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D);
793:
794: bfgsUpdates = 1;
795: ++nlsP->grad;
796: stepType = NLS_GRADIENT;
797: break;
798: }
799: VecScale(nlsP->D, -1.0);
801: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason);
802: TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full);
803: TaoAddLineSearchCounts(tao);
804: }
806: if (ls_reason != TAOLINESEARCH_SUCCESS &&
807: ls_reason != TAOLINESEARCH_SUCCESS_USER) {
808: /* Failed to find an improving point */
809: f = fold;
810: VecCopy(nlsP->Xold, tao->solution);
811: VecCopy(nlsP->Gold, tao->gradient);
812: step = 0.0;
813: reason = TAO_DIVERGED_LS_FAILURE;
814: tao->reason = TAO_DIVERGED_LS_FAILURE;
815: break;
816: }
818: /* Update trust region radius */
819: if (NLS_KSP_NASH == nlsP->ksp_type ||
820: NLS_KSP_STCG == nlsP->ksp_type ||
821: NLS_KSP_GLTR == nlsP->ksp_type) {
822: switch(nlsP->update_type) {
823: case NLS_UPDATE_STEP:
824: if (stepType == NLS_NEWTON) {
825: if (step < nlsP->nu1) {
826: /* Very bad step taken; reduce radius */
827: tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
828: }
829: else if (step < nlsP->nu2) {
830: /* Reasonably bad step taken; reduce radius */
831: tao->trust = nlsP->omega2 * PetscMin(norm_d, tao->trust);
832: }
833: else if (step < nlsP->nu3) {
834: /* Reasonable step was taken; leave radius alone */
835: if (nlsP->omega3 < 1.0) {
836: tao->trust = nlsP->omega3 * PetscMin(norm_d, tao->trust);
837: }
838: else if (nlsP->omega3 > 1.0) {
839: tao->trust = PetscMax(nlsP->omega3 * norm_d, tao->trust);
840: }
841: }
842: else if (step < nlsP->nu4) {
843: /* Full step taken; increase the radius */
844: tao->trust = PetscMax(nlsP->omega4 * norm_d, tao->trust);
845: }
846: else {
847: /* More than full step taken; increase the radius */
848: tao->trust = PetscMax(nlsP->omega5 * norm_d, tao->trust);
849: }
850: }
851: else {
852: /* Newton step was not good; reduce the radius */
853: tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
854: }
855: break;
857: case NLS_UPDATE_REDUCTION:
858: if (stepType == NLS_NEWTON) {
859: /* Get predicted reduction */
861: if (NLS_KSP_STCG == nlsP->ksp_type) {
862: KSPSTCGGetObjFcn(tao->ksp,&prered);
863: } else if (NLS_KSP_NASH == nlsP->ksp_type) {
864: KSPNASHGetObjFcn(tao->ksp,&prered);
865: } else {
866: KSPGLTRGetObjFcn(tao->ksp,&prered);
867: }
868:
869:
870:
872: if (prered >= 0.0) {
873: /* The predicted reduction has the wrong sign. This cannot */
874: /* happen in infinite precision arithmetic. Step should */
875: /* be rejected! */
876: tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
877: }
878: else {
879: if (PetscIsInfOrNanReal(f_full)) {
880: tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
881: }
882: else {
883: /* Compute and actual reduction */
884: actred = fold - f_full;
885: prered = -prered;
886: if ((PetscAbsScalar(actred) <= nlsP->epsilon) &&
887: (PetscAbsScalar(prered) <= nlsP->epsilon)) {
888: kappa = 1.0;
889: }
890: else {
891: kappa = actred / prered;
892: }
893:
894: /* Accept of reject the step and update radius */
895: if (kappa < nlsP->eta1) {
896: /* Very bad step */
897: tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
898: }
899: else if (kappa < nlsP->eta2) {
900: /* Marginal bad step */
901: tao->trust = nlsP->alpha2 * PetscMin(tao->trust, norm_d);
902: }
903: else if (kappa < nlsP->eta3) {
904: /* Reasonable step */
905: if (nlsP->alpha3 < 1.0) {
906: tao->trust = nlsP->alpha3 * PetscMin(norm_d, tao->trust);
907: }
908: else if (nlsP->alpha3 > 1.0) {
909: tao->trust = PetscMax(nlsP->alpha3 * norm_d, tao->trust);
910: }
911: }
912: else if (kappa < nlsP->eta4) {
913: /* Good step */
914: tao->trust = PetscMax(nlsP->alpha4 * norm_d, tao->trust);
915: }
916: else {
917: /* Very good step */
918: tao->trust = PetscMax(nlsP->alpha5 * norm_d, tao->trust);
919: }
920: }
921: }
922: }
923: else {
924: /* Newton step was not good; reduce the radius */
925: tao->trust = nlsP->alpha1 * PetscMin(norm_d, tao->trust);
926: }
927: break;
929: default:
930: if (stepType == NLS_NEWTON) {
932: if (NLS_KSP_STCG == nlsP->ksp_type) {
933: KSPSTCGGetObjFcn(tao->ksp,&prered);
934: } else if (NLS_KSP_NASH == nlsP->ksp_type) {
935: KSPNASHGetObjFcn(tao->ksp,&prered);
936: } else {
937: KSPGLTRGetObjFcn(tao->ksp,&prered);
938: }
939: if (prered >= 0.0) {
940: /* The predicted reduction has the wrong sign. This cannot */
941: /* happen in infinite precision arithmetic. Step should */
942: /* be rejected! */
943: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
944: }
945: else {
946: if (PetscIsInfOrNanReal(f_full)) {
947: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
948: }
949: else {
950: actred = fold - f_full;
951: prered = -prered;
952: if ((PetscAbsScalar(actred) <= nlsP->epsilon) &&
953: (PetscAbsScalar(prered) <= nlsP->epsilon)) {
954: kappa = 1.0;
955: }
956: else {
957: kappa = actred / prered;
958: }
960: tau_1 = nlsP->theta * gdx / (nlsP->theta * gdx - (1.0 - nlsP->theta) * prered + actred);
961: tau_2 = nlsP->theta * gdx / (nlsP->theta * gdx + (1.0 + nlsP->theta) * prered - actred);
962: tau_min = PetscMin(tau_1, tau_2);
963: tau_max = PetscMax(tau_1, tau_2);
965: if (kappa >= 1.0 - nlsP->mu1) {
966: /* Great agreement */
967: if (tau_max < 1.0) {
968: tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
969: }
970: else if (tau_max > nlsP->gamma4) {
971: tao->trust = PetscMax(tao->trust, nlsP->gamma4 * norm_d);
972: }
973: else {
974: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
975: }
976: }
977: else if (kappa >= 1.0 - nlsP->mu2) {
978: /* Good agreement */
980: if (tau_max < nlsP->gamma2) {
981: tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
982: }
983: else if (tau_max > nlsP->gamma3) {
984: tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
985: }
986: else if (tau_max < 1.0) {
987: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
988: }
989: else {
990: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
991: }
992: }
993: else {
994: /* Not good agreement */
995: if (tau_min > 1.0) {
996: tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
997: }
998: else if (tau_max < nlsP->gamma1) {
999: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
1000: }
1001: else if ((tau_min < nlsP->gamma1) && (tau_max >= 1.0)) {
1002: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
1003: }
1004: else if ((tau_1 >= nlsP->gamma1) && (tau_1 < 1.0) &&
1005: ((tau_2 < nlsP->gamma1) || (tau_2 >= 1.0))) {
1006: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
1007: }
1008: else if ((tau_2 >= nlsP->gamma1) && (tau_2 < 1.0) &&
1009: ((tau_1 < nlsP->gamma1) || (tau_2 >= 1.0))) {
1010: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
1011: }
1012: else {
1013: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
1014: }
1015: }
1016: }
1017: }
1018: }
1019: else {
1020: /* Newton step was not good; reduce the radius */
1021: tao->trust = nlsP->gamma1 * PetscMin(norm_d, tao->trust);
1022: }
1023: break;
1024: }
1026: /* The radius may have been increased; modify if it is too large */
1027: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
1028: }
1030: /* Check for termination */
1031: VecNorm(tao->gradient, NORM_2, &gnorm);
1032: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
1033: SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number");
1034: }
1035: needH = 1;
1037: TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);
1038: }
1039: return(0);
1040: }
1042: /* ---------------------------------------------------------- */
1045: static PetscErrorCode TaoSetUp_NLS(TaoSolver tao)
1046: {
1047: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
1052: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
1053: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
1054: if (!nlsP->W) {VecDuplicate(tao->solution,&nlsP->W); }
1055: if (!nlsP->D) {VecDuplicate(tao->solution,&nlsP->D); }
1056: if (!nlsP->Xold) {VecDuplicate(tao->solution,&nlsP->Xold); }
1057: if (!nlsP->Gold) {VecDuplicate(tao->solution,&nlsP->Gold); }
1059: nlsP->Diag = 0;
1060: nlsP->M = 0;
1062: return(0);
1063: }
1065: /*------------------------------------------------------------*/
1068: static PetscErrorCode TaoDestroy_NLS(TaoSolver tao)
1069: {
1070: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
1074: if (tao->setupcalled) {
1075: VecDestroy(&nlsP->D);
1076: VecDestroy(&nlsP->W);
1077: VecDestroy(&nlsP->Xold);
1078: VecDestroy(&nlsP->Gold);
1079: }
1080: if (nlsP->Diag) {
1081: VecDestroy(&nlsP->Diag);
1082: }
1083: if (nlsP->M) {
1084: MatDestroy(&nlsP->M);
1085: }
1087: PetscFree(tao->data);
1088: tao->data = PETSC_NULL;
1090: return(0);
1091: }
1093: /*------------------------------------------------------------*/
1096: static PetscErrorCode TaoSetFromOptions_NLS(TaoSolver tao)
1097: {
1098: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
1102: PetscOptionsHead("Newton line search method for unconstrained optimization");
1103: PetscOptionsEList("-tao_nls_ksp_type", "ksp type", "", NLS_KSP, NLS_KSP_TYPES, NLS_KSP[nlsP->ksp_type], &nlsP->ksp_type, 0);
1104: PetscOptionsEList("-tao_nls_pc_type", "pc type", "", NLS_PC, NLS_PC_TYPES, NLS_PC[nlsP->pc_type], &nlsP->pc_type, 0);
1105: 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);
1106: PetscOptionsEList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[nlsP->init_type], &nlsP->init_type, 0);
1107: PetscOptionsEList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[nlsP->update_type], &nlsP->update_type, 0);
1108: PetscOptionsReal("-tao_nls_sval", "perturbation starting value", "", nlsP->sval, &nlsP->sval, 0);
1109: PetscOptionsReal("-tao_nls_imin", "minimum initial perturbation", "", nlsP->imin, &nlsP->imin, 0);
1110: PetscOptionsReal("-tao_nls_imax", "maximum initial perturbation", "", nlsP->imax, &nlsP->imax, 0);
1111: PetscOptionsReal("-tao_nls_imfac", "initial merit factor", "", nlsP->imfac, &nlsP->imfac, 0);
1112: PetscOptionsReal("-tao_nls_pmin", "minimum perturbation", "", nlsP->pmin, &nlsP->pmin, 0);
1113: PetscOptionsReal("-tao_nls_pmax", "maximum perturbation", "", nlsP->pmax, &nlsP->pmax, 0);
1114: PetscOptionsReal("-tao_nls_pgfac", "growth factor", "", nlsP->pgfac, &nlsP->pgfac, 0);
1115: PetscOptionsReal("-tao_nls_psfac", "shrink factor", "", nlsP->psfac, &nlsP->psfac, 0);
1116: PetscOptionsReal("-tao_nls_pmgfac", "merit growth factor", "", nlsP->pmgfac, &nlsP->pmgfac, 0);
1117: PetscOptionsReal("-tao_nls_pmsfac", "merit shrink factor", "", nlsP->pmsfac, &nlsP->pmsfac, 0);
1118: PetscOptionsReal("-tao_nls_eta1", "poor steplength; reduce radius", "", nlsP->eta1, &nlsP->eta1, 0);
1119: PetscOptionsReal("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", nlsP->eta2, &nlsP->eta2, 0);
1120: PetscOptionsReal("-tao_nls_eta3", "good steplength; increase radius", "", nlsP->eta3, &nlsP->eta3, 0);
1121: PetscOptionsReal("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", nlsP->eta4, &nlsP->eta4, 0);
1122: PetscOptionsReal("-tao_nls_alpha1", "", "", nlsP->alpha1, &nlsP->alpha1, 0);
1123: PetscOptionsReal("-tao_nls_alpha2", "", "", nlsP->alpha2, &nlsP->alpha2, 0);
1124: PetscOptionsReal("-tao_nls_alpha3", "", "", nlsP->alpha3, &nlsP->alpha3, 0);
1125: PetscOptionsReal("-tao_nls_alpha4", "", "", nlsP->alpha4, &nlsP->alpha4, 0);
1126: PetscOptionsReal("-tao_nls_alpha5", "", "", nlsP->alpha5, &nlsP->alpha5, 0);
1127: PetscOptionsReal("-tao_nls_nu1", "poor steplength; reduce radius", "", nlsP->nu1, &nlsP->nu1, 0);
1128: PetscOptionsReal("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", nlsP->nu2, &nlsP->nu2, 0);
1129: PetscOptionsReal("-tao_nls_nu3", "good steplength; increase radius", "", nlsP->nu3, &nlsP->nu3, 0);
1130: PetscOptionsReal("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", nlsP->nu4, &nlsP->nu4, 0);
1131: PetscOptionsReal("-tao_nls_omega1", "", "", nlsP->omega1, &nlsP->omega1, 0);
1132: PetscOptionsReal("-tao_nls_omega2", "", "", nlsP->omega2, &nlsP->omega2, 0);
1133: PetscOptionsReal("-tao_nls_omega3", "", "", nlsP->omega3, &nlsP->omega3, 0);
1134: PetscOptionsReal("-tao_nls_omega4", "", "", nlsP->omega4, &nlsP->omega4, 0);
1135: PetscOptionsReal("-tao_nls_omega5", "", "", nlsP->omega5, &nlsP->omega5, 0);
1136: PetscOptionsReal("-tao_nls_mu1_i", "", "", nlsP->mu1_i, &nlsP->mu1_i, 0);
1137: PetscOptionsReal("-tao_nls_mu2_i", "", "", nlsP->mu2_i, &nlsP->mu2_i, 0);
1138: PetscOptionsReal("-tao_nls_gamma1_i", "", "", nlsP->gamma1_i, &nlsP->gamma1_i, 0);
1139: PetscOptionsReal("-tao_nls_gamma2_i", "", "", nlsP->gamma2_i, &nlsP->gamma2_i, 0);
1140: PetscOptionsReal("-tao_nls_gamma3_i", "", "", nlsP->gamma3_i, &nlsP->gamma3_i, 0);
1141: PetscOptionsReal("-tao_nls_gamma4_i", "", "", nlsP->gamma4_i, &nlsP->gamma4_i, 0);
1142: PetscOptionsReal("-tao_nls_theta_i", "", "", nlsP->theta_i, &nlsP->theta_i, 0);
1143: PetscOptionsReal("-tao_nls_mu1", "", "", nlsP->mu1, &nlsP->mu1, 0);
1144: PetscOptionsReal("-tao_nls_mu2", "", "", nlsP->mu2, &nlsP->mu2, 0);
1145: PetscOptionsReal("-tao_nls_gamma1", "", "", nlsP->gamma1, &nlsP->gamma1, 0);
1146: PetscOptionsReal("-tao_nls_gamma2", "", "", nlsP->gamma2, &nlsP->gamma2, 0);
1147: PetscOptionsReal("-tao_nls_gamma3", "", "", nlsP->gamma3, &nlsP->gamma3, 0);
1148: PetscOptionsReal("-tao_nls_gamma4", "", "", nlsP->gamma4, &nlsP->gamma4, 0);
1149: PetscOptionsReal("-tao_nls_theta", "", "", nlsP->theta, &nlsP->theta, 0);
1150: PetscOptionsReal("-tao_nls_min_radius", "lower bound on initial radius", "", nlsP->min_radius, &nlsP->min_radius, 0);
1151: PetscOptionsReal("-tao_nls_max_radius", "upper bound on radius", "", nlsP->max_radius, &nlsP->max_radius, 0);
1152: PetscOptionsReal("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", nlsP->epsilon, &nlsP->epsilon, 0);
1153: PetscOptionsTail();
1154: TaoLineSearchSetFromOptions(tao->linesearch);
1155: KSPSetFromOptions(tao->ksp);
1156: return(0);
1157: }
1160: /*------------------------------------------------------------*/
1163: static PetscErrorCode TaoView_NLS(TaoSolver tao, PetscViewer viewer)
1164: {
1165: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
1166: PetscInt nrejects;
1167: PetscBool isascii;
1172: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1173: if (isascii) {
1174: PetscViewerASCIIPushTab(viewer);
1175: if (NLS_PC_BFGS == nlsP->pc_type && nlsP->M) {
1176: MatLMVMGetRejects(nlsP->M,&nrejects);
1177: PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n",nrejects);
1178: }
1179: PetscViewerASCIIPrintf(viewer, "Newton steps: %D\n", nlsP->newt);
1180: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", nlsP->bfgs);
1181: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", nlsP->sgrad);
1182: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", nlsP->grad);
1184: PetscViewerASCIIPrintf(viewer, "nls ksp atol: %D\n", nlsP->ksp_atol);
1185: PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %D\n", nlsP->ksp_rtol);
1186: PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %D\n", nlsP->ksp_ctol);
1187: PetscViewerASCIIPrintf(viewer, "nls ksp negc: %D\n", nlsP->ksp_negc);
1188: PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %D\n", nlsP->ksp_dtol);
1189: PetscViewerASCIIPrintf(viewer, "nls ksp iter: %D\n", nlsP->ksp_iter);
1190: PetscViewerASCIIPrintf(viewer, "nls ksp othr: %D\n", nlsP->ksp_othr);
1191: PetscViewerASCIIPopTab(viewer);
1192: } else {
1193: SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO NLS",((PetscObject)viewer)->type_name);
1194: }
1195: return(0);
1196: }
1198: /* ---------------------------------------------------------- */
1202: PetscErrorCode TaoCreate_NLS(TaoSolver tao)
1203: {
1204: TAO_NLS *nlsP;
1205: const char *morethuente_type = TAOLINESEARCH_MT;
1209: PetscNewLog(tao,TAO_NLS,&nlsP);
1211: tao->ops->setup = TaoSetUp_NLS;
1212: tao->ops->solve = TaoSolve_NLS;
1213: tao->ops->view = TaoView_NLS;
1214: tao->ops->setfromoptions = TaoSetFromOptions_NLS;
1215: tao->ops->destroy = TaoDestroy_NLS;
1217: tao->max_it = 50;
1218: tao->fatol = 1e-10;
1219: tao->frtol = 1e-10;
1220: tao->data = (void*)nlsP;
1221: tao->trust0 = 100.0;
1223: nlsP->sval = 0.0;
1224: nlsP->imin = 1.0e-4;
1225: nlsP->imax = 1.0e+2;
1226: nlsP->imfac = 1.0e-1;
1228: nlsP->pmin = 1.0e-12;
1229: nlsP->pmax = 1.0e+2;
1230: nlsP->pgfac = 1.0e+1;
1231: nlsP->psfac = 4.0e-1;
1232: nlsP->pmgfac = 1.0e-1;
1233: nlsP->pmsfac = 1.0e-1;
1235: /* Default values for trust-region radius update based on steplength */
1236: nlsP->nu1 = 0.25;
1237: nlsP->nu2 = 0.50;
1238: nlsP->nu3 = 1.00;
1239: nlsP->nu4 = 1.25;
1241: nlsP->omega1 = 0.25;
1242: nlsP->omega2 = 0.50;
1243: nlsP->omega3 = 1.00;
1244: nlsP->omega4 = 2.00;
1245: nlsP->omega5 = 4.00;
1247: /* Default values for trust-region radius update based on reduction */
1248: nlsP->eta1 = 1.0e-4;
1249: nlsP->eta2 = 0.25;
1250: nlsP->eta3 = 0.50;
1251: nlsP->eta4 = 0.90;
1253: nlsP->alpha1 = 0.25;
1254: nlsP->alpha2 = 0.50;
1255: nlsP->alpha3 = 1.00;
1256: nlsP->alpha4 = 2.00;
1257: nlsP->alpha5 = 4.00;
1259: /* Default values for trust-region radius update based on interpolation */
1260: nlsP->mu1 = 0.10;
1261: nlsP->mu2 = 0.50;
1263: nlsP->gamma1 = 0.25;
1264: nlsP->gamma2 = 0.50;
1265: nlsP->gamma3 = 2.00;
1266: nlsP->gamma4 = 4.00;
1268: nlsP->theta = 0.05;
1270: /* Default values for trust region initialization based on interpolation */
1271: nlsP->mu1_i = 0.35;
1272: nlsP->mu2_i = 0.50;
1274: nlsP->gamma1_i = 0.0625;
1275: nlsP->gamma2_i = 0.5;
1276: nlsP->gamma3_i = 2.0;
1277: nlsP->gamma4_i = 5.0;
1278:
1279: nlsP->theta_i = 0.25;
1281: /* Remaining parameters */
1282: nlsP->min_radius = 1.0e-10;
1283: nlsP->max_radius = 1.0e10;
1284: nlsP->epsilon = 1.0e-6;
1286: nlsP->ksp_type = NLS_KSP_STCG;
1287: nlsP->pc_type = NLS_PC_BFGS;
1288: nlsP->bfgs_scale_type = BFGS_SCALE_PHESS;
1289: nlsP->init_type = NLS_INIT_INTERPOLATION;
1290: nlsP->update_type = NLS_UPDATE_STEP;
1292: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
1293: TaoLineSearchSetType(tao->linesearch,morethuente_type);
1294: TaoLineSearchUseTaoSolverRoutines(tao->linesearch,tao);
1296: /* Set linear solver to default for symmetric matrices */
1297: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
1299: return(0);
1300: }
1309: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x)
1310: {
1312: Mat M;
1317: PCShellGetContext(pc,(void**)&M);
1318: MatLMVMSolve(M, b, x);
1319: return(0);
1320: }