Actual source code: nls.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  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:   PetscErrorCode               ierr;
 41:   TAO_NLS                      *nlsP = (TAO_NLS *)tao->data;
 42:   KSPType                      ksp_type;
 43:   PetscBool                    is_nash,is_stcg,is_gltr,is_bfgs,is_jacobi,is_symmetric,sym_set;
 44:   KSPConvergedReason           ksp_reason;
 45:   PC                           pc;
 46:   TaoLineSearchConvergedReason ls_reason;

 48:   PetscReal                    fmin, ftrial, f_full, prered, actred, kappa, sigma;
 49:   PetscReal                    tau, tau_1, tau_2, tau_max, tau_min, max_radius;
 50:   PetscReal                    f, fold, gdx, gnorm, pert;
 51:   PetscReal                    step = 1.0;
 52:   PetscReal                    norm_d = 0.0, e_min;

 54:   PetscInt                     stepType;
 55:   PetscInt                     bfgsUpdates = 0;
 56:   PetscInt                     n,N,kspits;
 57:   PetscInt                     needH = 1;

 59:   PetscInt                     i_max = 5;
 60:   PetscInt                     j_max = 1;
 61:   PetscInt                     i, j;

 64:   if (tao->XL || tao->XU || tao->ops->computebounds) {
 65:     PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by nls algorithm\n");
 66:   }

 68:   /* Initialized variables */
 69:   pert = nlsP->sval;

 71:   /* Number of times ksp stopped because of these reasons */
 72:   nlsP->ksp_atol = 0;
 73:   nlsP->ksp_rtol = 0;
 74:   nlsP->ksp_dtol = 0;
 75:   nlsP->ksp_ctol = 0;
 76:   nlsP->ksp_negc = 0;
 77:   nlsP->ksp_iter = 0;
 78:   nlsP->ksp_othr = 0;

 80:   /* Initialize trust-region radius when using nash, stcg, or gltr
 81:      Command automatically ignored for other methods
 82:      Will be reset during the first iteration
 83:   */
 84:   KSPGetType(tao->ksp,&ksp_type);
 85:   PetscStrcmp(ksp_type,KSPNASH,&is_nash);
 86:   PetscStrcmp(ksp_type,KSPSTCG,&is_stcg);
 87:   PetscStrcmp(ksp_type,KSPGLTR,&is_gltr);

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

 91:   if (is_nash || is_stcg || is_gltr) {
 92:     if (tao->trust0 < 0.0) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_OUTOFRANGE,"Initial radius negative");
 93:     tao->trust = tao->trust0;
 94:     tao->trust = PetscMax(tao->trust, nlsP->min_radius);
 95:     tao->trust = PetscMin(tao->trust, nlsP->max_radius);
 96:   }

 98:   /* Check convergence criteria */
 99:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
100:   TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
101:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

103:   tao->reason = TAO_CONTINUE_ITERATING;
104:   TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
105:   TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
106:   (*tao->ops->convergencetest)(tao,tao->cnvP);
107:   if (tao->reason != TAO_CONTINUE_ITERATING) return(0);

109:   /* Allocate the vectors needed for the BFGS approximation */
110:   KSPGetPC(tao->ksp, &pc);
111:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);
112:   PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);
113:   if (is_bfgs) {
114:     nlsP->bfgs_pre = pc;
115:     PCLMVMGetMatLMVM(nlsP->bfgs_pre, &nlsP->M);
116:     VecGetLocalSize(tao->solution, &n);
117:     VecGetSize(tao->solution, &N);
118:     MatSetSizes(nlsP->M, n, n, N, N);
119:     MatLMVMAllocate(nlsP->M, tao->solution, tao->gradient);
120:     MatIsSymmetricKnown(nlsP->M, &sym_set, &is_symmetric);
121:     if (!sym_set || !is_symmetric) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
122:   } else if (is_jacobi) {
123:     PCJacobiSetUseAbs(pc,PETSC_TRUE);
124:   }

126:   /* Initialize trust-region radius.  The initialization is only performed
127:      when we are using Nash, Steihaug-Toint or the Generalized Lanczos method. */
128:   if (is_nash || is_stcg || is_gltr) {
129:     switch (nlsP->init_type) {
130:     case NLS_INIT_CONSTANT:
131:       /* Use the initial radius specified */
132:       break;

134:     case NLS_INIT_INTERPOLATION:
135:       /* Use the initial radius specified */
136:       max_radius = 0.0;

138:       for (j = 0; j < j_max; ++j) {
139:         fmin = f;
140:         sigma = 0.0;

142:         if (needH) {
143:           TaoComputeHessian(tao, tao->solution,tao->hessian,tao->hessian_pre);
144:           needH = 0;
145:         }

147:         for (i = 0; i < i_max; ++i) {
148:           VecCopy(tao->solution,nlsP->W);
149:           VecAXPY(nlsP->W,-tao->trust/gnorm,tao->gradient);
150:           TaoComputeObjective(tao, nlsP->W, &ftrial);
151:           if (PetscIsInfOrNanReal(ftrial)) {
152:             tau = nlsP->gamma1_i;
153:           } else {
154:             if (ftrial < fmin) {
155:               fmin = ftrial;
156:               sigma = -tao->trust / gnorm;
157:             }

159:             MatMult(tao->hessian, tao->gradient, nlsP->D);
160:             VecDot(tao->gradient, nlsP->D, &prered);

162:             prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
163:             actred = f - ftrial;
164:             if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) {
165:               kappa = 1.0;
166:             } else {
167:               kappa = actred / prered;
168:             }

170:             tau_1 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust + (1.0 - nlsP->theta_i) * prered - actred);
171:             tau_2 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust - (1.0 + nlsP->theta_i) * prered + actred);
172:             tau_min = PetscMin(tau_1, tau_2);
173:             tau_max = PetscMax(tau_1, tau_2);

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

179:               if (tau_max < 1.0) {
180:                 tau = nlsP->gamma3_i;
181:               } else if (tau_max > nlsP->gamma4_i) {
182:                 tau = nlsP->gamma4_i;
183:               } else if (tau_1 >= 1.0 && tau_1 <= nlsP->gamma4_i && tau_2 < 1.0) {
184:                 tau = tau_1;
185:               } else if (tau_2 >= 1.0 && tau_2 <= nlsP->gamma4_i && tau_1 < 1.0) {
186:                 tau = tau_2;
187:               } else {
188:                 tau = tau_max;
189:               }
190:             } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= nlsP->mu2_i) {
191:               /* Good agreement */
192:               max_radius = PetscMax(max_radius, tao->trust);

194:               if (tau_max < nlsP->gamma2_i) {
195:                 tau = nlsP->gamma2_i;
196:               } else if (tau_max > nlsP->gamma3_i) {
197:                 tau = nlsP->gamma3_i;
198:               } else {
199:                 tau = tau_max;
200:               }
201:             } else {
202:               /* Not good agreement */
203:               if (tau_min > 1.0) {
204:                 tau = nlsP->gamma2_i;
205:               } else if (tau_max < nlsP->gamma1_i) {
206:                 tau = nlsP->gamma1_i;
207:               } else if ((tau_min < nlsP->gamma1_i) && (tau_max >= 1.0)) {
208:                 tau = nlsP->gamma1_i;
209:               } else if ((tau_1 >= nlsP->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < nlsP->gamma1_i) || (tau_2 >= 1.0))) {
210:                 tau = tau_1;
211:               } else if ((tau_2 >= nlsP->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < nlsP->gamma1_i) || (tau_2 >= 1.0))) {
212:                 tau = tau_2;
213:               } else {
214:                 tau = tau_max;
215:               }
216:             }
217:           }
218:           tao->trust = tau * tao->trust;
219:         }

221:         if (fmin < f) {
222:           f = fmin;
223:           VecAXPY(tao->solution,sigma,tao->gradient);
224:           TaoComputeGradient(tao,tao->solution,tao->gradient);

226:           TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
227:           if (PetscIsInfOrNanReal(gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute gradient generated Inf or NaN");
228:           needH = 1;

230:           TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
231:           TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
232:           (*tao->ops->convergencetest)(tao,tao->cnvP);
233:           if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
234:         }
235:       }
236:       tao->trust = PetscMax(tao->trust, max_radius);

238:       /* Modify the radius if it is too large or small */
239:       tao->trust = PetscMax(tao->trust, nlsP->min_radius);
240:       tao->trust = PetscMin(tao->trust, nlsP->max_radius);
241:       break;

243:     default:
244:       /* Norm of the first direction will initialize radius */
245:       tao->trust = 0.0;
246:       break;
247:     }
248:   }

250:   /* Set counter for gradient/reset steps*/
251:   nlsP->newt = 0;
252:   nlsP->bfgs = 0;
253:   nlsP->grad = 0;

255:   /* Have not converged; continue with Newton method */
256:   while (tao->reason == TAO_CONTINUE_ITERATING) {
257:     /* Call general purpose update function */
258:     if (tao->ops->update) {
259:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
260:     }
261:     ++tao->niter;
262:     tao->ksp_its=0;

264:     /* Compute the Hessian */
265:     if (needH) {
266:       TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
267:     }

269:     /* Shift the Hessian matrix */
270:     if (pert > 0) {
271:       MatShift(tao->hessian, pert);
272:       if (tao->hessian != tao->hessian_pre) {
273:         MatShift(tao->hessian_pre, pert);
274:       }
275:     }

277:     if (nlsP->bfgs_pre) {
278:       MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
279:       ++bfgsUpdates;
280:     }

282:     /* Solve the Newton system of equations */
283:     KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre);
284:     if (is_nash || is_stcg || is_gltr) {
285:       KSPCGSetRadius(tao->ksp,nlsP->max_radius);
286:       KSPSolve(tao->ksp, tao->gradient, nlsP->D);
287:       KSPGetIterationNumber(tao->ksp,&kspits);
288:       tao->ksp_its+=kspits;
289:       tao->ksp_tot_its+=kspits;
290:       KSPCGGetNormD(tao->ksp,&norm_d);

292:       if (0.0 == tao->trust) {
293:         /* Radius was uninitialized; use the norm of the direction */
294:         if (norm_d > 0.0) {
295:           tao->trust = norm_d;

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:         } else {
301:           /* The direction was bad; set radius to default value and re-solve
302:              the trust-region subproblem to get a direction */
303:           tao->trust = tao->trust0;

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

309:           KSPCGSetRadius(tao->ksp,nlsP->max_radius);
310:           KSPSolve(tao->ksp, tao->gradient, nlsP->D);
311:           KSPGetIterationNumber(tao->ksp,&kspits);
312:           tao->ksp_its+=kspits;
313:           tao->ksp_tot_its+=kspits;
314:           KSPCGGetNormD(tao->ksp,&norm_d);

316:           if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Initial direction zero");
317:         }
318:       }
319:     } else {
320:       KSPSolve(tao->ksp, tao->gradient, nlsP->D);
321:       KSPGetIterationNumber(tao->ksp, &kspits);
322:       tao->ksp_its += kspits;
323:       tao->ksp_tot_its+=kspits;
324:     }
325:     VecScale(nlsP->D, -1.0);
326:     KSPGetConvergedReason(tao->ksp, &ksp_reason);
327:     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&  (nlsP->bfgs_pre)) {
328:       /* Preconditioner is numerically indefinite; reset the
329:          approximate if using BFGS preconditioning. */
330:       MatLMVMReset(nlsP->M, PETSC_FALSE);
331:       MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
332:       bfgsUpdates = 1;
333:     }

335:     if (KSP_CONVERGED_ATOL == ksp_reason) {
336:       ++nlsP->ksp_atol;
337:     } else if (KSP_CONVERGED_RTOL == ksp_reason) {
338:       ++nlsP->ksp_rtol;
339:     } else if (KSP_CONVERGED_CG_CONSTRAINED == ksp_reason) {
340:       ++nlsP->ksp_ctol;
341:     } else if (KSP_CONVERGED_CG_NEG_CURVE == ksp_reason) {
342:       ++nlsP->ksp_negc;
343:     } else if (KSP_DIVERGED_DTOL == ksp_reason) {
344:       ++nlsP->ksp_dtol;
345:     } else if (KSP_DIVERGED_ITS == ksp_reason) {
346:       ++nlsP->ksp_iter;
347:     } else {
348:       ++nlsP->ksp_othr;
349:     }

351:     /* Check for success (descent direction) */
352:     VecDot(nlsP->D, tao->gradient, &gdx);
353:     if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
354:       /* Newton step is not descent or direction produced Inf or NaN
355:          Update the perturbation for next time */
356:       if (pert <= 0.0) {
357:         /* Initialize the perturbation */
358:         pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
359:         if (is_gltr) {
360:           KSPGLTRGetMinEig(tao->ksp,&e_min);
361:           pert = PetscMax(pert, -e_min);
362:         }
363:       } else {
364:         /* Increase the perturbation */
365:         pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
366:       }

368:       if (!nlsP->bfgs_pre) {
369:         /* We don't have the bfgs matrix around and updated
370:            Must use gradient direction in this case */
371:         VecCopy(tao->gradient, nlsP->D);
372:         VecScale(nlsP->D, -1.0);
373:         ++nlsP->grad;
374:         stepType = NLS_GRADIENT;
375:       } else {
376:         /* Attempt to use the BFGS direction */
377:         MatSolve(nlsP->M, tao->gradient, nlsP->D);
378:         VecScale(nlsP->D, -1.0);

380:         /* Check for success (descent direction) */
381:         VecDot(tao->gradient, nlsP->D, &gdx);
382:         if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
383:           /* BFGS direction is not descent or direction produced not a number
384:              We can assert bfgsUpdates > 1 in this case because
385:              the first solve produces the scaled gradient direction,
386:              which is guaranteed to be descent */

388:           /* Use steepest descent direction (scaled) */
389:           MatLMVMReset(nlsP->M, PETSC_FALSE);
390:           MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
391:           MatSolve(nlsP->M, tao->gradient, nlsP->D);
392:           VecScale(nlsP->D, -1.0);

394:           bfgsUpdates = 1;
395:           ++nlsP->grad;
396:           stepType = NLS_GRADIENT;
397:         } else {
398:           MatLMVMGetUpdateCount(nlsP->M, &bfgsUpdates);
399:           if (1 == bfgsUpdates) {
400:             /* The first BFGS direction is always the scaled gradient */
401:             ++nlsP->grad;
402:             stepType = NLS_GRADIENT;
403:           } else {
404:             ++nlsP->bfgs;
405:             stepType = NLS_BFGS;
406:           }
407:         }
408:       }
409:     } else {
410:       /* Computed Newton step is descent */
411:       switch (ksp_reason) {
412:       case KSP_DIVERGED_NANORINF:
413:       case KSP_DIVERGED_BREAKDOWN:
414:       case KSP_DIVERGED_INDEFINITE_MAT:
415:       case KSP_DIVERGED_INDEFINITE_PC:
416:       case KSP_CONVERGED_CG_NEG_CURVE:
417:         /* Matrix or preconditioner is indefinite; increase perturbation */
418:         if (pert <= 0.0) {
419:           /* Initialize the perturbation */
420:           pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
421:           if (is_gltr) {
422:             KSPGLTRGetMinEig(tao->ksp, &e_min);
423:             pert = PetscMax(pert, -e_min);
424:           }
425:         } else {
426:           /* Increase the perturbation */
427:           pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
428:         }
429:         break;

431:       default:
432:         /* Newton step computation is good; decrease perturbation */
433:         pert = PetscMin(nlsP->psfac * pert, nlsP->pmsfac * gnorm);
434:         if (pert < nlsP->pmin) {
435:           pert = 0.0;
436:         }
437:         break;
438:       }

440:       ++nlsP->newt;
441:       stepType = NLS_NEWTON;
442:     }

444:     /* Perform the linesearch */
445:     fold = f;
446:     VecCopy(tao->solution, nlsP->Xold);
447:     VecCopy(tao->gradient, nlsP->Gold);

449:     TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason);
450:     TaoAddLineSearchCounts(tao);

452:     while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NLS_GRADIENT) {      /* Linesearch failed */
453:       f = fold;
454:       VecCopy(nlsP->Xold, tao->solution);
455:       VecCopy(nlsP->Gold, tao->gradient);

457:       switch(stepType) {
458:       case NLS_NEWTON:
459:         /* Failed to obtain acceptable iterate with Newton 1step
460:            Update the perturbation for next time */
461:         if (pert <= 0.0) {
462:           /* Initialize the perturbation */
463:           pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
464:           if (is_gltr) {
465:             KSPGLTRGetMinEig(tao->ksp,&e_min);
466:             pert = PetscMax(pert, -e_min);
467:           }
468:         } else {
469:           /* Increase the perturbation */
470:           pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
471:         }

473:         if (!nlsP->bfgs_pre) {
474:           /* We don't have the bfgs matrix around and being updated
475:              Must use gradient direction in this case */
476:           VecCopy(tao->gradient, nlsP->D);
477:           ++nlsP->grad;
478:           stepType = NLS_GRADIENT;
479:         } else {
480:           /* Attempt to use the BFGS direction */
481:           MatSolve(nlsP->M, tao->gradient, nlsP->D);
482:           /* Check for success (descent direction) */
483:           VecDot(tao->solution, nlsP->D, &gdx);
484:           if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
485:             /* BFGS direction is not descent or direction produced not a number
486:                We can assert bfgsUpdates > 1 in this case
487:                Use steepest descent direction (scaled) */
488:             MatLMVMReset(nlsP->M, PETSC_FALSE);
489:             MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
490:             MatSolve(nlsP->M, tao->gradient, nlsP->D);

492:             bfgsUpdates = 1;
493:             ++nlsP->grad;
494:             stepType = NLS_GRADIENT;
495:           } else {
496:             if (1 == bfgsUpdates) {
497:               /* The first BFGS direction is always the scaled gradient */
498:               ++nlsP->grad;
499:               stepType = NLS_GRADIENT;
500:             } else {
501:               ++nlsP->bfgs;
502:               stepType = NLS_BFGS;
503:             }
504:           }
505:         }
506:         break;

508:       case NLS_BFGS:
509:         /* Can only enter if pc_type == NLS_PC_BFGS
510:            Failed to obtain acceptable iterate with BFGS step
511:            Attempt to use the scaled gradient direction */
512:         MatLMVMReset(nlsP->M, PETSC_FALSE);
513:         MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient);
514:         MatSolve(nlsP->M, tao->gradient, nlsP->D);

516:         bfgsUpdates = 1;
517:         ++nlsP->grad;
518:         stepType = NLS_GRADIENT;
519:         break;
520:       }
521:       VecScale(nlsP->D, -1.0);

523:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason);
524:       TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full);
525:       TaoAddLineSearchCounts(tao);
526:     }

528:     if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
529:       /* Failed to find an improving point */
530:       f = fold;
531:       VecCopy(nlsP->Xold, tao->solution);
532:       VecCopy(nlsP->Gold, tao->gradient);
533:       step = 0.0;
534:       tao->reason = TAO_DIVERGED_LS_FAILURE;
535:       break;
536:     }

538:     /* Update trust region radius */
539:     if (is_nash || is_stcg || is_gltr) {
540:       switch(nlsP->update_type) {
541:       case NLS_UPDATE_STEP:
542:         if (stepType == NLS_NEWTON) {
543:           if (step < nlsP->nu1) {
544:             /* Very bad step taken; reduce radius */
545:             tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
546:           } else if (step < nlsP->nu2) {
547:             /* Reasonably bad step taken; reduce radius */
548:             tao->trust = nlsP->omega2 * PetscMin(norm_d, tao->trust);
549:           } else if (step < nlsP->nu3) {
550:             /*  Reasonable step was taken; leave radius alone */
551:             if (nlsP->omega3 < 1.0) {
552:               tao->trust = nlsP->omega3 * PetscMin(norm_d, tao->trust);
553:             } else if (nlsP->omega3 > 1.0) {
554:               tao->trust = PetscMax(nlsP->omega3 * norm_d, tao->trust);
555:             }
556:           } else if (step < nlsP->nu4) {
557:             /*  Full step taken; increase the radius */
558:             tao->trust = PetscMax(nlsP->omega4 * norm_d, tao->trust);
559:           } else {
560:             /*  More than full step taken; increase the radius */
561:             tao->trust = PetscMax(nlsP->omega5 * norm_d, tao->trust);
562:           }
563:         } else {
564:           /*  Newton step was not good; reduce the radius */
565:           tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
566:         }
567:         break;

569:       case NLS_UPDATE_REDUCTION:
570:         if (stepType == NLS_NEWTON) {
571:           /*  Get predicted reduction */

573:           KSPCGGetObjFcn(tao->ksp,&prered);
574:           if (prered >= 0.0) {
575:             /*  The predicted reduction has the wrong sign.  This cannot */
576:             /*  happen in infinite precision arithmetic.  Step should */
577:             /*  be rejected! */
578:             tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
579:           } else {
580:             if (PetscIsInfOrNanReal(f_full)) {
581:               tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
582:             } else {
583:               /*  Compute and actual reduction */
584:               actred = fold - f_full;
585:               prered = -prered;
586:               if ((PetscAbsScalar(actred) <= nlsP->epsilon) &&
587:                   (PetscAbsScalar(prered) <= nlsP->epsilon)) {
588:                 kappa = 1.0;
589:               } else {
590:                 kappa = actred / prered;
591:               }

593:               /*  Accept of reject the step and update radius */
594:               if (kappa < nlsP->eta1) {
595:                 /*  Very bad step */
596:                 tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
597:               } else if (kappa < nlsP->eta2) {
598:                 /*  Marginal bad step */
599:                 tao->trust = nlsP->alpha2 * PetscMin(tao->trust, norm_d);
600:               } else if (kappa < nlsP->eta3) {
601:                 /*  Reasonable step */
602:                 if (nlsP->alpha3 < 1.0) {
603:                   tao->trust = nlsP->alpha3 * PetscMin(norm_d, tao->trust);
604:                 } else if (nlsP->alpha3 > 1.0) {
605:                   tao->trust = PetscMax(nlsP->alpha3 * norm_d, tao->trust);
606:                 }
607:               } else if (kappa < nlsP->eta4) {
608:                 /*  Good step */
609:                 tao->trust = PetscMax(nlsP->alpha4 * norm_d, tao->trust);
610:               } else {
611:                 /*  Very good step */
612:                 tao->trust = PetscMax(nlsP->alpha5 * norm_d, tao->trust);
613:               }
614:             }
615:           }
616:         } else {
617:           /*  Newton step was not good; reduce the radius */
618:           tao->trust = nlsP->alpha1 * PetscMin(norm_d, tao->trust);
619:         }
620:         break;

622:       default:
623:         if (stepType == NLS_NEWTON) {
624:           KSPCGGetObjFcn(tao->ksp,&prered);
625:           if (prered >= 0.0) {
626:             /*  The predicted reduction has the wrong sign.  This cannot */
627:             /*  happen in infinite precision arithmetic.  Step should */
628:             /*  be rejected! */
629:             tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
630:           } else {
631:             if (PetscIsInfOrNanReal(f_full)) {
632:               tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
633:             } else {
634:               actred = fold - f_full;
635:               prered = -prered;
636:               if ((PetscAbsScalar(actred) <= nlsP->epsilon) && (PetscAbsScalar(prered) <= nlsP->epsilon)) {
637:                 kappa = 1.0;
638:               } else {
639:                 kappa = actred / prered;
640:               }

642:               tau_1 = nlsP->theta * gdx / (nlsP->theta * gdx - (1.0 - nlsP->theta) * prered + actred);
643:               tau_2 = nlsP->theta * gdx / (nlsP->theta * gdx + (1.0 + nlsP->theta) * prered - actred);
644:               tau_min = PetscMin(tau_1, tau_2);
645:               tau_max = PetscMax(tau_1, tau_2);

647:               if (kappa >= 1.0 - nlsP->mu1) {
648:                 /*  Great agreement */
649:                 if (tau_max < 1.0) {
650:                   tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
651:                 } else if (tau_max > nlsP->gamma4) {
652:                   tao->trust = PetscMax(tao->trust, nlsP->gamma4 * norm_d);
653:                 } else {
654:                   tao->trust = PetscMax(tao->trust, tau_max * norm_d);
655:                 }
656:               } else if (kappa >= 1.0 - nlsP->mu2) {
657:                 /*  Good agreement */

659:                 if (tau_max < nlsP->gamma2) {
660:                   tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
661:                 } else if (tau_max > nlsP->gamma3) {
662:                   tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
663:                 } else if (tau_max < 1.0) {
664:                   tao->trust = tau_max * PetscMin(tao->trust, norm_d);
665:                 } else {
666:                   tao->trust = PetscMax(tao->trust, tau_max * norm_d);
667:                 }
668:               } else {
669:                 /*  Not good agreement */
670:                 if (tau_min > 1.0) {
671:                   tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
672:                 } else if (tau_max < nlsP->gamma1) {
673:                   tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
674:                 } else if ((tau_min < nlsP->gamma1) && (tau_max >= 1.0)) {
675:                   tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
676:                 } else if ((tau_1 >= nlsP->gamma1) && (tau_1 < 1.0) && ((tau_2 < nlsP->gamma1) || (tau_2 >= 1.0))) {
677:                   tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
678:                 } else if ((tau_2 >= nlsP->gamma1) && (tau_2 < 1.0) && ((tau_1 < nlsP->gamma1) || (tau_2 >= 1.0))) {
679:                   tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
680:                 } else {
681:                   tao->trust = tau_max * PetscMin(tao->trust, norm_d);
682:                 }
683:               }
684:             }
685:           }
686:         } else {
687:           /*  Newton step was not good; reduce the radius */
688:           tao->trust = nlsP->gamma1 * PetscMin(norm_d, tao->trust);
689:         }
690:         break;
691:       }

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

697:     /*  Check for termination */
698:     TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
699:     if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER,"User provided compute function generated Not-a-Number");
700:     needH = 1;
701:     TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
702:     TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
703:     (*tao->ops->convergencetest)(tao,tao->cnvP);
704:   }
705:   return(0);
706: }

708: /* ---------------------------------------------------------- */
709: static PetscErrorCode TaoSetUp_NLS(Tao tao)
710: {
711:   TAO_NLS        *nlsP = (TAO_NLS *)tao->data;

715:   if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);}
716:   if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection);}
717:   if (!nlsP->W) {VecDuplicate(tao->solution,&nlsP->W);}
718:   if (!nlsP->D) {VecDuplicate(tao->solution,&nlsP->D);}
719:   if (!nlsP->Xold) {VecDuplicate(tao->solution,&nlsP->Xold);}
720:   if (!nlsP->Gold) {VecDuplicate(tao->solution,&nlsP->Gold);}
721:   nlsP->M = NULL;
722:   nlsP->bfgs_pre = NULL;
723:   return(0);
724: }

726: /*------------------------------------------------------------*/
727: static PetscErrorCode TaoDestroy_NLS(Tao tao)
728: {
729:   TAO_NLS        *nlsP = (TAO_NLS *)tao->data;

733:   if (tao->setupcalled) {
734:     VecDestroy(&nlsP->D);
735:     VecDestroy(&nlsP->W);
736:     VecDestroy(&nlsP->Xold);
737:     VecDestroy(&nlsP->Gold);
738:   }
739:   PetscFree(tao->data);
740:   return(0);
741: }

743: /*------------------------------------------------------------*/
744: static PetscErrorCode TaoSetFromOptions_NLS(PetscOptionItems *PetscOptionsObject,Tao tao)
745: {
746:   TAO_NLS        *nlsP = (TAO_NLS *)tao->data;

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


805: /*------------------------------------------------------------*/
806: static PetscErrorCode TaoView_NLS(Tao tao, PetscViewer viewer)
807: {
808:   TAO_NLS        *nlsP = (TAO_NLS *)tao->data;
809:   PetscBool      isascii;

813:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
814:   if (isascii) {
815:     PetscViewerASCIIPushTab(viewer);
816:     PetscViewerASCIIPrintf(viewer, "Newton steps: %D\n", nlsP->newt);
817:     PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", nlsP->bfgs);
818:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", nlsP->grad);

820:     PetscViewerASCIIPrintf(viewer, "nls ksp atol: %D\n", nlsP->ksp_atol);
821:     PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %D\n", nlsP->ksp_rtol);
822:     PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %D\n", nlsP->ksp_ctol);
823:     PetscViewerASCIIPrintf(viewer, "nls ksp negc: %D\n", nlsP->ksp_negc);
824:     PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %D\n", nlsP->ksp_dtol);
825:     PetscViewerASCIIPrintf(viewer, "nls ksp iter: %D\n", nlsP->ksp_iter);
826:     PetscViewerASCIIPrintf(viewer, "nls ksp othr: %D\n", nlsP->ksp_othr);
827:     PetscViewerASCIIPopTab(viewer);
828:   }
829:   return(0);
830: }

832: /* ---------------------------------------------------------- */
833: /*MC
834:   TAONLS - Newton's method with linesearch for unconstrained minimization.
835:   At each iteration, the Newton line search method solves the symmetric
836:   system of equations to obtain the step diretion dk:
837:               Hk dk = -gk
838:   a More-Thuente line search is applied on the direction dk to approximately
839:   solve
840:               min_t f(xk + t d_k)

842:     Options Database Keys:
843: + -tao_nls_init_type - "constant","direction","interpolation"
844: . -tao_nls_update_type - "step","direction","interpolation"
845: . -tao_nls_sval - perturbation starting value
846: . -tao_nls_imin - minimum initial perturbation
847: . -tao_nls_imax - maximum initial perturbation
848: . -tao_nls_pmin - minimum perturbation
849: . -tao_nls_pmax - maximum perturbation
850: . -tao_nls_pgfac - growth factor
851: . -tao_nls_psfac - shrink factor
852: . -tao_nls_imfac - initial merit factor
853: . -tao_nls_pmgfac - merit growth factor
854: . -tao_nls_pmsfac - merit shrink factor
855: . -tao_nls_eta1 - poor steplength; reduce radius
856: . -tao_nls_eta2 - reasonable steplength; leave radius
857: . -tao_nls_eta3 - good steplength; increase readius
858: . -tao_nls_eta4 - excellent steplength; greatly increase radius
859: . -tao_nls_alpha1 - alpha1 reduction
860: . -tao_nls_alpha2 - alpha2 reduction
861: . -tao_nls_alpha3 - alpha3 reduction
862: . -tao_nls_alpha4 - alpha4 reduction
863: . -tao_nls_alpha - alpha5 reduction
864: . -tao_nls_mu1 - mu1 interpolation update
865: . -tao_nls_mu2 - mu2 interpolation update
866: . -tao_nls_gamma1 - gamma1 interpolation update
867: . -tao_nls_gamma2 - gamma2 interpolation update
868: . -tao_nls_gamma3 - gamma3 interpolation update
869: . -tao_nls_gamma4 - gamma4 interpolation update
870: . -tao_nls_theta - theta interpolation update
871: . -tao_nls_omega1 - omega1 step update
872: . -tao_nls_omega2 - omega2 step update
873: . -tao_nls_omega3 - omega3 step update
874: . -tao_nls_omega4 - omega4 step update
875: . -tao_nls_omega5 - omega5 step update
876: . -tao_nls_mu1_i -  mu1 interpolation init factor
877: . -tao_nls_mu2_i -  mu2 interpolation init factor
878: . -tao_nls_gamma1_i -  gamma1 interpolation init factor
879: . -tao_nls_gamma2_i -  gamma2 interpolation init factor
880: . -tao_nls_gamma3_i -  gamma3 interpolation init factor
881: . -tao_nls_gamma4_i -  gamma4 interpolation init factor
882: - -tao_nls_theta_i -  theta interpolation init factor

884:   Level: beginner
885: M*/

887: PETSC_EXTERN PetscErrorCode TaoCreate_NLS(Tao tao)
888: {
889:   TAO_NLS        *nlsP;
890:   const char     *morethuente_type = TAOLINESEARCHMT;

894:   PetscNewLog(tao,&nlsP);

896:   tao->ops->setup = TaoSetUp_NLS;
897:   tao->ops->solve = TaoSolve_NLS;
898:   tao->ops->view = TaoView_NLS;
899:   tao->ops->setfromoptions = TaoSetFromOptions_NLS;
900:   tao->ops->destroy = TaoDestroy_NLS;

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

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

908:   nlsP->sval   = 0.0;
909:   nlsP->imin   = 1.0e-4;
910:   nlsP->imax   = 1.0e+2;
911:   nlsP->imfac  = 1.0e-1;

913:   nlsP->pmin   = 1.0e-12;
914:   nlsP->pmax   = 1.0e+2;
915:   nlsP->pgfac  = 1.0e+1;
916:   nlsP->psfac  = 4.0e-1;
917:   nlsP->pmgfac = 1.0e-1;
918:   nlsP->pmsfac = 1.0e-1;

920:   /*  Default values for trust-region radius update based on steplength */
921:   nlsP->nu1 = 0.25;
922:   nlsP->nu2 = 0.50;
923:   nlsP->nu3 = 1.00;
924:   nlsP->nu4 = 1.25;

926:   nlsP->omega1 = 0.25;
927:   nlsP->omega2 = 0.50;
928:   nlsP->omega3 = 1.00;
929:   nlsP->omega4 = 2.00;
930:   nlsP->omega5 = 4.00;

932:   /*  Default values for trust-region radius update based on reduction */
933:   nlsP->eta1 = 1.0e-4;
934:   nlsP->eta2 = 0.25;
935:   nlsP->eta3 = 0.50;
936:   nlsP->eta4 = 0.90;

938:   nlsP->alpha1 = 0.25;
939:   nlsP->alpha2 = 0.50;
940:   nlsP->alpha3 = 1.00;
941:   nlsP->alpha4 = 2.00;
942:   nlsP->alpha5 = 4.00;

944:   /*  Default values for trust-region radius update based on interpolation */
945:   nlsP->mu1 = 0.10;
946:   nlsP->mu2 = 0.50;

948:   nlsP->gamma1 = 0.25;
949:   nlsP->gamma2 = 0.50;
950:   nlsP->gamma3 = 2.00;
951:   nlsP->gamma4 = 4.00;

953:   nlsP->theta = 0.05;

955:   /*  Default values for trust region initialization based on interpolation */
956:   nlsP->mu1_i = 0.35;
957:   nlsP->mu2_i = 0.50;

959:   nlsP->gamma1_i = 0.0625;
960:   nlsP->gamma2_i = 0.5;
961:   nlsP->gamma3_i = 2.0;
962:   nlsP->gamma4_i = 5.0;

964:   nlsP->theta_i = 0.25;

966:   /*  Remaining parameters */
967:   nlsP->min_radius = 1.0e-10;
968:   nlsP->max_radius = 1.0e10;
969:   nlsP->epsilon = 1.0e-6;

971:   nlsP->init_type       = NLS_INIT_INTERPOLATION;
972:   nlsP->update_type     = NLS_UPDATE_STEP;

974:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
975:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch,(PetscObject)tao,1);
976:   TaoLineSearchSetType(tao->linesearch,morethuente_type);
977:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
978:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);

980:   /*  Set linear solver to default for symmetric matrices */
981:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
982:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1);
983:   KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
984:   KSPAppendOptionsPrefix(tao->ksp,"tao_nls_");
985:   KSPSetType(tao->ksp,KSPSTCG);
986:   return(0);
987: }