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