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