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

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

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

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

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

 89:   if (is_nash || is_stcg || is_gltr) {
 91:     tao->trust = tao->trust0;
 92:     tao->trust = PetscMax(tao->trust, nlsP->min_radius);
 93:     tao->trust = PetscMin(tao->trust, nlsP->max_radius);
 94:   }

 96:   /* Check convergence criteria */
 97:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 98:   TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);

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

107:   /* Allocate the vectors needed for the BFGS approximation */
108:   KSPGetPC(tao->ksp, &pc);
109:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);
110:   PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);
111:   if (is_bfgs) {
112:     nlsP->bfgs_pre = pc;
113:     PCLMVMGetMatLMVM(nlsP->bfgs_pre, &nlsP->M);
114:     VecGetLocalSize(tao->solution, &n);
115:     VecGetSize(tao->solution, &N);
116:     MatSetSizes(nlsP->M, n, n, N, N);
117:     MatLMVMAllocate(nlsP->M, tao->solution, tao->gradient);
118:     MatIsSymmetricKnown(nlsP->M, &sym_set, &is_symmetric);
120:   } else if (is_jacobi) {
121:     PCJacobiSetUseAbs(pc,PETSC_TRUE);
122:   }

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

132:     case NLS_INIT_INTERPOLATION:
133:       /* Use the initial radius specified */
134:       max_radius = 0.0;

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

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

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

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

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

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

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

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

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

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

224:           TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
226:           needH = 1;

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

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

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

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

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

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

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

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

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

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

295:           /* Modify the radius if it is too large or small */
296:           tao->trust = PetscMax(tao->trust, nlsP->min_radius);
297:           tao->trust = PetscMin(tao->trust, nlsP->max_radius);
298:         } else {
299:           /* The direction was bad; set radius to default value and re-solve
300:              the trust-region subproblem to get a direction */
301:           tao->trust = tao->trust0;

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

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

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

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

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

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

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

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

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

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

438:       ++nlsP->newt;
439:       stepType = NLS_NEWTON;
440:     }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

695:     /*  Check for termination */
696:     TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
698:     needH = 1;
699:     TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
700:     TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
701:     (*tao->ops->convergencetest)(tao,tao->cnvP);
702:   }
703:   return 0;
704: }

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

711:   if (!tao->gradient) VecDuplicate(tao->solution,&tao->gradient);
712:   if (!tao->stepdirection) VecDuplicate(tao->solution,&tao->stepdirection);
713:   if (!nlsP->W) VecDuplicate(tao->solution,&nlsP->W);
714:   if (!nlsP->D) VecDuplicate(tao->solution,&nlsP->D);
715:   if (!nlsP->Xold) VecDuplicate(tao->solution,&nlsP->Xold);
716:   if (!nlsP->Gold) VecDuplicate(tao->solution,&nlsP->Gold);
717:   nlsP->M = NULL;
718:   nlsP->bfgs_pre = NULL;
719:   return 0;
720: }

722: /*------------------------------------------------------------*/
723: static PetscErrorCode TaoDestroy_NLS(Tao tao)
724: {
725:   TAO_NLS        *nlsP = (TAO_NLS *)tao->data;

727:   if (tao->setupcalled) {
728:     VecDestroy(&nlsP->D);
729:     VecDestroy(&nlsP->W);
730:     VecDestroy(&nlsP->Xold);
731:     VecDestroy(&nlsP->Gold);
732:   }
733:   PetscFree(tao->data);
734:   return 0;
735: }

737: /*------------------------------------------------------------*/
738: static PetscErrorCode TaoSetFromOptions_NLS(PetscOptionItems *PetscOptionsObject,Tao tao)
739: {
740:   TAO_NLS        *nlsP = (TAO_NLS *)tao->data;

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

796: /*------------------------------------------------------------*/
797: static PetscErrorCode TaoView_NLS(Tao tao, PetscViewer viewer)
798: {
799:   TAO_NLS        *nlsP = (TAO_NLS *)tao->data;
800:   PetscBool      isascii;

802:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
803:   if (isascii) {
804:     PetscViewerASCIIPushTab(viewer);
805:     PetscViewerASCIIPrintf(viewer, "Newton steps: %D\n", nlsP->newt);
806:     PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", nlsP->bfgs);
807:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", nlsP->grad);

809:     PetscViewerASCIIPrintf(viewer, "nls ksp atol: %D\n", nlsP->ksp_atol);
810:     PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %D\n", nlsP->ksp_rtol);
811:     PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %D\n", nlsP->ksp_ctol);
812:     PetscViewerASCIIPrintf(viewer, "nls ksp negc: %D\n", nlsP->ksp_negc);
813:     PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %D\n", nlsP->ksp_dtol);
814:     PetscViewerASCIIPrintf(viewer, "nls ksp iter: %D\n", nlsP->ksp_iter);
815:     PetscViewerASCIIPrintf(viewer, "nls ksp othr: %D\n", nlsP->ksp_othr);
816:     PetscViewerASCIIPopTab(viewer);
817:   }
818:   return 0;
819: }

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

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

873:   Level: beginner
874: M*/

876: PETSC_EXTERN PetscErrorCode TaoCreate_NLS(Tao tao)
877: {
878:   TAO_NLS        *nlsP;
879:   const char     *morethuente_type = TAOLINESEARCHMT;

881:   PetscNewLog(tao,&nlsP);

883:   tao->ops->setup = TaoSetUp_NLS;
884:   tao->ops->solve = TaoSolve_NLS;
885:   tao->ops->view = TaoView_NLS;
886:   tao->ops->setfromoptions = TaoSetFromOptions_NLS;
887:   tao->ops->destroy = TaoDestroy_NLS;

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

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

895:   nlsP->sval   = 0.0;
896:   nlsP->imin   = 1.0e-4;
897:   nlsP->imax   = 1.0e+2;
898:   nlsP->imfac  = 1.0e-1;

900:   nlsP->pmin   = 1.0e-12;
901:   nlsP->pmax   = 1.0e+2;
902:   nlsP->pgfac  = 1.0e+1;
903:   nlsP->psfac  = 4.0e-1;
904:   nlsP->pmgfac = 1.0e-1;
905:   nlsP->pmsfac = 1.0e-1;

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

913:   nlsP->omega1 = 0.25;
914:   nlsP->omega2 = 0.50;
915:   nlsP->omega3 = 1.00;
916:   nlsP->omega4 = 2.00;
917:   nlsP->omega5 = 4.00;

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

925:   nlsP->alpha1 = 0.25;
926:   nlsP->alpha2 = 0.50;
927:   nlsP->alpha3 = 1.00;
928:   nlsP->alpha4 = 2.00;
929:   nlsP->alpha5 = 4.00;

931:   /*  Default values for trust-region radius update based on interpolation */
932:   nlsP->mu1 = 0.10;
933:   nlsP->mu2 = 0.50;

935:   nlsP->gamma1 = 0.25;
936:   nlsP->gamma2 = 0.50;
937:   nlsP->gamma3 = 2.00;
938:   nlsP->gamma4 = 4.00;

940:   nlsP->theta = 0.05;

942:   /*  Default values for trust region initialization based on interpolation */
943:   nlsP->mu1_i = 0.35;
944:   nlsP->mu2_i = 0.50;

946:   nlsP->gamma1_i = 0.0625;
947:   nlsP->gamma2_i = 0.5;
948:   nlsP->gamma3_i = 2.0;
949:   nlsP->gamma4_i = 5.0;

951:   nlsP->theta_i = 0.25;

953:   /*  Remaining parameters */
954:   nlsP->min_radius = 1.0e-10;
955:   nlsP->max_radius = 1.0e10;
956:   nlsP->epsilon = 1.0e-6;

958:   nlsP->init_type       = NLS_INIT_INTERPOLATION;
959:   nlsP->update_type     = NLS_UPDATE_STEP;

961:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
962:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch,(PetscObject)tao,1);
963:   TaoLineSearchSetType(tao->linesearch,morethuente_type);
964:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
965:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);

967:   /*  Set linear solver to default for symmetric matrices */
968:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
969:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1);
970:   KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
971:   KSPAppendOptionsPrefix(tao->ksp,"tao_nls_");
972:   KSPSetType(tao->ksp,KSPSTCG);
973:   return 0;
974: }