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:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
254:     ++tao->niter;
255:     tao->ksp_its = 0;

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

260:     /* Shift the Hessian matrix */
261:     if (pert > 0) {
262:       PetscCall(MatShift(tao->hessian, pert));
263:       if (tao->hessian != tao->hessian_pre) PetscCall(MatShift(tao->hessian_pre, pert));
264:     }

266:     if (nlsP->bfgs_pre) {
267:       PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient));
268:       ++bfgsUpdates;
269:     }

271:     /* Solve the Newton system of equations */
272:     PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));
273:     if (is_nash || is_stcg || is_gltr) {
274:       PetscCall(KSPCGSetRadius(tao->ksp, nlsP->max_radius));
275:       PetscCall(KSPSolve(tao->ksp, tao->gradient, nlsP->D));
276:       PetscCall(KSPGetIterationNumber(tao->ksp, &kspits));
277:       tao->ksp_its += kspits;
278:       tao->ksp_tot_its += kspits;
279:       PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));

281:       if (0.0 == tao->trust) {
282:         /* Radius was uninitialized; use the norm of the direction */
283:         if (norm_d > 0.0) {
284:           tao->trust = norm_d;

286:           /* Modify the radius if it is too large or small */
287:           tao->trust = PetscMax(tao->trust, nlsP->min_radius);
288:           tao->trust = PetscMin(tao->trust, nlsP->max_radius);
289:         } else {
290:           /* The direction was bad; set radius to default value and re-solve
291:              the trust-region subproblem to get a direction */
292:           tao->trust = tao->trust0;

294:           /* Modify the radius if it is too large or small */
295:           tao->trust = PetscMax(tao->trust, nlsP->min_radius);
296:           tao->trust = PetscMin(tao->trust, nlsP->max_radius);

298:           PetscCall(KSPCGSetRadius(tao->ksp, nlsP->max_radius));
299:           PetscCall(KSPSolve(tao->ksp, tao->gradient, nlsP->D));
300:           PetscCall(KSPGetIterationNumber(tao->ksp, &kspits));
301:           tao->ksp_its += kspits;
302:           tao->ksp_tot_its += kspits;
303:           PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));

305:           PetscCheck(norm_d != 0.0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Initial direction zero");
306:         }
307:       }
308:     } else {
309:       PetscCall(KSPSolve(tao->ksp, tao->gradient, nlsP->D));
310:       PetscCall(KSPGetIterationNumber(tao->ksp, &kspits));
311:       tao->ksp_its += kspits;
312:       tao->ksp_tot_its += kspits;
313:     }
314:     PetscCall(VecScale(nlsP->D, -1.0));
315:     PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason));
316:     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (nlsP->bfgs_pre)) {
317:       /* Preconditioner is numerically indefinite; reset the
318:          approximate if using BFGS preconditioning. */
319:       PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE));
320:       PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient));
321:       bfgsUpdates = 1;
322:     }

324:     if (KSP_CONVERGED_ATOL == ksp_reason) {
325:       ++nlsP->ksp_atol;
326:     } else if (KSP_CONVERGED_RTOL == ksp_reason) {
327:       ++nlsP->ksp_rtol;
328:     } else if (KSP_CONVERGED_STEP_LENGTH == ksp_reason) {
329:       ++nlsP->ksp_ctol;
330:     } else if (KSP_CONVERGED_NEG_CURVE == ksp_reason) {
331:       ++nlsP->ksp_negc;
332:     } else if (KSP_DIVERGED_DTOL == ksp_reason) {
333:       ++nlsP->ksp_dtol;
334:     } else if (KSP_DIVERGED_ITS == ksp_reason) {
335:       ++nlsP->ksp_iter;
336:     } else {
337:       ++nlsP->ksp_othr;
338:     }

340:     /* Check for success (descent direction) */
341:     PetscCall(VecDot(nlsP->D, tao->gradient, &gdx));
342:     if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
343:       /* Newton step is not descent or direction produced Inf or NaN
344:          Update the perturbation for next time */
345:       if (pert <= 0.0) {
346:         /* Initialize the perturbation */
347:         pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
348:         if (is_gltr) {
349:           PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min));
350:           pert = PetscMax(pert, -e_min);
351:         }
352:       } else {
353:         /* Increase the perturbation */
354:         pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
355:       }

357:       if (!nlsP->bfgs_pre) {
358:         /* We don't have the bfgs matrix around and updated
359:            Must use gradient direction in this case */
360:         PetscCall(VecCopy(tao->gradient, nlsP->D));
361:         PetscCall(VecScale(nlsP->D, -1.0));
362:         ++nlsP->grad;
363:         stepType = NLS_GRADIENT;
364:       } else {
365:         /* Attempt to use the BFGS direction */
366:         PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D));
367:         PetscCall(VecScale(nlsP->D, -1.0));

369:         /* Check for success (descent direction) */
370:         PetscCall(VecDot(tao->gradient, nlsP->D, &gdx));
371:         if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
372:           /* BFGS direction is not descent or direction produced not a number
373:              We can assert bfgsUpdates > 1 in this case because
374:              the first solve produces the scaled gradient direction,
375:              which is guaranteed to be descent */

377:           /* Use steepest descent direction (scaled) */
378:           PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE));
379:           PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient));
380:           PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D));
381:           PetscCall(VecScale(nlsP->D, -1.0));

383:           bfgsUpdates = 1;
384:           ++nlsP->grad;
385:           stepType = NLS_GRADIENT;
386:         } else {
387:           PetscCall(MatLMVMGetUpdateCount(nlsP->M, &bfgsUpdates));
388:           if (1 == bfgsUpdates) {
389:             /* The first BFGS direction is always the scaled gradient */
390:             ++nlsP->grad;
391:             stepType = NLS_GRADIENT;
392:           } else {
393:             ++nlsP->bfgs;
394:             stepType = NLS_BFGS;
395:           }
396:         }
397:       }
398:     } else {
399:       /* Computed Newton step is descent */
400:       switch (ksp_reason) {
401:       case KSP_DIVERGED_NANORINF:
402:       case KSP_DIVERGED_BREAKDOWN:
403:       case KSP_DIVERGED_INDEFINITE_MAT:
404:       case KSP_DIVERGED_INDEFINITE_PC:
405:       case KSP_CONVERGED_NEG_CURVE:
406:         /* Matrix or preconditioner is indefinite; increase perturbation */
407:         if (pert <= 0.0) {
408:           /* Initialize the perturbation */
409:           pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
410:           if (is_gltr) {
411:             PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min));
412:             pert = PetscMax(pert, -e_min);
413:           }
414:         } else {
415:           /* Increase the perturbation */
416:           pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
417:         }
418:         break;

420:       default:
421:         /* Newton step computation is good; decrease perturbation */
422:         pert = PetscMin(nlsP->psfac * pert, nlsP->pmsfac * gnorm);
423:         if (pert < nlsP->pmin) pert = 0.0;
424:         break;
425:       }

427:       ++nlsP->newt;
428:       stepType = NLS_NEWTON;
429:     }

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

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

439:     while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NLS_GRADIENT) { /* Linesearch failed */
440:       f = fold;
441:       PetscCall(VecCopy(nlsP->Xold, tao->solution));
442:       PetscCall(VecCopy(nlsP->Gold, tao->gradient));

444:       switch (stepType) {
445:       case NLS_NEWTON:
446:         /* Failed to obtain acceptable iterate with Newton 1step
447:            Update the perturbation for next time */
448:         if (pert <= 0.0) {
449:           /* Initialize the perturbation */
450:           pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
451:           if (is_gltr) {
452:             PetscCall(KSPGLTRGetMinEig(tao->ksp, &e_min));
453:             pert = PetscMax(pert, -e_min);
454:           }
455:         } else {
456:           /* Increase the perturbation */
457:           pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
458:         }

460:         if (!nlsP->bfgs_pre) {
461:           /* We don't have the bfgs matrix around and being updated
462:              Must use gradient direction in this case */
463:           PetscCall(VecCopy(tao->gradient, nlsP->D));
464:           ++nlsP->grad;
465:           stepType = NLS_GRADIENT;
466:         } else {
467:           /* Attempt to use the BFGS direction */
468:           PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D));
469:           /* Check for success (descent direction) */
470:           PetscCall(VecDot(tao->solution, nlsP->D, &gdx));
471:           if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
472:             /* BFGS direction is not descent or direction produced not a number
473:                We can assert bfgsUpdates > 1 in this case
474:                Use steepest descent direction (scaled) */
475:             PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE));
476:             PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient));
477:             PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D));

479:             bfgsUpdates = 1;
480:             ++nlsP->grad;
481:             stepType = NLS_GRADIENT;
482:           } else {
483:             if (1 == bfgsUpdates) {
484:               /* The first BFGS direction is always the scaled gradient */
485:               ++nlsP->grad;
486:               stepType = NLS_GRADIENT;
487:             } else {
488:               ++nlsP->bfgs;
489:               stepType = NLS_BFGS;
490:             }
491:           }
492:         }
493:         break;

495:       case NLS_BFGS:
496:         /* Can only enter if pc_type == NLS_PC_BFGS
497:            Failed to obtain acceptable iterate with BFGS step
498:            Attempt to use the scaled gradient direction */
499:         PetscCall(MatLMVMReset(nlsP->M, PETSC_FALSE));
500:         PetscCall(MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient));
501:         PetscCall(MatSolve(nlsP->M, tao->gradient, nlsP->D));

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

510:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason));
511:       PetscCall(TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full));
512:       PetscCall(TaoAddLineSearchCounts(tao));
513:     }

515:     if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
516:       /* Failed to find an improving point */
517:       f = fold;
518:       PetscCall(VecCopy(nlsP->Xold, tao->solution));
519:       PetscCall(VecCopy(nlsP->Gold, tao->gradient));
520:       step        = 0.0;
521:       tao->reason = TAO_DIVERGED_LS_FAILURE;
522:       break;
523:     }

525:     /* Update trust region radius */
526:     if (is_nash || is_stcg || is_gltr) {
527:       switch (nlsP->update_type) {
528:       case NLS_UPDATE_STEP:
529:         if (stepType == NLS_NEWTON) {
530:           if (step < nlsP->nu1) {
531:             /* Very bad step taken; reduce radius */
532:             tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
533:           } else if (step < nlsP->nu2) {
534:             /* Reasonably bad step taken; reduce radius */
535:             tao->trust = nlsP->omega2 * PetscMin(norm_d, tao->trust);
536:           } else if (step < nlsP->nu3) {
537:             /*  Reasonable step was taken; leave radius alone */
538:             if (nlsP->omega3 < 1.0) {
539:               tao->trust = nlsP->omega3 * PetscMin(norm_d, tao->trust);
540:             } else if (nlsP->omega3 > 1.0) {
541:               tao->trust = PetscMax(nlsP->omega3 * norm_d, tao->trust);
542:             }
543:           } else if (step < nlsP->nu4) {
544:             /*  Full step taken; increase the radius */
545:             tao->trust = PetscMax(nlsP->omega4 * norm_d, tao->trust);
546:           } else {
547:             /*  More than full step taken; increase the radius */
548:             tao->trust = PetscMax(nlsP->omega5 * norm_d, tao->trust);
549:           }
550:         } else {
551:           /*  Newton step was not good; reduce the radius */
552:           tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
553:         }
554:         break;

556:       case NLS_UPDATE_REDUCTION:
557:         if (stepType == NLS_NEWTON) {
558:           /*  Get predicted reduction */

560:           PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
561:           if (prered >= 0.0) {
562:             /*  The predicted reduction has the wrong sign.  This cannot */
563:             /*  happen in infinite precision arithmetic.  Step should */
564:             /*  be rejected! */
565:             tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
566:           } else {
567:             if (PetscIsInfOrNanReal(f_full)) {
568:               tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
569:             } else {
570:               /*  Compute and actual reduction */
571:               actred = fold - f_full;
572:               prered = -prered;
573:               if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) {
574:                 kappa = 1.0;
575:               } else {
576:                 kappa = actred / prered;
577:               }

579:               /*  Accept of reject the step and update radius */
580:               if (kappa < nlsP->eta1) {
581:                 /*  Very bad step */
582:                 tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
583:               } else if (kappa < nlsP->eta2) {
584:                 /*  Marginal bad step */
585:                 tao->trust = nlsP->alpha2 * PetscMin(tao->trust, norm_d);
586:               } else if (kappa < nlsP->eta3) {
587:                 /*  Reasonable step */
588:                 if (nlsP->alpha3 < 1.0) {
589:                   tao->trust = nlsP->alpha3 * PetscMin(norm_d, tao->trust);
590:                 } else if (nlsP->alpha3 > 1.0) {
591:                   tao->trust = PetscMax(nlsP->alpha3 * norm_d, tao->trust);
592:                 }
593:               } else if (kappa < nlsP->eta4) {
594:                 /*  Good step */
595:                 tao->trust = PetscMax(nlsP->alpha4 * norm_d, tao->trust);
596:               } else {
597:                 /*  Very good step */
598:                 tao->trust = PetscMax(nlsP->alpha5 * norm_d, tao->trust);
599:               }
600:             }
601:           }
602:         } else {
603:           /*  Newton step was not good; reduce the radius */
604:           tao->trust = nlsP->alpha1 * PetscMin(norm_d, tao->trust);
605:         }
606:         break;

608:       default:
609:         if (stepType == NLS_NEWTON) {
610:           PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
611:           if (prered >= 0.0) {
612:             /*  The predicted reduction has the wrong sign.  This cannot */
613:             /*  happen in infinite precision arithmetic.  Step should */
614:             /*  be rejected! */
615:             tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
616:           } else {
617:             if (PetscIsInfOrNanReal(f_full)) {
618:               tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
619:             } else {
620:               actred = fold - f_full;
621:               prered = -prered;
622:               if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) {
623:                 kappa = 1.0;
624:               } else {
625:                 kappa = actred / prered;
626:               }

628:               tau_1   = nlsP->theta * gdx / (nlsP->theta * gdx - (1.0 - nlsP->theta) * prered + actred);
629:               tau_2   = nlsP->theta * gdx / (nlsP->theta * gdx + (1.0 + nlsP->theta) * prered - actred);
630:               tau_min = PetscMin(tau_1, tau_2);
631:               tau_max = PetscMax(tau_1, tau_2);

633:               if (kappa >= 1.0 - nlsP->mu1) {
634:                 /*  Great agreement */
635:                 if (tau_max < 1.0) {
636:                   tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
637:                 } else if (tau_max > nlsP->gamma4) {
638:                   tao->trust = PetscMax(tao->trust, nlsP->gamma4 * norm_d);
639:                 } else {
640:                   tao->trust = PetscMax(tao->trust, tau_max * norm_d);
641:                 }
642:               } else if (kappa >= 1.0 - nlsP->mu2) {
643:                 /*  Good agreement */

645:                 if (tau_max < nlsP->gamma2) {
646:                   tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
647:                 } else if (tau_max > nlsP->gamma3) {
648:                   tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
649:                 } else if (tau_max < 1.0) {
650:                   tao->trust = tau_max * PetscMin(tao->trust, norm_d);
651:                 } else {
652:                   tao->trust = PetscMax(tao->trust, tau_max * norm_d);
653:                 }
654:               } else {
655:                 /*  Not good agreement */
656:                 if (tau_min > 1.0) {
657:                   tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
658:                 } else if (tau_max < nlsP->gamma1) {
659:                   tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
660:                 } else if ((tau_min < nlsP->gamma1) && (tau_max >= 1.0)) {
661:                   tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
662:                 } else if ((tau_1 >= nlsP->gamma1) && (tau_1 < 1.0) && ((tau_2 < nlsP->gamma1) || (tau_2 >= 1.0))) {
663:                   tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
664:                 } else if ((tau_2 >= nlsP->gamma1) && (tau_2 < 1.0) && ((tau_1 < nlsP->gamma1) || (tau_2 >= 1.0))) {
665:                   tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
666:                 } else {
667:                   tao->trust = tau_max * PetscMin(tao->trust, norm_d);
668:                 }
669:               }
670:             }
671:           }
672:         } else {
673:           /*  Newton step was not good; reduce the radius */
674:           tao->trust = nlsP->gamma1 * PetscMin(norm_d, tao->trust);
675:         }
676:         break;
677:       }

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

683:     /*  Check for termination */
684:     PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
685:     PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
686:     needH = 1;
687:     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
688:     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
689:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
690:   }
691:   PetscFunctionReturn(PETSC_SUCCESS);
692: }

694: /* ---------------------------------------------------------- */
695: static PetscErrorCode TaoSetUp_NLS(Tao tao)
696: {
697:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

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

711: /*------------------------------------------------------------*/
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: /*------------------------------------------------------------*/
729: static PetscErrorCode TaoSetFromOptions_NLS(Tao tao, PetscOptionItems *PetscOptionsObject)
730: {
731:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

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

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

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

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

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

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

866:   Level: beginner
867: M*/

869: PETSC_EXTERN PetscErrorCode TaoCreate_NLS(Tao tao)
870: {
871:   TAO_NLS    *nlsP;
872:   const char *morethuente_type = TAOLINESEARCHMT;

874:   PetscFunctionBegin;
875:   PetscCall(PetscNew(&nlsP));

877:   tao->ops->setup          = TaoSetUp_NLS;
878:   tao->ops->solve          = TaoSolve_NLS;
879:   tao->ops->view           = TaoView_NLS;
880:   tao->ops->setfromoptions = TaoSetFromOptions_NLS;
881:   tao->ops->destroy        = TaoDestroy_NLS;

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

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

889:   nlsP->sval  = 0.0;
890:   nlsP->imin  = 1.0e-4;
891:   nlsP->imax  = 1.0e+2;
892:   nlsP->imfac = 1.0e-1;

894:   nlsP->pmin   = 1.0e-12;
895:   nlsP->pmax   = 1.0e+2;
896:   nlsP->pgfac  = 1.0e+1;
897:   nlsP->psfac  = 4.0e-1;
898:   nlsP->pmgfac = 1.0e-1;
899:   nlsP->pmsfac = 1.0e-1;

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

907:   nlsP->omega1 = 0.25;
908:   nlsP->omega2 = 0.50;
909:   nlsP->omega3 = 1.00;
910:   nlsP->omega4 = 2.00;
911:   nlsP->omega5 = 4.00;

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

919:   nlsP->alpha1 = 0.25;
920:   nlsP->alpha2 = 0.50;
921:   nlsP->alpha3 = 1.00;
922:   nlsP->alpha4 = 2.00;
923:   nlsP->alpha5 = 4.00;

925:   /*  Default values for trust-region radius update based on interpolation */
926:   nlsP->mu1 = 0.10;
927:   nlsP->mu2 = 0.50;

929:   nlsP->gamma1 = 0.25;
930:   nlsP->gamma2 = 0.50;
931:   nlsP->gamma3 = 2.00;
932:   nlsP->gamma4 = 4.00;

934:   nlsP->theta = 0.05;

936:   /*  Default values for trust region initialization based on interpolation */
937:   nlsP->mu1_i = 0.35;
938:   nlsP->mu2_i = 0.50;

940:   nlsP->gamma1_i = 0.0625;
941:   nlsP->gamma2_i = 0.5;
942:   nlsP->gamma3_i = 2.0;
943:   nlsP->gamma4_i = 5.0;

945:   nlsP->theta_i = 0.25;

947:   /*  Remaining parameters */
948:   nlsP->min_radius = 1.0e-10;
949:   nlsP->max_radius = 1.0e10;
950:   nlsP->epsilon    = 1.0e-6;

952:   nlsP->init_type   = NLS_INIT_INTERPOLATION;
953:   nlsP->update_type = NLS_UPDATE_STEP;

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

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