Actual source code: nls.c
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/unconstrained/impls/nls/nlsimpl.h>
4: #include <petscksp.h>
6: #define NLS_INIT_CONSTANT 0
7: #define NLS_INIT_DIRECTION 1
8: #define NLS_INIT_INTERPOLATION 2
9: #define NLS_INIT_TYPES 3
11: #define NLS_UPDATE_STEP 0
12: #define NLS_UPDATE_REDUCTION 1
13: #define NLS_UPDATE_INTERPOLATION 2
14: #define NLS_UPDATE_TYPES 3
16: static const char *NLS_INIT[64] = {"constant", "direction", "interpolation"};
18: static const char *NLS_UPDATE[64] = {"step", "reduction", "interpolation"};
20: /*
21: Implements Newton's Method with a line search approach for solving
22: unconstrained minimization problems. A More'-Thuente line search
23: is used to guarantee that the bfgs preconditioner remains positive
24: definite.
26: The method can shift the Hessian matrix. The shifting procedure is
27: adapted from the PATH algorithm for solving complementarity
28: problems.
30: The linear system solve should be done with a conjugate gradient
31: method, although any method can be used.
32: */
34: #define NLS_NEWTON 0
35: #define NLS_BFGS 1
36: #define NLS_GRADIENT 2
38: static PetscErrorCode TaoSolve_NLS(Tao tao)
39: {
40: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
41: KSPType ksp_type;
42: PetscBool is_nash,is_stcg,is_gltr,is_bfgs,is_jacobi,is_symmetric,sym_set;
43: KSPConvergedReason ksp_reason;
44: PC pc;
45: TaoLineSearchConvergedReason ls_reason;
47: PetscReal fmin, ftrial, f_full, prered, actred, kappa, sigma;
48: PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
49: PetscReal f, fold, gdx, gnorm, pert;
50: PetscReal step = 1.0;
51: PetscReal norm_d = 0.0, e_min;
53: PetscInt stepType;
54: PetscInt bfgsUpdates = 0;
55: PetscInt n,N,kspits;
56: PetscInt needH = 1;
58: PetscInt i_max = 5;
59: PetscInt j_max = 1;
60: PetscInt i, j;
62: if (tao->XL || tao->XU || tao->ops->computebounds) {
63: PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by nls algorithm\n");
64: }
66: /* Initialized variables */
67: pert = nlsP->sval;
69: /* Number of times ksp stopped because of these reasons */
70: nlsP->ksp_atol = 0;
71: nlsP->ksp_rtol = 0;
72: nlsP->ksp_dtol = 0;
73: nlsP->ksp_ctol = 0;
74: nlsP->ksp_negc = 0;
75: nlsP->ksp_iter = 0;
76: nlsP->ksp_othr = 0;
78: /* Initialize trust-region radius when using nash, stcg, or gltr
79: Command automatically ignored for other methods
80: Will be reset during the first iteration
81: */
82: KSPGetType(tao->ksp,&ksp_type);
83: PetscStrcmp(ksp_type,KSPNASH,&is_nash);
84: PetscStrcmp(ksp_type,KSPSTCG,&is_stcg);
85: PetscStrcmp(ksp_type,KSPGLTR,&is_gltr);
87: KSPCGSetRadius(tao->ksp,nlsP->max_radius);
89: if (is_nash || is_stcg || is_gltr) {
91: tao->trust = tao->trust0;
92: tao->trust = PetscMax(tao->trust, nlsP->min_radius);
93: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
94: }
96: /* Check convergence criteria */
97: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
98: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
101: tao->reason = TAO_CONTINUE_ITERATING;
102: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
103: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
104: (*tao->ops->convergencetest)(tao,tao->cnvP);
105: if (tao->reason != TAO_CONTINUE_ITERATING) return 0;
107: /* Allocate the vectors needed for the BFGS approximation */
108: KSPGetPC(tao->ksp, &pc);
109: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);
110: PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);
111: if (is_bfgs) {
112: nlsP->bfgs_pre = pc;
113: PCLMVMGetMatLMVM(nlsP->bfgs_pre, &nlsP->M);
114: VecGetLocalSize(tao->solution, &n);
115: VecGetSize(tao->solution, &N);
116: MatSetSizes(nlsP->M, n, n, N, N);
117: MatLMVMAllocate(nlsP->M, tao->solution, tao->gradient);
118: MatIsSymmetricKnown(nlsP->M, &sym_set, &is_symmetric);
120: } else if (is_jacobi) {
121: PCJacobiSetUseAbs(pc,PETSC_TRUE);
122: }
124: /* Initialize trust-region radius. The initialization is only performed
125: when we are using Nash, Steihaug-Toint or the Generalized Lanczos method. */
126: if (is_nash || is_stcg || is_gltr) {
127: switch (nlsP->init_type) {
128: case NLS_INIT_CONSTANT:
129: /* Use the initial radius specified */
130: break;
132: case NLS_INIT_INTERPOLATION:
133: /* Use the initial radius specified */
134: max_radius = 0.0;
136: for (j = 0; j < j_max; ++j) {
137: fmin = f;
138: sigma = 0.0;
140: if (needH) {
141: TaoComputeHessian(tao, tao->solution,tao->hessian,tao->hessian_pre);
142: needH = 0;
143: }
145: for (i = 0; i < i_max; ++i) {
146: VecCopy(tao->solution,nlsP->W);
147: VecAXPY(nlsP->W,-tao->trust/gnorm,tao->gradient);
148: TaoComputeObjective(tao, nlsP->W, &ftrial);
149: if (PetscIsInfOrNanReal(ftrial)) {
150: tau = nlsP->gamma1_i;
151: } else {
152: if (ftrial < fmin) {
153: fmin = ftrial;
154: sigma = -tao->trust / gnorm;
155: }
157: MatMult(tao->hessian, tao->gradient, nlsP->D);
158: VecDot(tao->gradient, nlsP->D, &prered);
160: prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
161: actred = f - ftrial;
162: if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) {
163: kappa = 1.0;
164: } else {
165: kappa = actred / prered;
166: }
168: tau_1 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust + (1.0 - nlsP->theta_i) * prered - actred);
169: tau_2 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust - (1.0 + nlsP->theta_i) * prered + actred);
170: tau_min = PetscMin(tau_1, tau_2);
171: tau_max = PetscMax(tau_1, tau_2);
173: if (PetscAbsScalar(kappa - (PetscReal)1.0) <= nlsP->mu1_i) {
174: /* Great agreement */
175: max_radius = PetscMax(max_radius, tao->trust);
177: if (tau_max < 1.0) {
178: tau = nlsP->gamma3_i;
179: } else if (tau_max > nlsP->gamma4_i) {
180: tau = nlsP->gamma4_i;
181: } else if (tau_1 >= 1.0 && tau_1 <= nlsP->gamma4_i && tau_2 < 1.0) {
182: tau = tau_1;
183: } else if (tau_2 >= 1.0 && tau_2 <= nlsP->gamma4_i && tau_1 < 1.0) {
184: tau = tau_2;
185: } else {
186: tau = tau_max;
187: }
188: } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= nlsP->mu2_i) {
189: /* Good agreement */
190: max_radius = PetscMax(max_radius, tao->trust);
192: if (tau_max < nlsP->gamma2_i) {
193: tau = nlsP->gamma2_i;
194: } else if (tau_max > nlsP->gamma3_i) {
195: tau = nlsP->gamma3_i;
196: } else {
197: tau = tau_max;
198: }
199: } else {
200: /* Not good agreement */
201: if (tau_min > 1.0) {
202: tau = nlsP->gamma2_i;
203: } else if (tau_max < nlsP->gamma1_i) {
204: tau = nlsP->gamma1_i;
205: } else if ((tau_min < nlsP->gamma1_i) && (tau_max >= 1.0)) {
206: tau = nlsP->gamma1_i;
207: } else if ((tau_1 >= nlsP->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < nlsP->gamma1_i) || (tau_2 >= 1.0))) {
208: tau = tau_1;
209: } else if ((tau_2 >= nlsP->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < nlsP->gamma1_i) || (tau_2 >= 1.0))) {
210: tau = tau_2;
211: } else {
212: tau = tau_max;
213: }
214: }
215: }
216: tao->trust = tau * tao->trust;
217: }
219: if (fmin < f) {
220: f = fmin;
221: VecAXPY(tao->solution,sigma,tao->gradient);
222: TaoComputeGradient(tao,tao->solution,tao->gradient);
224: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
226: needH = 1;
228: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
229: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
230: (*tao->ops->convergencetest)(tao,tao->cnvP);
231: if (tao->reason != TAO_CONTINUE_ITERATING) return 0;
232: }
233: }
234: tao->trust = PetscMax(tao->trust, max_radius);
236: /* Modify the radius if it is too large or small */
237: tao->trust = PetscMax(tao->trust, nlsP->min_radius);
238: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
239: break;
241: default:
242: /* Norm of the first direction will initialize radius */
243: tao->trust = 0.0;
244: break;
245: }
246: }
248: /* Set counter for gradient/reset steps*/
249: nlsP->newt = 0;
250: nlsP->bfgs = 0;
251: nlsP->grad = 0;
253: /* Have not converged; continue with Newton method */
254: while (tao->reason == TAO_CONTINUE_ITERATING) {
255: /* Call general purpose update function */
256: if (tao->ops->update) {
257: (*tao->ops->update)(tao, tao->niter, tao->user_update);
258: }
259: ++tao->niter;
260: tao->ksp_its = 0;
262: /* Compute the Hessian */
263: if (needH) {
264: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
265: }
267: /* Shift the Hessian matrix */
268: if (pert > 0) {
269: MatShift(tao->hessian, pert);
270: if (tao->hessian != tao->hessian_pre) {
271: MatShift(tao->hessian_pre, pert);
272: }
273: }
275: if (nlsP->bfgs_pre) {
276: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
277: ++bfgsUpdates;
278: }
280: /* Solve the Newton system of equations */
281: KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre);
282: if (is_nash || is_stcg || is_gltr) {
283: KSPCGSetRadius(tao->ksp,nlsP->max_radius);
284: KSPSolve(tao->ksp, tao->gradient, nlsP->D);
285: KSPGetIterationNumber(tao->ksp,&kspits);
286: tao->ksp_its += kspits;
287: tao->ksp_tot_its += kspits;
288: KSPCGGetNormD(tao->ksp,&norm_d);
290: if (0.0 == tao->trust) {
291: /* Radius was uninitialized; use the norm of the direction */
292: if (norm_d > 0.0) {
293: tao->trust = norm_d;
295: /* Modify the radius if it is too large or small */
296: tao->trust = PetscMax(tao->trust, nlsP->min_radius);
297: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
298: } else {
299: /* The direction was bad; set radius to default value and re-solve
300: the trust-region subproblem to get a direction */
301: tao->trust = tao->trust0;
303: /* Modify the radius if it is too large or small */
304: tao->trust = PetscMax(tao->trust, nlsP->min_radius);
305: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
307: KSPCGSetRadius(tao->ksp,nlsP->max_radius);
308: KSPSolve(tao->ksp, tao->gradient, nlsP->D);
309: KSPGetIterationNumber(tao->ksp,&kspits);
310: tao->ksp_its += kspits;
311: tao->ksp_tot_its += kspits;
312: KSPCGGetNormD(tao->ksp,&norm_d);
315: }
316: }
317: } else {
318: KSPSolve(tao->ksp, tao->gradient, nlsP->D);
319: KSPGetIterationNumber(tao->ksp, &kspits);
320: tao->ksp_its += kspits;
321: tao->ksp_tot_its += kspits;
322: }
323: VecScale(nlsP->D, -1.0);
324: KSPGetConvergedReason(tao->ksp, &ksp_reason);
325: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (nlsP->bfgs_pre)) {
326: /* Preconditioner is numerically indefinite; reset the
327: approximate if using BFGS preconditioning. */
328: MatLMVMReset(nlsP->M, PETSC_FALSE);
329: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
330: bfgsUpdates = 1;
331: }
333: if (KSP_CONVERGED_ATOL == ksp_reason) {
334: ++nlsP->ksp_atol;
335: } else if (KSP_CONVERGED_RTOL == ksp_reason) {
336: ++nlsP->ksp_rtol;
337: } else if (KSP_CONVERGED_CG_CONSTRAINED == ksp_reason) {
338: ++nlsP->ksp_ctol;
339: } else if (KSP_CONVERGED_CG_NEG_CURVE == ksp_reason) {
340: ++nlsP->ksp_negc;
341: } else if (KSP_DIVERGED_DTOL == ksp_reason) {
342: ++nlsP->ksp_dtol;
343: } else if (KSP_DIVERGED_ITS == ksp_reason) {
344: ++nlsP->ksp_iter;
345: } else {
346: ++nlsP->ksp_othr;
347: }
349: /* Check for success (descent direction) */
350: VecDot(nlsP->D, tao->gradient, &gdx);
351: if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
352: /* Newton step is not descent or direction produced Inf or NaN
353: Update the perturbation for next time */
354: if (pert <= 0.0) {
355: /* Initialize the perturbation */
356: pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
357: if (is_gltr) {
358: KSPGLTRGetMinEig(tao->ksp,&e_min);
359: pert = PetscMax(pert, -e_min);
360: }
361: } else {
362: /* Increase the perturbation */
363: pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
364: }
366: if (!nlsP->bfgs_pre) {
367: /* We don't have the bfgs matrix around and updated
368: Must use gradient direction in this case */
369: VecCopy(tao->gradient, nlsP->D);
370: VecScale(nlsP->D, -1.0);
371: ++nlsP->grad;
372: stepType = NLS_GRADIENT;
373: } else {
374: /* Attempt to use the BFGS direction */
375: MatSolve(nlsP->M, tao->gradient, nlsP->D);
376: VecScale(nlsP->D, -1.0);
378: /* Check for success (descent direction) */
379: VecDot(tao->gradient, nlsP->D, &gdx);
380: if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
381: /* BFGS direction is not descent or direction produced not a number
382: We can assert bfgsUpdates > 1 in this case because
383: the first solve produces the scaled gradient direction,
384: which is guaranteed to be descent */
386: /* Use steepest descent direction (scaled) */
387: MatLMVMReset(nlsP->M, PETSC_FALSE);
388: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
389: MatSolve(nlsP->M, tao->gradient, nlsP->D);
390: VecScale(nlsP->D, -1.0);
392: bfgsUpdates = 1;
393: ++nlsP->grad;
394: stepType = NLS_GRADIENT;
395: } else {
396: MatLMVMGetUpdateCount(nlsP->M, &bfgsUpdates);
397: if (1 == bfgsUpdates) {
398: /* The first BFGS direction is always the scaled gradient */
399: ++nlsP->grad;
400: stepType = NLS_GRADIENT;
401: } else {
402: ++nlsP->bfgs;
403: stepType = NLS_BFGS;
404: }
405: }
406: }
407: } else {
408: /* Computed Newton step is descent */
409: switch (ksp_reason) {
410: case KSP_DIVERGED_NANORINF:
411: case KSP_DIVERGED_BREAKDOWN:
412: case KSP_DIVERGED_INDEFINITE_MAT:
413: case KSP_DIVERGED_INDEFINITE_PC:
414: case KSP_CONVERGED_CG_NEG_CURVE:
415: /* Matrix or preconditioner is indefinite; increase perturbation */
416: if (pert <= 0.0) {
417: /* Initialize the perturbation */
418: pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
419: if (is_gltr) {
420: KSPGLTRGetMinEig(tao->ksp, &e_min);
421: pert = PetscMax(pert, -e_min);
422: }
423: } else {
424: /* Increase the perturbation */
425: pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
426: }
427: break;
429: default:
430: /* Newton step computation is good; decrease perturbation */
431: pert = PetscMin(nlsP->psfac * pert, nlsP->pmsfac * gnorm);
432: if (pert < nlsP->pmin) {
433: pert = 0.0;
434: }
435: break;
436: }
438: ++nlsP->newt;
439: stepType = NLS_NEWTON;
440: }
442: /* Perform the linesearch */
443: fold = f;
444: VecCopy(tao->solution, nlsP->Xold);
445: VecCopy(tao->gradient, nlsP->Gold);
447: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason);
448: TaoAddLineSearchCounts(tao);
450: while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NLS_GRADIENT) { /* Linesearch failed */
451: f = fold;
452: VecCopy(nlsP->Xold, tao->solution);
453: VecCopy(nlsP->Gold, tao->gradient);
455: switch(stepType) {
456: case NLS_NEWTON:
457: /* Failed to obtain acceptable iterate with Newton 1step
458: Update the perturbation for next time */
459: if (pert <= 0.0) {
460: /* Initialize the perturbation */
461: pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
462: if (is_gltr) {
463: KSPGLTRGetMinEig(tao->ksp,&e_min);
464: pert = PetscMax(pert, -e_min);
465: }
466: } else {
467: /* Increase the perturbation */
468: pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
469: }
471: if (!nlsP->bfgs_pre) {
472: /* We don't have the bfgs matrix around and being updated
473: Must use gradient direction in this case */
474: VecCopy(tao->gradient, nlsP->D);
475: ++nlsP->grad;
476: stepType = NLS_GRADIENT;
477: } else {
478: /* Attempt to use the BFGS direction */
479: MatSolve(nlsP->M, tao->gradient, nlsP->D);
480: /* Check for success (descent direction) */
481: VecDot(tao->solution, nlsP->D, &gdx);
482: if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
483: /* BFGS direction is not descent or direction produced not a number
484: We can assert bfgsUpdates > 1 in this case
485: Use steepest descent direction (scaled) */
486: MatLMVMReset(nlsP->M, PETSC_FALSE);
487: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
488: MatSolve(nlsP->M, tao->gradient, nlsP->D);
490: bfgsUpdates = 1;
491: ++nlsP->grad;
492: stepType = NLS_GRADIENT;
493: } else {
494: if (1 == bfgsUpdates) {
495: /* The first BFGS direction is always the scaled gradient */
496: ++nlsP->grad;
497: stepType = NLS_GRADIENT;
498: } else {
499: ++nlsP->bfgs;
500: stepType = NLS_BFGS;
501: }
502: }
503: }
504: break;
506: case NLS_BFGS:
507: /* Can only enter if pc_type == NLS_PC_BFGS
508: Failed to obtain acceptable iterate with BFGS step
509: Attempt to use the scaled gradient direction */
510: MatLMVMReset(nlsP->M, PETSC_FALSE);
511: MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
512: MatSolve(nlsP->M, tao->gradient, nlsP->D);
514: bfgsUpdates = 1;
515: ++nlsP->grad;
516: stepType = NLS_GRADIENT;
517: break;
518: }
519: VecScale(nlsP->D, -1.0);
521: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason);
522: TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full);
523: TaoAddLineSearchCounts(tao);
524: }
526: if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
527: /* Failed to find an improving point */
528: f = fold;
529: VecCopy(nlsP->Xold, tao->solution);
530: VecCopy(nlsP->Gold, tao->gradient);
531: step = 0.0;
532: tao->reason = TAO_DIVERGED_LS_FAILURE;
533: break;
534: }
536: /* Update trust region radius */
537: if (is_nash || is_stcg || is_gltr) {
538: switch(nlsP->update_type) {
539: case NLS_UPDATE_STEP:
540: if (stepType == NLS_NEWTON) {
541: if (step < nlsP->nu1) {
542: /* Very bad step taken; reduce radius */
543: tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
544: } else if (step < nlsP->nu2) {
545: /* Reasonably bad step taken; reduce radius */
546: tao->trust = nlsP->omega2 * PetscMin(norm_d, tao->trust);
547: } else if (step < nlsP->nu3) {
548: /* Reasonable step was taken; leave radius alone */
549: if (nlsP->omega3 < 1.0) {
550: tao->trust = nlsP->omega3 * PetscMin(norm_d, tao->trust);
551: } else if (nlsP->omega3 > 1.0) {
552: tao->trust = PetscMax(nlsP->omega3 * norm_d, tao->trust);
553: }
554: } else if (step < nlsP->nu4) {
555: /* Full step taken; increase the radius */
556: tao->trust = PetscMax(nlsP->omega4 * norm_d, tao->trust);
557: } else {
558: /* More than full step taken; increase the radius */
559: tao->trust = PetscMax(nlsP->omega5 * norm_d, tao->trust);
560: }
561: } else {
562: /* Newton step was not good; reduce the radius */
563: tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
564: }
565: break;
567: case NLS_UPDATE_REDUCTION:
568: if (stepType == NLS_NEWTON) {
569: /* Get predicted reduction */
571: KSPCGGetObjFcn(tao->ksp,&prered);
572: if (prered >= 0.0) {
573: /* The predicted reduction has the wrong sign. This cannot */
574: /* happen in infinite precision arithmetic. Step should */
575: /* be rejected! */
576: tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
577: } else {
578: if (PetscIsInfOrNanReal(f_full)) {
579: tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
580: } else {
581: /* Compute and actual reduction */
582: actred = fold - f_full;
583: prered = -prered;
584: if ((PetscAbsScalar(actred) <= nlsP->epsilon) &&
585: (PetscAbsScalar(prered) <= nlsP->epsilon)) {
586: kappa = 1.0;
587: } else {
588: kappa = actred / prered;
589: }
591: /* Accept of reject the step and update radius */
592: if (kappa < nlsP->eta1) {
593: /* Very bad step */
594: tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
595: } else if (kappa < nlsP->eta2) {
596: /* Marginal bad step */
597: tao->trust = nlsP->alpha2 * PetscMin(tao->trust, norm_d);
598: } else if (kappa < nlsP->eta3) {
599: /* Reasonable step */
600: if (nlsP->alpha3 < 1.0) {
601: tao->trust = nlsP->alpha3 * PetscMin(norm_d, tao->trust);
602: } else if (nlsP->alpha3 > 1.0) {
603: tao->trust = PetscMax(nlsP->alpha3 * norm_d, tao->trust);
604: }
605: } else if (kappa < nlsP->eta4) {
606: /* Good step */
607: tao->trust = PetscMax(nlsP->alpha4 * norm_d, tao->trust);
608: } else {
609: /* Very good step */
610: tao->trust = PetscMax(nlsP->alpha5 * norm_d, tao->trust);
611: }
612: }
613: }
614: } else {
615: /* Newton step was not good; reduce the radius */
616: tao->trust = nlsP->alpha1 * PetscMin(norm_d, tao->trust);
617: }
618: break;
620: default:
621: if (stepType == NLS_NEWTON) {
622: KSPCGGetObjFcn(tao->ksp,&prered);
623: if (prered >= 0.0) {
624: /* The predicted reduction has the wrong sign. This cannot */
625: /* happen in infinite precision arithmetic. Step should */
626: /* be rejected! */
627: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
628: } else {
629: if (PetscIsInfOrNanReal(f_full)) {
630: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
631: } else {
632: actred = fold - f_full;
633: prered = -prered;
634: if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) {
635: kappa = 1.0;
636: } else {
637: kappa = actred / prered;
638: }
640: tau_1 = nlsP->theta * gdx / (nlsP->theta * gdx - (1.0 - nlsP->theta) * prered + actred);
641: tau_2 = nlsP->theta * gdx / (nlsP->theta * gdx + (1.0 + nlsP->theta) * prered - actred);
642: tau_min = PetscMin(tau_1, tau_2);
643: tau_max = PetscMax(tau_1, tau_2);
645: if (kappa >= 1.0 - nlsP->mu1) {
646: /* Great agreement */
647: if (tau_max < 1.0) {
648: tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
649: } else if (tau_max > nlsP->gamma4) {
650: tao->trust = PetscMax(tao->trust, nlsP->gamma4 * norm_d);
651: } else {
652: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
653: }
654: } else if (kappa >= 1.0 - nlsP->mu2) {
655: /* Good agreement */
657: if (tau_max < nlsP->gamma2) {
658: tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
659: } else if (tau_max > nlsP->gamma3) {
660: tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
661: } else if (tau_max < 1.0) {
662: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
663: } else {
664: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
665: }
666: } else {
667: /* Not good agreement */
668: if (tau_min > 1.0) {
669: tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
670: } else if (tau_max < nlsP->gamma1) {
671: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
672: } else if ((tau_min < nlsP->gamma1) && (tau_max >= 1.0)) {
673: tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
674: } else if ((tau_1 >= nlsP->gamma1) && (tau_1 < 1.0) && ((tau_2 < nlsP->gamma1) || (tau_2 >= 1.0))) {
675: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
676: } else if ((tau_2 >= nlsP->gamma1) && (tau_2 < 1.0) && ((tau_1 < nlsP->gamma1) || (tau_2 >= 1.0))) {
677: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
678: } else {
679: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
680: }
681: }
682: }
683: }
684: } else {
685: /* Newton step was not good; reduce the radius */
686: tao->trust = nlsP->gamma1 * PetscMin(norm_d, tao->trust);
687: }
688: break;
689: }
691: /* The radius may have been increased; modify if it is too large */
692: tao->trust = PetscMin(tao->trust, nlsP->max_radius);
693: }
695: /* Check for termination */
696: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
698: needH = 1;
699: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
700: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
701: (*tao->ops->convergencetest)(tao,tao->cnvP);
702: }
703: return 0;
704: }
706: /* ---------------------------------------------------------- */
707: static PetscErrorCode TaoSetUp_NLS(Tao tao)
708: {
709: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
711: if (!tao->gradient) VecDuplicate(tao->solution,&tao->gradient);
712: if (!tao->stepdirection) VecDuplicate(tao->solution,&tao->stepdirection);
713: if (!nlsP->W) VecDuplicate(tao->solution,&nlsP->W);
714: if (!nlsP->D) VecDuplicate(tao->solution,&nlsP->D);
715: if (!nlsP->Xold) VecDuplicate(tao->solution,&nlsP->Xold);
716: if (!nlsP->Gold) VecDuplicate(tao->solution,&nlsP->Gold);
717: nlsP->M = NULL;
718: nlsP->bfgs_pre = NULL;
719: return 0;
720: }
722: /*------------------------------------------------------------*/
723: static PetscErrorCode TaoDestroy_NLS(Tao tao)
724: {
725: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
727: if (tao->setupcalled) {
728: VecDestroy(&nlsP->D);
729: VecDestroy(&nlsP->W);
730: VecDestroy(&nlsP->Xold);
731: VecDestroy(&nlsP->Gold);
732: }
733: PetscFree(tao->data);
734: return 0;
735: }
737: /*------------------------------------------------------------*/
738: static PetscErrorCode TaoSetFromOptions_NLS(PetscOptionItems *PetscOptionsObject,Tao tao)
739: {
740: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
742: PetscOptionsHead(PetscOptionsObject,"Newton line search method for unconstrained optimization");
743: PetscOptionsEList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[nlsP->init_type], &nlsP->init_type, NULL);
744: PetscOptionsEList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[nlsP->update_type], &nlsP->update_type, NULL);
745: PetscOptionsReal("-tao_nls_sval", "perturbation starting value", "", nlsP->sval, &nlsP->sval,NULL);
746: PetscOptionsReal("-tao_nls_imin", "minimum initial perturbation", "", nlsP->imin, &nlsP->imin,NULL);
747: PetscOptionsReal("-tao_nls_imax", "maximum initial perturbation", "", nlsP->imax, &nlsP->imax,NULL);
748: PetscOptionsReal("-tao_nls_imfac", "initial merit factor", "", nlsP->imfac, &nlsP->imfac,NULL);
749: PetscOptionsReal("-tao_nls_pmin", "minimum perturbation", "", nlsP->pmin, &nlsP->pmin,NULL);
750: PetscOptionsReal("-tao_nls_pmax", "maximum perturbation", "", nlsP->pmax, &nlsP->pmax,NULL);
751: PetscOptionsReal("-tao_nls_pgfac", "growth factor", "", nlsP->pgfac, &nlsP->pgfac,NULL);
752: PetscOptionsReal("-tao_nls_psfac", "shrink factor", "", nlsP->psfac, &nlsP->psfac,NULL);
753: PetscOptionsReal("-tao_nls_pmgfac", "merit growth factor", "", nlsP->pmgfac, &nlsP->pmgfac,NULL);
754: PetscOptionsReal("-tao_nls_pmsfac", "merit shrink factor", "", nlsP->pmsfac, &nlsP->pmsfac,NULL);
755: PetscOptionsReal("-tao_nls_eta1", "poor steplength; reduce radius", "", nlsP->eta1, &nlsP->eta1,NULL);
756: PetscOptionsReal("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", nlsP->eta2, &nlsP->eta2,NULL);
757: PetscOptionsReal("-tao_nls_eta3", "good steplength; increase radius", "", nlsP->eta3, &nlsP->eta3,NULL);
758: PetscOptionsReal("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", nlsP->eta4, &nlsP->eta4,NULL);
759: PetscOptionsReal("-tao_nls_alpha1", "", "", nlsP->alpha1, &nlsP->alpha1,NULL);
760: PetscOptionsReal("-tao_nls_alpha2", "", "", nlsP->alpha2, &nlsP->alpha2,NULL);
761: PetscOptionsReal("-tao_nls_alpha3", "", "", nlsP->alpha3, &nlsP->alpha3,NULL);
762: PetscOptionsReal("-tao_nls_alpha4", "", "", nlsP->alpha4, &nlsP->alpha4,NULL);
763: PetscOptionsReal("-tao_nls_alpha5", "", "", nlsP->alpha5, &nlsP->alpha5,NULL);
764: PetscOptionsReal("-tao_nls_nu1", "poor steplength; reduce radius", "", nlsP->nu1, &nlsP->nu1,NULL);
765: PetscOptionsReal("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", nlsP->nu2, &nlsP->nu2,NULL);
766: PetscOptionsReal("-tao_nls_nu3", "good steplength; increase radius", "", nlsP->nu3, &nlsP->nu3,NULL);
767: PetscOptionsReal("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", nlsP->nu4, &nlsP->nu4,NULL);
768: PetscOptionsReal("-tao_nls_omega1", "", "", nlsP->omega1, &nlsP->omega1,NULL);
769: PetscOptionsReal("-tao_nls_omega2", "", "", nlsP->omega2, &nlsP->omega2,NULL);
770: PetscOptionsReal("-tao_nls_omega3", "", "", nlsP->omega3, &nlsP->omega3,NULL);
771: PetscOptionsReal("-tao_nls_omega4", "", "", nlsP->omega4, &nlsP->omega4,NULL);
772: PetscOptionsReal("-tao_nls_omega5", "", "", nlsP->omega5, &nlsP->omega5,NULL);
773: PetscOptionsReal("-tao_nls_mu1_i", "", "", nlsP->mu1_i, &nlsP->mu1_i,NULL);
774: PetscOptionsReal("-tao_nls_mu2_i", "", "", nlsP->mu2_i, &nlsP->mu2_i,NULL);
775: PetscOptionsReal("-tao_nls_gamma1_i", "", "", nlsP->gamma1_i, &nlsP->gamma1_i,NULL);
776: PetscOptionsReal("-tao_nls_gamma2_i", "", "", nlsP->gamma2_i, &nlsP->gamma2_i,NULL);
777: PetscOptionsReal("-tao_nls_gamma3_i", "", "", nlsP->gamma3_i, &nlsP->gamma3_i,NULL);
778: PetscOptionsReal("-tao_nls_gamma4_i", "", "", nlsP->gamma4_i, &nlsP->gamma4_i,NULL);
779: PetscOptionsReal("-tao_nls_theta_i", "", "", nlsP->theta_i, &nlsP->theta_i,NULL);
780: PetscOptionsReal("-tao_nls_mu1", "", "", nlsP->mu1, &nlsP->mu1,NULL);
781: PetscOptionsReal("-tao_nls_mu2", "", "", nlsP->mu2, &nlsP->mu2,NULL);
782: PetscOptionsReal("-tao_nls_gamma1", "", "", nlsP->gamma1, &nlsP->gamma1,NULL);
783: PetscOptionsReal("-tao_nls_gamma2", "", "", nlsP->gamma2, &nlsP->gamma2,NULL);
784: PetscOptionsReal("-tao_nls_gamma3", "", "", nlsP->gamma3, &nlsP->gamma3,NULL);
785: PetscOptionsReal("-tao_nls_gamma4", "", "", nlsP->gamma4, &nlsP->gamma4,NULL);
786: PetscOptionsReal("-tao_nls_theta", "", "", nlsP->theta, &nlsP->theta,NULL);
787: PetscOptionsReal("-tao_nls_min_radius", "lower bound on initial radius", "", nlsP->min_radius, &nlsP->min_radius,NULL);
788: PetscOptionsReal("-tao_nls_max_radius", "upper bound on radius", "", nlsP->max_radius, &nlsP->max_radius,NULL);
789: PetscOptionsReal("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", nlsP->epsilon, &nlsP->epsilon,NULL);
790: PetscOptionsTail();
791: TaoLineSearchSetFromOptions(tao->linesearch);
792: KSPSetFromOptions(tao->ksp);
793: return 0;
794: }
796: /*------------------------------------------------------------*/
797: static PetscErrorCode TaoView_NLS(Tao tao, PetscViewer viewer)
798: {
799: TAO_NLS *nlsP = (TAO_NLS *)tao->data;
800: PetscBool isascii;
802: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
803: if (isascii) {
804: PetscViewerASCIIPushTab(viewer);
805: PetscViewerASCIIPrintf(viewer, "Newton steps: %D\n", nlsP->newt);
806: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", nlsP->bfgs);
807: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", nlsP->grad);
809: PetscViewerASCIIPrintf(viewer, "nls ksp atol: %D\n", nlsP->ksp_atol);
810: PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %D\n", nlsP->ksp_rtol);
811: PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %D\n", nlsP->ksp_ctol);
812: PetscViewerASCIIPrintf(viewer, "nls ksp negc: %D\n", nlsP->ksp_negc);
813: PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %D\n", nlsP->ksp_dtol);
814: PetscViewerASCIIPrintf(viewer, "nls ksp iter: %D\n", nlsP->ksp_iter);
815: PetscViewerASCIIPrintf(viewer, "nls ksp othr: %D\n", nlsP->ksp_othr);
816: PetscViewerASCIIPopTab(viewer);
817: }
818: return 0;
819: }
821: /* ---------------------------------------------------------- */
822: /*MC
823: TAONLS - Newton's method with linesearch for unconstrained minimization.
824: At each iteration, the Newton line search method solves the symmetric
825: system of equations to obtain the step diretion dk:
826: Hk dk = -gk
827: a More-Thuente line search is applied on the direction dk to approximately
828: solve
829: min_t f(xk + t d_k)
831: Options Database Keys:
832: + -tao_nls_init_type - "constant","direction","interpolation"
833: . -tao_nls_update_type - "step","direction","interpolation"
834: . -tao_nls_sval - perturbation starting value
835: . -tao_nls_imin - minimum initial perturbation
836: . -tao_nls_imax - maximum initial perturbation
837: . -tao_nls_pmin - minimum perturbation
838: . -tao_nls_pmax - maximum perturbation
839: . -tao_nls_pgfac - growth factor
840: . -tao_nls_psfac - shrink factor
841: . -tao_nls_imfac - initial merit factor
842: . -tao_nls_pmgfac - merit growth factor
843: . -tao_nls_pmsfac - merit shrink factor
844: . -tao_nls_eta1 - poor steplength; reduce radius
845: . -tao_nls_eta2 - reasonable steplength; leave radius
846: . -tao_nls_eta3 - good steplength; increase readius
847: . -tao_nls_eta4 - excellent steplength; greatly increase radius
848: . -tao_nls_alpha1 - alpha1 reduction
849: . -tao_nls_alpha2 - alpha2 reduction
850: . -tao_nls_alpha3 - alpha3 reduction
851: . -tao_nls_alpha4 - alpha4 reduction
852: . -tao_nls_alpha - alpha5 reduction
853: . -tao_nls_mu1 - mu1 interpolation update
854: . -tao_nls_mu2 - mu2 interpolation update
855: . -tao_nls_gamma1 - gamma1 interpolation update
856: . -tao_nls_gamma2 - gamma2 interpolation update
857: . -tao_nls_gamma3 - gamma3 interpolation update
858: . -tao_nls_gamma4 - gamma4 interpolation update
859: . -tao_nls_theta - theta interpolation update
860: . -tao_nls_omega1 - omega1 step update
861: . -tao_nls_omega2 - omega2 step update
862: . -tao_nls_omega3 - omega3 step update
863: . -tao_nls_omega4 - omega4 step update
864: . -tao_nls_omega5 - omega5 step update
865: . -tao_nls_mu1_i - mu1 interpolation init factor
866: . -tao_nls_mu2_i - mu2 interpolation init factor
867: . -tao_nls_gamma1_i - gamma1 interpolation init factor
868: . -tao_nls_gamma2_i - gamma2 interpolation init factor
869: . -tao_nls_gamma3_i - gamma3 interpolation init factor
870: . -tao_nls_gamma4_i - gamma4 interpolation init factor
871: - -tao_nls_theta_i - theta interpolation init factor
873: Level: beginner
874: M*/
876: PETSC_EXTERN PetscErrorCode TaoCreate_NLS(Tao tao)
877: {
878: TAO_NLS *nlsP;
879: const char *morethuente_type = TAOLINESEARCHMT;
881: PetscNewLog(tao,&nlsP);
883: tao->ops->setup = TaoSetUp_NLS;
884: tao->ops->solve = TaoSolve_NLS;
885: tao->ops->view = TaoView_NLS;
886: tao->ops->setfromoptions = TaoSetFromOptions_NLS;
887: tao->ops->destroy = TaoDestroy_NLS;
889: /* Override default settings (unless already changed) */
890: if (!tao->max_it_changed) tao->max_it = 50;
891: if (!tao->trust0_changed) tao->trust0 = 100.0;
893: tao->data = (void*)nlsP;
895: nlsP->sval = 0.0;
896: nlsP->imin = 1.0e-4;
897: nlsP->imax = 1.0e+2;
898: nlsP->imfac = 1.0e-1;
900: nlsP->pmin = 1.0e-12;
901: nlsP->pmax = 1.0e+2;
902: nlsP->pgfac = 1.0e+1;
903: nlsP->psfac = 4.0e-1;
904: nlsP->pmgfac = 1.0e-1;
905: nlsP->pmsfac = 1.0e-1;
907: /* Default values for trust-region radius update based on steplength */
908: nlsP->nu1 = 0.25;
909: nlsP->nu2 = 0.50;
910: nlsP->nu3 = 1.00;
911: nlsP->nu4 = 1.25;
913: nlsP->omega1 = 0.25;
914: nlsP->omega2 = 0.50;
915: nlsP->omega3 = 1.00;
916: nlsP->omega4 = 2.00;
917: nlsP->omega5 = 4.00;
919: /* Default values for trust-region radius update based on reduction */
920: nlsP->eta1 = 1.0e-4;
921: nlsP->eta2 = 0.25;
922: nlsP->eta3 = 0.50;
923: nlsP->eta4 = 0.90;
925: nlsP->alpha1 = 0.25;
926: nlsP->alpha2 = 0.50;
927: nlsP->alpha3 = 1.00;
928: nlsP->alpha4 = 2.00;
929: nlsP->alpha5 = 4.00;
931: /* Default values for trust-region radius update based on interpolation */
932: nlsP->mu1 = 0.10;
933: nlsP->mu2 = 0.50;
935: nlsP->gamma1 = 0.25;
936: nlsP->gamma2 = 0.50;
937: nlsP->gamma3 = 2.00;
938: nlsP->gamma4 = 4.00;
940: nlsP->theta = 0.05;
942: /* Default values for trust region initialization based on interpolation */
943: nlsP->mu1_i = 0.35;
944: nlsP->mu2_i = 0.50;
946: nlsP->gamma1_i = 0.0625;
947: nlsP->gamma2_i = 0.5;
948: nlsP->gamma3_i = 2.0;
949: nlsP->gamma4_i = 5.0;
951: nlsP->theta_i = 0.25;
953: /* Remaining parameters */
954: nlsP->min_radius = 1.0e-10;
955: nlsP->max_radius = 1.0e10;
956: nlsP->epsilon = 1.0e-6;
958: nlsP->init_type = NLS_INIT_INTERPOLATION;
959: nlsP->update_type = NLS_UPDATE_STEP;
961: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
962: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch,(PetscObject)tao,1);
963: TaoLineSearchSetType(tao->linesearch,morethuente_type);
964: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
965: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
967: /* Set linear solver to default for symmetric matrices */
968: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
969: PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1);
970: KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
971: KSPAppendOptionsPrefix(tao->ksp,"tao_nls_");
972: KSPSetType(tao->ksp,KSPSTCG);
973: return 0;
974: }