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:   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"));

 65:   /* Initialized variables */
 66:   pert = nlsP->sval;

 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;

 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));

 86:   PetscCall(KSPCGSetRadius(tao->ksp, nlsP->max_radius));

 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:   }

 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 Inf or NaN");

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);

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));

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;

129:     case NLS_INIT_INTERPOLATION:
130:       /* Use the initial radius specified */
131:       max_radius = 0.0;

133:       for (j = 0; j < j_max; ++j) {
134:         fmin  = f;
135:         sigma = 0.0;

137:         if (needH) {
138:           PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
139:           needH = 0;
140:         }

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:             }

154:             PetscCall(MatMult(tao->hessian, tao->gradient, nlsP->D));
155:             PetscCall(VecDot(tao->gradient, nlsP->D, &prered));

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:             }

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);

170:             if (PetscAbsScalar(kappa - (PetscReal)1.0) <= nlsP->mu1_i) {
171:               /* Great agreement */
172:               max_radius = PetscMax(max_radius, tao->trust);

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);

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:         }

216:         if (fmin < f) {
217:           f = fmin;
218:           PetscCall(VecAXPY(tao->solution, sigma, tao->gradient));
219:           PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));

221:           PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
222:           PetscCheck(!PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute gradient generated Inf or NaN");
223:           needH = 1;

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);

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;

238:     default:
239:       /* Norm of the first direction will initialize radius */
240:       tao->trust = 0.0;
241:       break;
242:     }
243:   }

245:   /* Set counter for gradient/reset steps*/
246:   nlsP->newt = 0;
247:   nlsP->bfgs = 0;
248:   nlsP->grad = 0;

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;

260:     /* Compute the Hessian */
261:     if (needH) PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));

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:     }

269:     if (nlsP->bfgs_pre) {
270:       PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient));
271:       ++bfgsUpdates;
272:     }

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));

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;

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;

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);

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));

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:     }

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:     }

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 Inf 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:       }

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));

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 */

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));

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;

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:       }

430:       ++nlsP->newt;
431:       stepType = NLS_NEWTON;
432:     }

434:     /* Perform the linesearch */
435:     fold = f;
436:     PetscCall(VecCopy(tao->solution, nlsP->Xold));
437:     PetscCall(VecCopy(tao->gradient, nlsP->Gold));

439:     PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason));
440:     PetscCall(TaoAddLineSearchCounts(tao));

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));

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:         }

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->solution, 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));

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;

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));

506:         bfgsUpdates = 1;
507:         ++nlsP->grad;
508:         stepType = NLS_GRADIENT;
509:         break;
510:       }
511:       PetscCall(VecScale(nlsP->D, -1.0));

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:     }

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:     }

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;

559:       case NLS_UPDATE_REDUCTION:
560:         if (stepType == NLS_NEWTON) {
561:           /*  Get predicted reduction */

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:               }

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;

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:               }

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);

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 */

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:       }

682:       /*  The radius may have been increased; modify if it is too large */
683:       tao->trust = PetscMin(tao->trust, nlsP->max_radius);
684:     }

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: }

697: /* ---------------------------------------------------------- */
698: static PetscErrorCode TaoSetUp_NLS(Tao tao)
699: {
700:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

702:   PetscFunctionBegin;
703:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
704:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
705:   if (!nlsP->W) PetscCall(VecDuplicate(tao->solution, &nlsP->W));
706:   if (!nlsP->D) PetscCall(VecDuplicate(tao->solution, &nlsP->D));
707:   if (!nlsP->Xold) PetscCall(VecDuplicate(tao->solution, &nlsP->Xold));
708:   if (!nlsP->Gold) PetscCall(VecDuplicate(tao->solution, &nlsP->Gold));
709:   nlsP->M        = NULL;
710:   nlsP->bfgs_pre = NULL;
711:   PetscFunctionReturn(PETSC_SUCCESS);
712: }

714: /*------------------------------------------------------------*/
715: static PetscErrorCode TaoDestroy_NLS(Tao tao)
716: {
717:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

719:   PetscFunctionBegin;
720:   if (tao->setupcalled) {
721:     PetscCall(VecDestroy(&nlsP->D));
722:     PetscCall(VecDestroy(&nlsP->W));
723:     PetscCall(VecDestroy(&nlsP->Xold));
724:     PetscCall(VecDestroy(&nlsP->Gold));
725:   }
726:   PetscCall(KSPDestroy(&tao->ksp));
727:   PetscCall(PetscFree(tao->data));
728:   PetscFunctionReturn(PETSC_SUCCESS);
729: }

731: /*------------------------------------------------------------*/
732: static PetscErrorCode TaoSetFromOptions_NLS(Tao tao, PetscOptionItems *PetscOptionsObject)
733: {
734:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

736:   PetscFunctionBegin;
737:   PetscOptionsHeadBegin(PetscOptionsObject, "Newton line search method for unconstrained optimization");
738:   PetscCall(PetscOptionsEList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[nlsP->init_type], &nlsP->init_type, NULL));
739:   PetscCall(PetscOptionsEList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[nlsP->update_type], &nlsP->update_type, NULL));
740:   PetscCall(PetscOptionsReal("-tao_nls_sval", "perturbation starting value", "", nlsP->sval, &nlsP->sval, NULL));
741:   PetscCall(PetscOptionsReal("-tao_nls_imin", "minimum initial perturbation", "", nlsP->imin, &nlsP->imin, NULL));
742:   PetscCall(PetscOptionsReal("-tao_nls_imax", "maximum initial perturbation", "", nlsP->imax, &nlsP->imax, NULL));
743:   PetscCall(PetscOptionsReal("-tao_nls_imfac", "initial merit factor", "", nlsP->imfac, &nlsP->imfac, NULL));
744:   PetscCall(PetscOptionsReal("-tao_nls_pmin", "minimum perturbation", "", nlsP->pmin, &nlsP->pmin, NULL));
745:   PetscCall(PetscOptionsReal("-tao_nls_pmax", "maximum perturbation", "", nlsP->pmax, &nlsP->pmax, NULL));
746:   PetscCall(PetscOptionsReal("-tao_nls_pgfac", "growth factor", "", nlsP->pgfac, &nlsP->pgfac, NULL));
747:   PetscCall(PetscOptionsReal("-tao_nls_psfac", "shrink factor", "", nlsP->psfac, &nlsP->psfac, NULL));
748:   PetscCall(PetscOptionsReal("-tao_nls_pmgfac", "merit growth factor", "", nlsP->pmgfac, &nlsP->pmgfac, NULL));
749:   PetscCall(PetscOptionsReal("-tao_nls_pmsfac", "merit shrink factor", "", nlsP->pmsfac, &nlsP->pmsfac, NULL));
750:   PetscCall(PetscOptionsReal("-tao_nls_eta1", "poor steplength; reduce radius", "", nlsP->eta1, &nlsP->eta1, NULL));
751:   PetscCall(PetscOptionsReal("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", nlsP->eta2, &nlsP->eta2, NULL));
752:   PetscCall(PetscOptionsReal("-tao_nls_eta3", "good steplength; increase radius", "", nlsP->eta3, &nlsP->eta3, NULL));
753:   PetscCall(PetscOptionsReal("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", nlsP->eta4, &nlsP->eta4, NULL));
754:   PetscCall(PetscOptionsReal("-tao_nls_alpha1", "", "", nlsP->alpha1, &nlsP->alpha1, NULL));
755:   PetscCall(PetscOptionsReal("-tao_nls_alpha2", "", "", nlsP->alpha2, &nlsP->alpha2, NULL));
756:   PetscCall(PetscOptionsReal("-tao_nls_alpha3", "", "", nlsP->alpha3, &nlsP->alpha3, NULL));
757:   PetscCall(PetscOptionsReal("-tao_nls_alpha4", "", "", nlsP->alpha4, &nlsP->alpha4, NULL));
758:   PetscCall(PetscOptionsReal("-tao_nls_alpha5", "", "", nlsP->alpha5, &nlsP->alpha5, NULL));
759:   PetscCall(PetscOptionsReal("-tao_nls_nu1", "poor steplength; reduce radius", "", nlsP->nu1, &nlsP->nu1, NULL));
760:   PetscCall(PetscOptionsReal("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", nlsP->nu2, &nlsP->nu2, NULL));
761:   PetscCall(PetscOptionsReal("-tao_nls_nu3", "good steplength; increase radius", "", nlsP->nu3, &nlsP->nu3, NULL));
762:   PetscCall(PetscOptionsReal("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", nlsP->nu4, &nlsP->nu4, NULL));
763:   PetscCall(PetscOptionsReal("-tao_nls_omega1", "", "", nlsP->omega1, &nlsP->omega1, NULL));
764:   PetscCall(PetscOptionsReal("-tao_nls_omega2", "", "", nlsP->omega2, &nlsP->omega2, NULL));
765:   PetscCall(PetscOptionsReal("-tao_nls_omega3", "", "", nlsP->omega3, &nlsP->omega3, NULL));
766:   PetscCall(PetscOptionsReal("-tao_nls_omega4", "", "", nlsP->omega4, &nlsP->omega4, NULL));
767:   PetscCall(PetscOptionsReal("-tao_nls_omega5", "", "", nlsP->omega5, &nlsP->omega5, NULL));
768:   PetscCall(PetscOptionsReal("-tao_nls_mu1_i", "", "", nlsP->mu1_i, &nlsP->mu1_i, NULL));
769:   PetscCall(PetscOptionsReal("-tao_nls_mu2_i", "", "", nlsP->mu2_i, &nlsP->mu2_i, NULL));
770:   PetscCall(PetscOptionsReal("-tao_nls_gamma1_i", "", "", nlsP->gamma1_i, &nlsP->gamma1_i, NULL));
771:   PetscCall(PetscOptionsReal("-tao_nls_gamma2_i", "", "", nlsP->gamma2_i, &nlsP->gamma2_i, NULL));
772:   PetscCall(PetscOptionsReal("-tao_nls_gamma3_i", "", "", nlsP->gamma3_i, &nlsP->gamma3_i, NULL));
773:   PetscCall(PetscOptionsReal("-tao_nls_gamma4_i", "", "", nlsP->gamma4_i, &nlsP->gamma4_i, NULL));
774:   PetscCall(PetscOptionsReal("-tao_nls_theta_i", "", "", nlsP->theta_i, &nlsP->theta_i, NULL));
775:   PetscCall(PetscOptionsReal("-tao_nls_mu1", "", "", nlsP->mu1, &nlsP->mu1, NULL));
776:   PetscCall(PetscOptionsReal("-tao_nls_mu2", "", "", nlsP->mu2, &nlsP->mu2, NULL));
777:   PetscCall(PetscOptionsReal("-tao_nls_gamma1", "", "", nlsP->gamma1, &nlsP->gamma1, NULL));
778:   PetscCall(PetscOptionsReal("-tao_nls_gamma2", "", "", nlsP->gamma2, &nlsP->gamma2, NULL));
779:   PetscCall(PetscOptionsReal("-tao_nls_gamma3", "", "", nlsP->gamma3, &nlsP->gamma3, NULL));
780:   PetscCall(PetscOptionsReal("-tao_nls_gamma4", "", "", nlsP->gamma4, &nlsP->gamma4, NULL));
781:   PetscCall(PetscOptionsReal("-tao_nls_theta", "", "", nlsP->theta, &nlsP->theta, NULL));
782:   PetscCall(PetscOptionsReal("-tao_nls_min_radius", "lower bound on initial radius", "", nlsP->min_radius, &nlsP->min_radius, NULL));
783:   PetscCall(PetscOptionsReal("-tao_nls_max_radius", "upper bound on radius", "", nlsP->max_radius, &nlsP->max_radius, NULL));
784:   PetscCall(PetscOptionsReal("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", nlsP->epsilon, &nlsP->epsilon, NULL));
785:   PetscOptionsHeadEnd();
786:   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
787:   PetscCall(KSPSetFromOptions(tao->ksp));
788:   PetscFunctionReturn(PETSC_SUCCESS);
789: }

791: /*------------------------------------------------------------*/
792: static PetscErrorCode TaoView_NLS(Tao tao, PetscViewer viewer)
793: {
794:   TAO_NLS  *nlsP = (TAO_NLS *)tao->data;
795:   PetscBool isascii;

797:   PetscFunctionBegin;
798:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
799:   if (isascii) {
800:     PetscCall(PetscViewerASCIIPushTab(viewer));
801:     PetscCall(PetscViewerASCIIPrintf(viewer, "Newton steps: %" PetscInt_FMT "\n", nlsP->newt));
802:     PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS steps: %" PetscInt_FMT "\n", nlsP->bfgs));
803:     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", nlsP->grad));

805:     PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp atol: %" PetscInt_FMT "\n", nlsP->ksp_atol));
806:     PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %" PetscInt_FMT "\n", nlsP->ksp_rtol));
807:     PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %" PetscInt_FMT "\n", nlsP->ksp_ctol));
808:     PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp negc: %" PetscInt_FMT "\n", nlsP->ksp_negc));
809:     PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %" PetscInt_FMT "\n", nlsP->ksp_dtol));
810:     PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp iter: %" PetscInt_FMT "\n", nlsP->ksp_iter));
811:     PetscCall(PetscViewerASCIIPrintf(viewer, "nls ksp othr: %" PetscInt_FMT "\n", nlsP->ksp_othr));
812:     PetscCall(PetscViewerASCIIPopTab(viewer));
813:   }
814:   PetscFunctionReturn(PETSC_SUCCESS);
815: }

817: /* ---------------------------------------------------------- */
818: /*MC
819:   TAONLS - Newton's method with linesearch for unconstrained minimization.
820:   At each iteration, the Newton line search method solves the symmetric
821:   system of equations to obtain the step direction dk:
822:               Hk dk = -gk
823:   a More-Thuente line search is applied on the direction dk to approximately
824:   solve
825:               min_t f(xk + t d_k)

827:     Options Database Keys:
828: + -tao_nls_init_type - "constant","direction","interpolation"
829: . -tao_nls_update_type - "step","direction","interpolation"
830: . -tao_nls_sval - perturbation starting value
831: . -tao_nls_imin - minimum initial perturbation
832: . -tao_nls_imax - maximum initial perturbation
833: . -tao_nls_pmin - minimum perturbation
834: . -tao_nls_pmax - maximum perturbation
835: . -tao_nls_pgfac - growth factor
836: . -tao_nls_psfac - shrink factor
837: . -tao_nls_imfac - initial merit factor
838: . -tao_nls_pmgfac - merit growth factor
839: . -tao_nls_pmsfac - merit shrink factor
840: . -tao_nls_eta1 - poor steplength; reduce radius
841: . -tao_nls_eta2 - reasonable steplength; leave radius
842: . -tao_nls_eta3 - good steplength; increase radius
843: . -tao_nls_eta4 - excellent steplength; greatly increase radius
844: . -tao_nls_alpha1 - alpha1 reduction
845: . -tao_nls_alpha2 - alpha2 reduction
846: . -tao_nls_alpha3 - alpha3 reduction
847: . -tao_nls_alpha4 - alpha4 reduction
848: . -tao_nls_alpha - alpha5 reduction
849: . -tao_nls_mu1 - mu1 interpolation update
850: . -tao_nls_mu2 - mu2 interpolation update
851: . -tao_nls_gamma1 - gamma1 interpolation update
852: . -tao_nls_gamma2 - gamma2 interpolation update
853: . -tao_nls_gamma3 - gamma3 interpolation update
854: . -tao_nls_gamma4 - gamma4 interpolation update
855: . -tao_nls_theta - theta interpolation update
856: . -tao_nls_omega1 - omega1 step update
857: . -tao_nls_omega2 - omega2 step update
858: . -tao_nls_omega3 - omega3 step update
859: . -tao_nls_omega4 - omega4 step update
860: . -tao_nls_omega5 - omega5 step update
861: . -tao_nls_mu1_i -  mu1 interpolation init factor
862: . -tao_nls_mu2_i -  mu2 interpolation init factor
863: . -tao_nls_gamma1_i -  gamma1 interpolation init factor
864: . -tao_nls_gamma2_i -  gamma2 interpolation init factor
865: . -tao_nls_gamma3_i -  gamma3 interpolation init factor
866: . -tao_nls_gamma4_i -  gamma4 interpolation init factor
867: - -tao_nls_theta_i -  theta interpolation init factor

869:   Level: beginner
870: M*/

872: PETSC_EXTERN PetscErrorCode TaoCreate_NLS(Tao tao)
873: {
874:   TAO_NLS    *nlsP;
875:   const char *morethuente_type = TAOLINESEARCHMT;

877:   PetscFunctionBegin;
878:   PetscCall(PetscNew(&nlsP));

880:   tao->ops->setup          = TaoSetUp_NLS;
881:   tao->ops->solve          = TaoSolve_NLS;
882:   tao->ops->view           = TaoView_NLS;
883:   tao->ops->setfromoptions = TaoSetFromOptions_NLS;
884:   tao->ops->destroy        = TaoDestroy_NLS;

886:   /* Override default settings (unless already changed) */
887:   if (!tao->max_it_changed) tao->max_it = 50;
888:   if (!tao->trust0_changed) tao->trust0 = 100.0;

890:   tao->data = (void *)nlsP;

892:   nlsP->sval  = 0.0;
893:   nlsP->imin  = 1.0e-4;
894:   nlsP->imax  = 1.0e+2;
895:   nlsP->imfac = 1.0e-1;

897:   nlsP->pmin   = 1.0e-12;
898:   nlsP->pmax   = 1.0e+2;
899:   nlsP->pgfac  = 1.0e+1;
900:   nlsP->psfac  = 4.0e-1;
901:   nlsP->pmgfac = 1.0e-1;
902:   nlsP->pmsfac = 1.0e-1;

904:   /*  Default values for trust-region radius update based on steplength */
905:   nlsP->nu1 = 0.25;
906:   nlsP->nu2 = 0.50;
907:   nlsP->nu3 = 1.00;
908:   nlsP->nu4 = 1.25;

910:   nlsP->omega1 = 0.25;
911:   nlsP->omega2 = 0.50;
912:   nlsP->omega3 = 1.00;
913:   nlsP->omega4 = 2.00;
914:   nlsP->omega5 = 4.00;

916:   /*  Default values for trust-region radius update based on reduction */
917:   nlsP->eta1 = 1.0e-4;
918:   nlsP->eta2 = 0.25;
919:   nlsP->eta3 = 0.50;
920:   nlsP->eta4 = 0.90;

922:   nlsP->alpha1 = 0.25;
923:   nlsP->alpha2 = 0.50;
924:   nlsP->alpha3 = 1.00;
925:   nlsP->alpha4 = 2.00;
926:   nlsP->alpha5 = 4.00;

928:   /*  Default values for trust-region radius update based on interpolation */
929:   nlsP->mu1 = 0.10;
930:   nlsP->mu2 = 0.50;

932:   nlsP->gamma1 = 0.25;
933:   nlsP->gamma2 = 0.50;
934:   nlsP->gamma3 = 2.00;
935:   nlsP->gamma4 = 4.00;

937:   nlsP->theta = 0.05;

939:   /*  Default values for trust region initialization based on interpolation */
940:   nlsP->mu1_i = 0.35;
941:   nlsP->mu2_i = 0.50;

943:   nlsP->gamma1_i = 0.0625;
944:   nlsP->gamma2_i = 0.5;
945:   nlsP->gamma3_i = 2.0;
946:   nlsP->gamma4_i = 5.0;

948:   nlsP->theta_i = 0.25;

950:   /*  Remaining parameters */
951:   nlsP->min_radius = 1.0e-10;
952:   nlsP->max_radius = 1.0e10;
953:   nlsP->epsilon    = 1.0e-6;

955:   nlsP->init_type   = NLS_INIT_INTERPOLATION;
956:   nlsP->update_type = NLS_UPDATE_STEP;

958:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
959:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
960:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
961:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
962:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));

964:   /*  Set linear solver to default for symmetric matrices */
965:   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
966:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
967:   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
968:   PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_nls_"));
969:   PetscCall(KSPSetType(tao->ksp, KSPSTCG));
970:   PetscFunctionReturn(PETSC_SUCCESS);
971: }