Actual source code: nls.c

  1: /*$Id$*/

  3: #include "nls.h"             /*I "tao_solver.h" I*/

  5: #ifdef TAO_USE_PETSC
  6: #include "petscksp.h"
  7: #include "petscpc.h"
  8: #include "src/petsctao/linearsolver/taolinearsolver_petsc.h"
  9: #include "src/petsctao/vector/taovec_petsc.h"

 11: #include "private/kspimpl.h"
 12: #include "private/pcimpl.h"

 14: #define NLS_KSP_CG        0
 15: #define NLS_KSP_NASH        1
 16: #define NLS_KSP_STCG        2
 17: #define NLS_KSP_GLTR        3
 18: #define NLS_KSP_PETSC        4
 19: #define NLS_KSP_TYPES        5

 21: #define NLS_PC_NONE        0
 22: #define NLS_PC_AHESS        1
 23: #define NLS_PC_BFGS        2
 24: #define NLS_PC_PETSC        3
 25: #define NLS_PC_TYPES        4

 27: #define BFGS_SCALE_AHESS        0
 28: #define BFGS_SCALE_PHESS        1
 29: #define BFGS_SCALE_BFGS                2
 30: #define BFGS_SCALE_TYPES        3

 32: #define NLS_INIT_CONSTANT         0
 33: #define NLS_INIT_DIRECTION        1
 34: #define NLS_INIT_INTERPOLATION    2
 35: #define NLS_INIT_TYPES            3

 37: #define NLS_UPDATE_STEP           0
 38: #define NLS_UPDATE_REDUCTION      1
 39: #define NLS_UPDATE_INTERPOLATION  2
 40: #define NLS_UPDATE_TYPES          3

 42: static const char *NLS_KSP[64] = {
 43:   "cg", "nash", "stcg", "gltr", "petsc"
 44: };

 46: static const char *NLS_PC[64] = {
 47:   "none", "ahess", "bfgs", "petsc"
 48: };

 50: static const char *BFGS_SCALE[64] = {
 51:   "ahess", "phess", "bfgs"
 52: };

 54: static const char *NLS_INIT[64] = {
 55:   "constant", "direction", "interpolation"
 56: };

 58: static const char *NLS_UPDATE[64] = {
 59:   "step", "reduction", "interpolation"
 60: };

 62: // Routine for BFGS preconditioner

 66: static PetscErrorCode bfgs_apply(PC pc, Vec xin, Vec xout)
 67: {
 68:   TaoLMVMMat *M ;
 69:   TaoVecPetsc Xin(xin);
 70:   TaoVecPetsc Xout(xout);
 71:   TaoTruth info2;
 72:   int info;

 75:   info = PCShellGetContext(pc,(void**)&M); CHKERRQ(info);
 76:   info = M->Solve(&Xin, &Xout, &info2); CHKERRQ(info);
 77:   return(0);
 78: }

 80: // Implements Newton's Method with a line search approach for solving
 81: // unconstrained minimization problems.  A More'-Thuente line search 
 82: // is used to guarantee that the bfgs preconditioner remains positive
 83: // definite.
 84: //
 85: // The method can shift the Hessian matrix.  The shifting procedure is
 86: // adapted from the PATH algorithm for solving complementarity
 87: // problems.
 88: //
 89: // The linear system solve should be done with a conjugate gradient
 90: // method, although any method can be used.

 92: #define NLS_NEWTON                 0
 93: #define NLS_BFGS                 1
 94: #define NLS_SCALED_GRADIENT         2
 95: #define NLS_GRADIENT                 3

 99: static int TaoSolve_NLS(TAO_SOLVER tao, void *solver)
100: {
101:   TAO_NLS *ls = (TAO_NLS *)solver;
102:   TaoVec *X, *G = ls->G, *D = ls->D, *W = ls->W;
103:   TaoVec *Xold = ls->Xold, *Gold = ls->Gold, *Diag = ls->Diag;
104:   TaoMat *H;
105:   TaoLMVMMat *M = ls->M;

107:   TaoLinearSolver *tls;
108:   TaoLinearSolverPetsc *pls;

110:   KSP pksp;
111:   PC ppc;

113:   KSPConvergedReason ksp_reason;
114:   TaoTerminateReason reason;
115:   TaoTruth success;
116:   
117:   double fmin, ftrial, f_full, prered, actred, kappa, sigma;
118:   double tau, tau_1, tau_2, tau_max, tau_min, max_radius;
119:   double f, fold, gdx, gnorm, pert;
120:   double step = 1.0;

122:   double delta;
123:   double radius, norm_d = 0.0, e_min;

125:   int info;
126:   TaoInt stepType;
127:   TaoInt iter = 0, status = 0;
128:   TaoInt bfgsUpdates = 0;
129:   TaoInt needH;

131:   TaoInt i_max = 5;
132:   TaoInt j_max = 1;
133:   TaoInt i, j;

135:   TaoFunctionBegin;
136:   // Initialized variables
137:   pert = ls->sval;

139:   ls->ksp_atol = 0;
140:   ls->ksp_rtol = 0;
141:   ls->ksp_dtol = 0;
142:   ls->ksp_ctol = 0;
143:   ls->ksp_negc = 0;
144:   ls->ksp_iter = 0;
145:   ls->ksp_othr = 0;

147:   // Initialize trust-region radius when using nash, stcg, or gltr
148:   // Will be reset during the first iteration
149:   if (NLS_KSP_NASH == ls->ksp_type ||
150:       NLS_KSP_STCG == ls->ksp_type || 
151:       NLS_KSP_GLTR == ls->ksp_type) {
152:     info = TaoGetInitialTrustRegionRadius(tao, &radius); CHKERRQ(info);
153:     if (radius < 0.0) {
154:       SETERRQ(1, "Initial radius negative");
155:     }

157:     // Modify the radius if it is too large or small
158:     radius = TaoMax(radius, ls->min_radius);
159:     radius = TaoMin(radius, ls->max_radius);
160:   }

162:   // Get vectors we will need
163:   info = TaoGetSolution(tao, &X); CHKERRQ(info);
164:   info = TaoGetHessian(tao, &H); CHKERRQ(info);

166:   if (NLS_PC_BFGS == ls->pc_type && !M) {
167:     ls->M = new TaoLMVMMat(X);
168:     M = ls->M;
169:   }

171:   // Check convergence criteria
172:   info = TaoComputeFunctionGradient(tao, X, &f, G); CHKERRQ(info);
173:   info = G->Norm2(&gnorm); CHKERRQ(info);
174:   if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
175:     SETERRQ(1, "User provided compute function generated Inf or NaN");
176:   }
177:   needH = 1;

179:   info = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); CHKERRQ(info);
180:   if (reason != TAO_CONTINUE_ITERATING) {
181:     TaoFunctionReturn(0);
182:   }

184:   // Create vectors for the limited memory preconditioner
185:   if ((NLS_PC_BFGS == ls->pc_type) && 
186:       (BFGS_SCALE_BFGS != ls->bfgs_scale_type)) {
187:     if (!Diag) {
188:       info = X->Clone(&ls->Diag); CHKERRQ(info);
189:       Diag = ls->Diag;
190:     }
191:   }

193:   // Modify the linear solver to a conjugate gradient method
194:   info = TaoGetLinearSolver(tao, &tls); CHKERRQ(info);
195:   pls  = dynamic_cast <TaoLinearSolverPetsc *> (tls);

197:   pksp = pls->GetKSP();
198:   switch(ls->ksp_type) {
199:   case NLS_KSP_CG:
200:     info = KSPSetType(pksp, KSPCG); CHKERRQ(info);
201:     if (pksp->ops->setfromoptions) {
202:       (*pksp->ops->setfromoptions)(pksp);
203:     }
204:     break;

206:   case NLS_KSP_NASH:
207:     info = KSPSetType(pksp, KSPNASH); CHKERRQ(info);
208:     if (pksp->ops->setfromoptions) {
209:       (*pksp->ops->setfromoptions)(pksp);
210:     }
211:     break;

213:   case NLS_KSP_STCG:
214:     info = KSPSetType(pksp, KSPSTCG); CHKERRQ(info);
215:     if (pksp->ops->setfromoptions) {
216:       (*pksp->ops->setfromoptions)(pksp);
217:     }
218:     break;

220:   case NLS_KSP_GLTR:
221:     info = KSPSetType(pksp, KSPGLTR); CHKERRQ(info);
222:     if (pksp->ops->setfromoptions) {
223:       (*pksp->ops->setfromoptions)(pksp);
224:     }
225:     break;

227:   default:
228:     // Use the method set by the ksp_type
229:     break;
230:   }

232:   // Modify the preconditioner to use the bfgs approximation
233:   info = KSPGetPC(pksp, &ppc); CHKERRQ(info);
234:   switch(ls->pc_type) {
235:   case NLS_PC_NONE:
236:     info = PCSetType(ppc, PCNONE); CHKERRQ(info);
237:     if (ppc->ops->setfromoptions) {
238:       (*ppc->ops->setfromoptions)(ppc);
239:     }
240:     break;

242:   case NLS_PC_AHESS:
243:     info = PCSetType(ppc, PCJACOBI); CHKERRQ(info);
244:     if (ppc->ops->setfromoptions) {
245:       (*ppc->ops->setfromoptions)(ppc);
246:     }
247:     info = PCJacobiSetUseAbs(ppc); CHKERRQ(info);
248:     break;

250:   case NLS_PC_BFGS:
251:     info = PCSetType(ppc, PCSHELL); CHKERRQ(info);
252:     if (ppc->ops->setfromoptions) {
253:       (*ppc->ops->setfromoptions)(ppc);
254:     }
255:     info = PCShellSetName(ppc, "bfgs"); CHKERRQ(info);
256:     info = PCShellSetContext(ppc, M); CHKERRQ(info);
257:     info = PCShellSetApply(ppc, bfgs_apply); CHKERRQ(info);
258:     break;

260:   default:
261:     // Use the pc method set by pc_type
262:     break;
263:   }

265:   // Initialize trust-region radius.  The initialization is only performed 
266:   // when we are using Steihaug-Toint or the Generalized Lanczos method.
267:   if (NLS_KSP_NASH == ls->ksp_type ||
268:       NLS_KSP_STCG == ls->ksp_type || 
269:       NLS_KSP_GLTR == ls->ksp_type) {
270:     switch(ls->init_type) {
271:     case NLS_INIT_CONSTANT:
272:       // Use the initial radius specified
273:       break;

275:     case NLS_INIT_INTERPOLATION:
276:       // Use the initial radius specified
277:       max_radius = 0.0;
278:   
279:       for (j = 0; j < j_max; ++j) {
280:         fmin = f;
281:         sigma = 0.0;
282:   
283:         if (needH) {
284:           info = TaoComputeHessian(tao, X, H); CHKERRQ(info);
285:           needH = 0;
286:         }
287:   
288:         for (i = 0; i < i_max; ++i) {
289:           info = W->Waxpby(1.0, X, -radius / gnorm, G); CHKERRQ(info);
290:   
291:           info = TaoComputeFunction(tao, W, &ftrial); CHKERRQ(info);
292:           if (TaoInfOrNaN(ftrial)) {
293:             tau = ls->gamma1_i;
294:           }
295:           else {
296:             if (ftrial < fmin) {
297:               fmin = ftrial;
298:               sigma = -radius / gnorm;
299:             }
300:   
301:             info = H->Multiply(G, D); CHKERRQ(info);
302:             info = D->Dot(G, &prered); CHKERRQ(info);
303:   
304:             prered = radius * (gnorm - 0.5 * radius * prered / (gnorm * gnorm));
305:             actred = f - ftrial;
306:             if ((fabs(actred) <= ls->epsilon) && 
307:                 (fabs(prered) <= ls->epsilon)) {
308:               kappa = 1.0;
309:             }
310:             else {
311:               kappa = actred / prered;
312:             }
313:   
314:             tau_1 = ls->theta_i * gnorm * radius / (ls->theta_i * gnorm * radius + (1.0 - ls->theta_i) * prered - actred);
315:             tau_2 = ls->theta_i * gnorm * radius / (ls->theta_i * gnorm * radius - (1.0 + ls->theta_i) * prered + actred);
316:             tau_min = TaoMin(tau_1, tau_2);
317:             tau_max = TaoMax(tau_1, tau_2);
318:   
319:             if (fabs(kappa - 1.0) <= ls->mu1_i) {
320:               // Great agreement
321:               max_radius = TaoMax(max_radius, radius);
322:   
323:               if (tau_max < 1.0) {
324:                 tau = ls->gamma3_i;
325:               }
326:               else if (tau_max > ls->gamma4_i) {
327:                 tau = ls->gamma4_i;
328:               }
329:               else if (tau_1 >= 1.0 && tau_1 <= ls->gamma4_i && tau_2 < 1.0) {
330:                 tau = tau_1;
331:               }
332:               else if (tau_2 >= 1.0 && tau_2 <= ls->gamma4_i && tau_1 < 1.0) {
333:                 tau = tau_2;
334:               }
335:               else {
336:                 tau = tau_max;
337:               }
338:             }
339:             else if (fabs(kappa - 1.0) <= ls->mu2_i) {
340:               // Good agreement
341:               max_radius = TaoMax(max_radius, radius);
342:   
343:               if (tau_max < ls->gamma2_i) {
344:                 tau = ls->gamma2_i;
345:               }
346:               else if (tau_max > ls->gamma3_i) {
347:                 tau = ls->gamma3_i;
348:               }
349:               else {
350:                 tau = tau_max;
351:               }
352:             }
353:             else {
354:               // Not good agreement
355:               if (tau_min > 1.0) {
356:                 tau = ls->gamma2_i;
357:               }
358:               else if (tau_max < ls->gamma1_i) {
359:                 tau = ls->gamma1_i;
360:               }
361:               else if ((tau_min < ls->gamma1_i) && (tau_max >= 1.0)) {
362:                 tau = ls->gamma1_i;
363:               }
364:               else if ((tau_1 >= ls->gamma1_i) && (tau_1 < 1.0) &&
365:                        ((tau_2 < ls->gamma1_i) || (tau_2 >= 1.0))) {
366:                 tau = tau_1;
367:               }
368:               else if ((tau_2 >= ls->gamma1_i) && (tau_2 < 1.0) &&
369:                        ((tau_1 < ls->gamma1_i) || (tau_2 >= 1.0))) {
370:                 tau = tau_2;
371:               }
372:               else {
373:                 tau = tau_max;
374:               }
375:             }
376:           }
377:           radius = tau * radius;
378:         }
379:   
380:         if (fmin < f) {
381:           f = fmin;
382:           info = X->Axpy(sigma, G); CHKERRQ(info);
383:           info = TaoComputeGradient(tao, X, G); CHKERRQ(info);
384:   
385:           info = G->Norm2(&gnorm); CHKERRQ(info);
386:           if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
387:             SETERRQ(1, "User provided compute function generated Inf or NaN");
388:           }
389:           needH = 1;
390:   
391:           info = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); CHKERRQ(info);
392:           if (reason != TAO_CONTINUE_ITERATING) {
393:             TaoFunctionReturn(0);
394:           }
395:         }
396:       }
397:       radius = TaoMax(radius, max_radius);

399:       // Modify the radius if it is too large or small
400:       radius = TaoMax(radius, ls->min_radius);
401:       radius = TaoMin(radius, ls->max_radius);
402:       break;

404:     default:
405:       // Norm of the first direction will initialize radius
406:       radius = 0.0;
407:       break;
408:     }
409:   } 

411:   // Set initial scaling for the BFGS preconditioner
412:   // This step is done after computing the initial trust-region radius
413:   // since the function value may have decreased
414:   if (NLS_PC_BFGS == ls->pc_type) {
415:     if (f != 0.0) {
416:       delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
417:     }
418:     else {
419:       delta = 2.0 / (gnorm*gnorm);
420:     }
421:     info = M->SetDelta(delta); CHKERRQ(info);
422:   }

424:   // Set counter for gradient/reset steps
425:   ls->newt = 0;
426:   ls->bfgs = 0;
427:   ls->sgrad = 0;
428:   ls->grad = 0;

430:   // Have not converged; continue with Newton method
431:   while (reason == TAO_CONTINUE_ITERATING) {
432:     ++iter;

434:     // Compute the Hessian
435:     if (needH) {
436:       info = TaoComputeHessian(tao, X, H); CHKERRQ(info);
437:       needH = 0;
438:     }

440:     if ((NLS_PC_BFGS == ls->pc_type) && 
441:         (BFGS_SCALE_AHESS == ls->bfgs_scale_type)) {
442:       // Obtain diagonal for the bfgs preconditioner 
443:       info = H->GetDiagonal(Diag); CHKERRQ(info);
444:       info = Diag->AbsoluteValue(); CHKERRQ(info);
445:       info = Diag->Reciprocal(); CHKERRQ(info);
446:       info = M->SetScale(Diag); CHKERRQ(info);
447:     }
448:  
449:     // Shift the Hessian matrix
450:     if (pert > 0) {
451:       info = H->ShiftDiagonal(pert); CHKERRQ(info);
452:     }
453:     
454:     if (NLS_PC_BFGS == ls->pc_type) {
455:       if (BFGS_SCALE_PHESS == ls->bfgs_scale_type) {
456:         // Obtain diagonal for the bfgs preconditioner 
457:         info = H->GetDiagonal(Diag); CHKERRQ(info);
458:         info = Diag->AbsoluteValue(); CHKERRQ(info);
459:         info = Diag->Reciprocal(); CHKERRQ(info);
460:         info = M->SetScale(Diag); CHKERRQ(info);
461:       }

463:       // Update the limited memory preconditioner
464:       info = M->Update(X, G); CHKERRQ(info);
465:       ++bfgsUpdates;
466:     }

468:     // Solve the Newton system of equations
469:     info = TaoPreLinearSolve(tao, H); CHKERRQ(info);
470:     if (NLS_KSP_NASH == ls->ksp_type ||
471:         NLS_KSP_STCG == ls->ksp_type || 
472:         NLS_KSP_GLTR == ls->ksp_type) {
473:       info = TaoLinearSolveTrustRegion(tao, H, G, D, radius, &success); CHKERRQ(info);
474:       info = pls->GetNormDirection(&norm_d); CHKERRQ(info);
475:       if (0.0 == radius) {
476:         // Radius was uninitialized; use the norm of the direction
477:         if (norm_d > 0.0) {
478:           radius = norm_d;

480:           // Modify the radius if it is too large or small
481:           radius = TaoMax(radius, ls->min_radius);
482:           radius = TaoMin(radius, ls->max_radius);
483:         }
484:         else {
485:           // The direction was bad; set radius to default value and re-solve 
486:           // the trust-region subproblem to get a direction
487:           info = TaoGetInitialTrustRegionRadius(tao, &radius); CHKERRQ(info);

489:           // Modify the radius if it is too large or small
490:           radius = TaoMax(radius, ls->min_radius);
491:           radius = TaoMin(radius, ls->max_radius);

493:           info = TaoLinearSolveTrustRegion(tao, H, G, D, radius, &success); CHKERRQ(info);
494:           info = pls->GetNormDirection(&norm_d); CHKERRQ(info);
495:           if (norm_d == 0.0) {
496:             SETERRQ(1, "Initial direction zero");
497:           }
498:         }
499:       }
500:     }
501:     else {
502:       info = TaoLinearSolve(tao, H, G, D, &success); CHKERRQ(info);
503:     }
504:     info = D->Negate(); CHKERRQ(info);

506:     info = KSPGetConvergedReason(pksp, &ksp_reason); CHKERRQ(info);
507:     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
508:         (NLS_PC_BFGS == ls->pc_type) && (bfgsUpdates > 1)) {
509:       // Preconditioner is numerically indefinite; reset the
510:       // approximate if using BFGS preconditioning.

512:       if (f != 0.0) {
513:         delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
514:       }
515:       else {
516:         delta = 2.0 / (gnorm*gnorm);
517:       }
518:       info = M->SetDelta(delta); CHKERRQ(info);
519:       info = M->Reset(); CHKERRQ(info);
520:       info = M->Update(X, G); CHKERRQ(info);
521:       bfgsUpdates = 1;
522:     }

524:     if (KSP_CONVERGED_ATOL == ksp_reason) {
525:       ++ls->ksp_atol;
526:     }
527:     else if (KSP_CONVERGED_RTOL == ksp_reason) {
528:       ++ls->ksp_rtol;
529:     }
530:     else if (KSP_CONVERGED_CG_CONSTRAINED == ksp_reason) {
531:       ++ls->ksp_ctol;
532:     }
533:     else if (KSP_CONVERGED_CG_NEG_CURVE == ksp_reason) {
534:       ++ls->ksp_negc;
535:     }
536:     else if (KSP_DIVERGED_DTOL == ksp_reason) {
537:       ++ls->ksp_dtol;
538:     }
539:     else if (KSP_DIVERGED_ITS == ksp_reason) {
540:       ++ls->ksp_iter;
541:     }
542:     else {
543:       ++ls->ksp_othr;
544:     } 

546:     // Check for success (descent direction)
547:     info = D->Dot(G, &gdx); CHKERRQ(info);
548:     if ((gdx >= 0.0) || TaoInfOrNaN(gdx)) {
549:       // Newton step is not descent or direction produced Inf or NaN
550:       // Update the perturbation for next time
551:       if (pert <= 0.0) {
552:         // Initialize the perturbation
553:         pert = TaoMin(ls->imax, TaoMax(ls->imin, ls->imfac * gnorm));
554:         if (NLS_KSP_GLTR == ls->ksp_type) {
555:           info = pls->GetMinEig(&e_min); CHKERRQ(info);
556:           pert = TaoMax(pert, -e_min);
557:         }
558:       }
559:       else {
560:         // Increase the perturbation
561:         pert = TaoMin(ls->pmax, TaoMax(ls->pgfac * pert, ls->pmgfac * gnorm));
562:       }

564:       if (NLS_PC_BFGS != ls->pc_type) {
565:         // We don't have the bfgs matrix around and updated
566:         // Must use gradient direction in this case
567:         info = D->CopyFrom(G); CHKERRQ(info);
568:         info = D->Negate(); CHKERRQ(info);

570:         ++ls->grad;
571:         stepType = NLS_GRADIENT;
572:       }
573:       else {
574:         // Attempt to use the BFGS direction
575:         info = M->Solve(G, D, &success); CHKERRQ(info);
576:         info = D->Negate(); CHKERRQ(info);

578:         // Check for success (descent direction)
579:         info = D->Dot(G, &gdx); CHKERRQ(info);
580:         if ((gdx >= 0) || TaoInfOrNaN(gdx)) {
581:           // BFGS direction is not descent or direction produced not a number
582:           // We can assert bfgsUpdates > 1 in this case because
583:           // the first solve produces the scaled gradient direction,
584:           // which is guaranteed to be descent
585:           //
586:           // Use steepest descent direction (scaled)

588:           if (f != 0.0) {
589:             delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
590:           }
591:           else {
592:             delta = 2.0 / (gnorm*gnorm);
593:           }
594:           info = M->SetDelta(delta); CHKERRQ(info);
595:           info = M->Reset(); CHKERRQ(info);
596:           info = M->Update(X, G); CHKERRQ(info);
597:           info = M->Solve(G, D, &success); CHKERRQ(info);
598:           info = D->Negate(); CHKERRQ(info);
599:   
600:           bfgsUpdates = 1;
601:           ++ls->sgrad;
602:           stepType = NLS_SCALED_GRADIENT;
603:         }
604:         else {
605:           if (1 == bfgsUpdates) {
606:             // The first BFGS direction is always the scaled gradient
607:             ++ls->sgrad;
608:             stepType = NLS_SCALED_GRADIENT;
609:           }
610:           else {
611:             ++ls->bfgs;
612:             stepType = NLS_BFGS;
613:           }
614:         }
615:       }
616:     }
617:     else {
618:       // Computed Newton step is descent
619:       switch (ksp_reason) {
620:       case KSP_DIVERGED_NAN:
621:       case KSP_DIVERGED_BREAKDOWN:
622:       case KSP_DIVERGED_INDEFINITE_MAT:
623:       case KSP_DIVERGED_INDEFINITE_PC:
624:       case KSP_CONVERGED_CG_NEG_CURVE:
625:         // Matrix or preconditioner is indefinite; increase perturbation
626:         if (pert <= 0.0) {
627:           // Initialize the perturbation
628:           pert = TaoMin(ls->imax, TaoMax(ls->imin, ls->imfac * gnorm));
629:           if (NLS_KSP_GLTR == ls->ksp_type) {
630:             info = pls->GetMinEig(&e_min); CHKERRQ(info);
631:             pert = TaoMax(pert, -e_min);
632:           }
633:         }
634:         else {
635:           // Increase the perturbation
636:           pert = TaoMin(ls->pmax, TaoMax(ls->pgfac * pert, ls->pmgfac * gnorm));
637:         }
638:         break;

640:       default:
641:         // Newton step computation is good; decrease perturbation
642:         pert = TaoMin(ls->psfac * pert, ls->pmsfac * gnorm);
643:         if (pert < ls->pmin) {
644:           pert = 0.0;
645:         }
646:         break; 
647:       }

649:       ++ls->newt;
650:       stepType = NLS_NEWTON;
651:     }

653:     // Perform the linesearch
654:     fold = f;
655:     info = Xold->CopyFrom(X); CHKERRQ(info);
656:     info = Gold->CopyFrom(G); CHKERRQ(info);

658:     step = 1.0;
659:     info = TaoLineSearchApply(tao, X, G, D, W, &f, &f_full, &step, &status); CHKERRQ(info);

661:     while (status && stepType != NLS_GRADIENT) {
662:       // Linesearch failed
663:       f = fold;
664:       info = X->CopyFrom(Xold); CHKERRQ(info);
665:       info = G->CopyFrom(Gold); CHKERRQ(info);

667:       switch(stepType) {
668:       case NLS_NEWTON:
669:         // Failed to obtain acceptable iterate with Newton step
670:         // Update the perturbation for next time
671:         if (pert <= 0.0) {
672:           // Initialize the perturbation
673:           pert = TaoMin(ls->imax, TaoMax(ls->imin, ls->imfac * gnorm));
674:           if (NLS_KSP_GLTR == ls->ksp_type) {
675:             info = pls->GetMinEig(&e_min); CHKERRQ(info);
676:             pert = TaoMax(pert, -e_min);
677:           }
678:         }
679:         else {
680:           // Increase the perturbation
681:           pert = TaoMin(ls->pmax, TaoMax(ls->pgfac * pert, ls->pmgfac * gnorm));
682:         }

684:         if (NLS_PC_BFGS != ls->pc_type) {
685:           // We don't have the bfgs matrix around and being updated
686:           // Must use gradient direction in this case
687:           info = D->CopyFrom(G); CHKERRQ(info);

689:           ++ls->grad;
690:           stepType = NLS_GRADIENT;
691:         }
692:         else {
693:           // Attempt to use the BFGS direction
694:           info = M->Solve(G, D, &success); CHKERRQ(info);

696:           // Check for success (descent direction)
697:           info = D->Dot(G, &gdx); CHKERRQ(info);
698:           if ((gdx <= 0) || TaoInfOrNaN(gdx)) {
699:             // BFGS direction is not descent or direction produced not a number
700:             // We can assert bfgsUpdates > 1 in this case
701:             // Use steepest descent direction (scaled)
702:     
703:             if (f != 0.0) {
704:               delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
705:             }
706:             else {
707:               delta = 2.0 / (gnorm*gnorm);
708:             }
709:             info = M->SetDelta(delta); CHKERRQ(info);
710:             info = M->Reset(); CHKERRQ(info);
711:             info = M->Update(X, G); CHKERRQ(info);
712:             info = M->Solve(G, D, &success); CHKERRQ(info);
713:   
714:             bfgsUpdates = 1;
715:             ++ls->sgrad;
716:             stepType = NLS_SCALED_GRADIENT;
717:           }
718:           else {
719:             if (1 == bfgsUpdates) {
720:               // The first BFGS direction is always the scaled gradient
721:               ++ls->sgrad;
722:               stepType = NLS_SCALED_GRADIENT;
723:             }
724:             else {
725:               ++ls->bfgs;
726:               stepType = NLS_BFGS;
727:             }
728:           }
729:         }
730:         break;

732:       case NLS_BFGS:
733:         // Can only enter if pc_type == NLS_PC_BFGS
734:         // Failed to obtain acceptable iterate with BFGS step
735:         // Attempt to use the scaled gradient direction

737:         if (f != 0.0) {
738:           delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
739:         }
740:         else {
741:           delta = 2.0 / (gnorm*gnorm);
742:         }
743:         info = M->SetDelta(delta); CHKERRQ(info);
744:         info = M->Reset(); CHKERRQ(info);
745:         info = M->Update(X, G); CHKERRQ(info);
746:         info = M->Solve(G, D, &success); CHKERRQ(info);

748:         bfgsUpdates = 1;
749:         ++ls->sgrad;
750:         stepType = NLS_SCALED_GRADIENT;
751:         break;

753:       case NLS_SCALED_GRADIENT:
754:         // Can only enter if pc_type == NLS_PC_BFGS
755:         // The scaled gradient step did not produce a new iterate;
756:         // attemp to use the gradient direction.
757:         // Need to make sure we are not using a different diagonal scaling
758:         info = M->SetScale(0); CHKERRQ(info);
759:         info = M->SetDelta(1.0); CHKERRQ(info);
760:         info = M->Reset(); CHKERRQ(info);
761:         info = M->Update(X, G); CHKERRQ(info);
762:         info = M->Solve(G, D, &success); CHKERRQ(info);

764:         bfgsUpdates = 1;
765:         ++ls->grad;
766:         stepType = NLS_GRADIENT;
767:         break;
768:       }
769:       info = D->Negate(); CHKERRQ(info);

771:       // This may be incorrect; linesearch has values for stepmax and stepmin
772:       // that should be reset.
773:       step = 1.0;
774:       info = TaoLineSearchApply(tao, X, G, D, W, &f, &f_full, &step, &status); CHKERRQ(info);
775:     }

777:     if (status) {
778:       // Failed to find an improving point
779:       f = fold;
780:       info = X->CopyFrom(Xold); CHKERRQ(info);
781:       info = G->CopyFrom(Gold); CHKERRQ(info);
782:       step = 0.0;
783:     }

785:     // Update trust region radius
786:     if (NLS_KSP_NASH == ls->ksp_type ||
787:         NLS_KSP_STCG == ls->ksp_type || 
788:         NLS_KSP_GLTR == ls->ksp_type) {
789:       switch(ls->update_type) {
790:       case NLS_UPDATE_STEP:
791:         if (stepType == NLS_NEWTON) {
792:           if (step < ls->nu1) {
793:             // Very bad step taken; reduce radius
794:             radius = ls->omega1 * TaoMin(norm_d, radius);
795:           }
796:           else if (step < ls->nu2) {
797:             // Reasonably bad step taken; reduce radius
798:             radius = ls->omega2 * TaoMin(norm_d, radius);
799:           }
800:           else if (step < ls->nu3) {
801:             // Reasonable step was taken; leave radius alone
802:             if (ls->omega3 < 1.0) {
803:               radius = ls->omega3 * TaoMin(norm_d, radius);
804:             }
805:             else if (ls->omega3 > 1.0) {
806:               radius = TaoMax(ls->omega3 * norm_d, radius);  
807:             }
808:           }
809:           else if (step < ls->nu4) {
810:             // Full step taken; increase the radius
811:             radius = TaoMax(ls->omega4 * norm_d, radius);  
812:           }
813:           else {
814:             // More than full step taken; increase the radius
815:             radius = TaoMax(ls->omega5 * norm_d, radius);  
816:           }
817:         }
818:         else {
819:           // Newton step was not good; reduce the radius
820:           radius = ls->omega1 * TaoMin(norm_d, radius);
821:         }
822:         break;

824:       case NLS_UPDATE_REDUCTION:
825:         if (stepType == NLS_NEWTON) {
826:           // Get predicted reduction
827:           info = pls->GetObjFcn(&prered); CHKERRQ(info);

829:           if (prered >= 0.0) {
830:             // The predicted reduction has the wrong sign.  This cannot
831:             // happen in infinite precision arithmetic.  Step should
832:             // be rejected!
833:             radius = ls->alpha1 * TaoMin(radius, norm_d);
834:           }
835:           else {
836:             if (TaoInfOrNaN(f_full)) {
837:               radius = ls->alpha1 * TaoMin(radius, norm_d);
838:             }
839:             else {
840:               // Compute and actual reduction
841:               actred = fold - f_full;
842:               prered = -prered;
843:               if ((fabs(actred) <= ls->epsilon) && 
844:                   (fabs(prered) <= ls->epsilon)) {
845:                 kappa = 1.0;
846:               }
847:               else {
848:                 kappa = actred / prered;
849:               }
850:   
851:               // Accept of reject the step and update radius
852:               if (kappa < ls->eta1) {
853:                 // Very bad step
854:                 radius = ls->alpha1 * TaoMin(radius, norm_d);
855:               }
856:               else if (kappa < ls->eta2) {
857:                 // Marginal bad step
858:                 radius = ls->alpha2 * TaoMin(radius, norm_d);
859:               }
860:               else if (kappa < ls->eta3) {
861:                 // Reasonable step
862:                 if (ls->alpha3 < 1.0) {
863:                   radius = ls->alpha3 * TaoMin(norm_d, radius);
864:                 }
865:                 else if (ls->alpha3 > 1.0) {
866:                   radius = TaoMax(ls->alpha3 * norm_d, radius);  
867:                 }
868:               }
869:               else if (kappa < ls->eta4) {
870:                 // Good step
871:                 radius = TaoMax(ls->alpha4 * norm_d, radius);
872:               }
873:               else {
874:                 // Very good step
875:                 radius = TaoMax(ls->alpha5 * norm_d, radius);
876:               }
877:             }
878:           }
879:         }
880:         else {
881:           // Newton step was not good; reduce the radius
882:           radius = ls->alpha1 * TaoMin(norm_d, radius);
883:         }
884:         break;

886:       default:
887:         if (stepType == NLS_NEWTON) {
888:           // Get predicted reduction
889:           info = pls->GetObjFcn(&prered); CHKERRQ(info);

891:           if (prered >= 0.0) {
892:             // The predicted reduction has the wrong sign.  This cannot
893:             // happen in infinite precision arithmetic.  Step should
894:             // be rejected!
895:             radius = ls->gamma1 * TaoMin(radius, norm_d);
896:           }
897:           else {
898:             if (TaoInfOrNaN(f_full)) {
899:               radius = ls->gamma1 * TaoMin(radius, norm_d);
900:             }
901:             else {
902:               actred = fold - f_full;
903:               prered = -prered;
904:               if ((fabs(actred) <= ls->epsilon) && 
905:                   (fabs(prered) <= ls->epsilon)) {
906:                 kappa = 1.0;
907:               }
908:               else {
909:                 kappa = actred / prered;
910:               }

912:               tau_1 = ls->theta * gdx / (ls->theta * gdx - (1.0 - ls->theta) * prered + actred);
913:               tau_2 = ls->theta * gdx / (ls->theta * gdx + (1.0 + ls->theta) * prered - actred);
914:               tau_min = TaoMin(tau_1, tau_2);
915:               tau_max = TaoMax(tau_1, tau_2);

917:               if (kappa >= 1.0 - ls->mu1) {
918:                 // Great agreement
919:                 if (tau_max < 1.0) {
920:                   radius = TaoMax(radius, ls->gamma3 * norm_d);
921:                 }
922:                 else if (tau_max > ls->gamma4) {
923:                   radius = TaoMax(radius, ls->gamma4 * norm_d);
924:                 }
925:                 else {
926:                   radius = TaoMax(radius, tau_max * norm_d);
927:                 }
928:               }
929:               else if (kappa >= 1.0 - ls->mu2) {
930:                 // Good agreement

932:                 if (tau_max < ls->gamma2) {
933:                   radius = ls->gamma2 * TaoMin(radius, norm_d);
934:                 }
935:                 else if (tau_max > ls->gamma3) {
936:                   radius = TaoMax(radius, ls->gamma3 * norm_d);
937:                 }
938:                 else if (tau_max < 1.0) {
939:                   radius = tau_max * TaoMin(radius, norm_d);
940:                 }
941:                 else {
942:                   radius = TaoMax(radius, tau_max * norm_d);
943:                 }
944:               }
945:               else {
946:                 // Not good agreement
947:                 if (tau_min > 1.0) {
948:                   radius = ls->gamma2 * TaoMin(radius, norm_d);
949:                 }
950:                 else if (tau_max < ls->gamma1) {
951:                   radius = ls->gamma1 * TaoMin(radius, norm_d);
952:                 }
953:                 else if ((tau_min < ls->gamma1) && (tau_max >= 1.0)) {
954:                   radius = ls->gamma1 * TaoMin(radius, norm_d);
955:                 }
956:                 else if ((tau_1 >= ls->gamma1) && (tau_1 < 1.0) &&
957:                          ((tau_2 < ls->gamma1) || (tau_2 >= 1.0))) {
958:                   radius = tau_1 * TaoMin(radius, norm_d);
959:                 }
960:                 else if ((tau_2 >= ls->gamma1) && (tau_2 < 1.0) &&
961:                          ((tau_1 < ls->gamma1) || (tau_2 >= 1.0))) {
962:                   radius = tau_2 * TaoMin(radius, norm_d);
963:                 }
964:                 else {
965:                   radius = tau_max * TaoMin(radius, norm_d);
966:                 }
967:               }
968:             } 
969:           }
970:         }
971:         else {
972:           // Newton step was not good; reduce the radius
973:           radius = ls->gamma1 * TaoMin(norm_d, radius);
974:         }
975:         break;
976:       }

978:       // The radius may have been increased; modify if it is too large
979:       radius = TaoMin(radius, ls->max_radius);
980:     }

982:     // Check for termination
983:     info = G->Norm2(&gnorm); CHKERRQ(info);
984:     if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
985:       SETERRQ(1,"User provided compute function generated Not-a-Number");
986:     }
987:     needH = 1;

989:     info = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason); CHKERRQ(info);
990:   }
991:   TaoFunctionReturn(0);
992: }

994: /* ---------------------------------------------------------- */
997: static int TaoSetUp_NLS(TAO_SOLVER tao, void *solver)
998: {
999:   TAO_NLS *ls = (TAO_NLS *)solver;
1000:   TaoVec *X;
1001:   TaoMat *H;
1002:   int info;

1004:   TaoFunctionBegin;

1006:   info = TaoGetSolution(tao, &X); CHKERRQ(info);
1007:   info = X->Clone(&ls->G); CHKERRQ(info);
1008:   info = X->Clone(&ls->D); CHKERRQ(info);
1009:   info = X->Clone(&ls->W); CHKERRQ(info);

1011:   info = X->Clone(&ls->Xold); CHKERRQ(info);
1012:   info = X->Clone(&ls->Gold); CHKERRQ(info);

1014:   ls->Diag = 0;
1015:   ls->M = 0;

1017:   info = TaoSetLagrangianGradientVector(tao, ls->G); CHKERRQ(info);
1018:   info = TaoSetStepDirectionVector(tao, ls->D); CHKERRQ(info);

1020:   // Set linear solver to default for symmetric matrices
1021:   info = TaoGetHessian(tao, &H); CHKERRQ(info);
1022:   info = TaoCreateLinearSolver(tao, H, 200, 0); CHKERRQ(info);

1024:   // Check sizes for compatability
1025:   info = TaoCheckFGH(tao); CHKERRQ(info);
1026:   TaoFunctionReturn(0);
1027: }

1029: /*------------------------------------------------------------*/
1032: static int TaoDestroy_NLS(TAO_SOLVER tao, void *solver)
1033: {
1034:   TAO_NLS *ls = (TAO_NLS *)solver;
1035:   int info;

1037:   TaoFunctionBegin;
1038:   info = TaoVecDestroy(ls->G); CHKERRQ(info);
1039:   info = TaoVecDestroy(ls->D); CHKERRQ(info);
1040:   info = TaoVecDestroy(ls->W); CHKERRQ(info);

1042:   info = TaoVecDestroy(ls->Xold); CHKERRQ(info);
1043:   info = TaoVecDestroy(ls->Gold); CHKERRQ(info);

1045:   info = TaoSetLagrangianGradientVector(tao, 0); CHKERRQ(info);
1046:   info = TaoSetStepDirectionVector(tao, 0); CHKERRQ(info);

1048:   if (ls->Diag) {
1049:     info = TaoVecDestroy(ls->Diag); CHKERRQ(info);
1050:     ls->Diag = 0;
1051:   }

1053:   if (ls->M) {
1054:     info = TaoMatDestroy(ls->M); CHKERRQ(info);
1055:     ls->M = 0;
1056:   }
1057:   TaoFunctionReturn(0);
1058: }

1060: /*------------------------------------------------------------*/
1063: static int TaoSetOptions_NLS(TAO_SOLVER tao, void *solver)
1064: {
1065:   TAO_NLS *ls = (TAO_NLS *)solver;
1066:   int info;

1068:   TaoFunctionBegin;
1069:   info = TaoOptionsHead("Newton line search method for unconstrained optimization"); CHKERRQ(info);
1070:   info = TaoOptionList("-tao_nls_ksp_type", "ksp type", "", NLS_KSP, NLS_KSP_TYPES, NLS_KSP[ls->ksp_type], &ls->ksp_type, 0); CHKERRQ(info);
1071:   info = TaoOptionList("-tao_nls_pc_type", "pc type", "", NLS_PC, NLS_PC_TYPES, NLS_PC[ls->pc_type], &ls->pc_type, 0); CHKERRQ(info);
1072:   info = TaoOptionList("-tao_nls_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[ls->bfgs_scale_type], &ls->bfgs_scale_type, 0); CHKERRQ(info);
1073:   info = TaoOptionList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[ls->init_type], &ls->init_type, 0); CHKERRQ(info);
1074:   info = TaoOptionList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[ls->update_type], &ls->update_type, 0); CHKERRQ(info);
1075:  info = TaoOptionDouble("-tao_nls_sval", "perturbation starting value", "", ls->sval, &ls->sval, 0); CHKERRQ(info);
1076:   info = TaoOptionDouble("-tao_nls_imin", "minimum initial perturbation", "", ls->imin, &ls->imin, 0); CHKERRQ(info);
1077:   info = TaoOptionDouble("-tao_nls_imax", "maximum initial perturbation", "", ls->imax, &ls->imax, 0); CHKERRQ(info);
1078:   info = TaoOptionDouble("-tao_nls_imfac", "initial merit factor", "", ls->imfac, &ls->imfac, 0); CHKERRQ(info);
1079:   info = TaoOptionDouble("-tao_nls_pmin", "minimum perturbation", "", ls->pmin, &ls->pmin, 0); CHKERRQ(info);
1080:   info = TaoOptionDouble("-tao_nls_pmax", "maximum perturbation", "", ls->pmax, &ls->pmax, 0); CHKERRQ(info);
1081:   info = TaoOptionDouble("-tao_nls_pgfac", "growth factor", "", ls->pgfac, &ls->pgfac, 0); CHKERRQ(info);
1082:   info = TaoOptionDouble("-tao_nls_psfac", "shrink factor", "", ls->psfac, &ls->psfac, 0); CHKERRQ(info);
1083:   info = TaoOptionDouble("-tao_nls_pmgfac", "merit growth factor", "", ls->pmgfac, &ls->pmgfac, 0); CHKERRQ(info);
1084:   info = TaoOptionDouble("-tao_nls_pmsfac", "merit shrink factor", "", ls->pmsfac, &ls->pmsfac, 0); CHKERRQ(info);
1085:   info = TaoOptionDouble("-tao_nls_eta1", "poor steplength; reduce radius", "", ls->eta1, &ls->eta1, 0); CHKERRQ(info);
1086:   info = TaoOptionDouble("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", ls->eta2, &ls->eta2, 0); CHKERRQ(info);
1087:   info = TaoOptionDouble("-tao_nls_eta3", "good steplength; increase radius", "", ls->eta3, &ls->eta3, 0); CHKERRQ(info);
1088:   info = TaoOptionDouble("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", ls->eta4, &ls->eta4, 0); CHKERRQ(info);
1089:   info = TaoOptionDouble("-tao_nls_alpha1", "", "", ls->alpha1, &ls->alpha1, 0); CHKERRQ(info);
1090:   info = TaoOptionDouble("-tao_nls_alpha2", "", "", ls->alpha2, &ls->alpha2, 0); CHKERRQ(info);
1091:   info = TaoOptionDouble("-tao_nls_alpha3", "", "", ls->alpha3, &ls->alpha3, 0); CHKERRQ(info);
1092:   info = TaoOptionDouble("-tao_nls_alpha4", "", "", ls->alpha4, &ls->alpha4, 0); CHKERRQ(info);
1093:   info = TaoOptionDouble("-tao_nls_alpha5", "", "", ls->alpha5, &ls->alpha5, 0); CHKERRQ(info);
1094:   info = TaoOptionDouble("-tao_nls_nu1", "poor steplength; reduce radius", "", ls->nu1, &ls->nu1, 0); CHKERRQ(info);
1095:   info = TaoOptionDouble("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", ls->nu2, &ls->nu2, 0); CHKERRQ(info);
1096:   info = TaoOptionDouble("-tao_nls_nu3", "good steplength; increase radius", "", ls->nu3, &ls->nu3, 0); CHKERRQ(info);
1097:   info = TaoOptionDouble("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", ls->nu4, &ls->nu4, 0); CHKERRQ(info);
1098:   info = TaoOptionDouble("-tao_nls_omega1", "", "", ls->omega1, &ls->omega1, 0); CHKERRQ(info);
1099:   info = TaoOptionDouble("-tao_nls_omega2", "", "", ls->omega2, &ls->omega2, 0); CHKERRQ(info);
1100:   info = TaoOptionDouble("-tao_nls_omega3", "", "", ls->omega3, &ls->omega3, 0); CHKERRQ(info);
1101:   info = TaoOptionDouble("-tao_nls_omega4", "", "", ls->omega4, &ls->omega4, 0); CHKERRQ(info);
1102:   info = TaoOptionDouble("-tao_nls_omega5", "", "", ls->omega5, &ls->omega5, 0); CHKERRQ(info);
1103:   info = TaoOptionDouble("-tao_nls_mu1_i", "", "", ls->mu1_i, &ls->mu1_i, 0); CHKERRQ(info);
1104:   info = TaoOptionDouble("-tao_nls_mu2_i", "", "", ls->mu2_i, &ls->mu2_i, 0); CHKERRQ(info);
1105:   info = TaoOptionDouble("-tao_nls_gamma1_i", "", "", ls->gamma1_i, &ls->gamma1_i, 0); CHKERRQ(info);
1106:   info = TaoOptionDouble("-tao_nls_gamma2_i", "", "", ls->gamma2_i, &ls->gamma2_i, 0); CHKERRQ(info);
1107:   info = TaoOptionDouble("-tao_nls_gamma3_i", "", "", ls->gamma3_i, &ls->gamma3_i, 0); CHKERRQ(info);
1108:   info = TaoOptionDouble("-tao_nls_gamma4_i", "", "", ls->gamma4_i, &ls->gamma4_i, 0); CHKERRQ(info);
1109:   info = TaoOptionDouble("-tao_nls_theta_i", "", "", ls->theta_i, &ls->theta_i, 0); CHKERRQ(info);
1110:   info = TaoOptionDouble("-tao_nls_mu1", "", "", ls->mu1, &ls->mu1, 0); CHKERRQ(info);
1111:   info = TaoOptionDouble("-tao_nls_mu2", "", "", ls->mu2, &ls->mu2, 0); CHKERRQ(info);
1112:   info = TaoOptionDouble("-tao_nls_gamma1", "", "", ls->gamma1, &ls->gamma1, 0); CHKERRQ(info);
1113:   info = TaoOptionDouble("-tao_nls_gamma2", "", "", ls->gamma2, &ls->gamma2, 0); CHKERRQ(info);
1114:   info = TaoOptionDouble("-tao_nls_gamma3", "", "", ls->gamma3, &ls->gamma3, 0); CHKERRQ(info);
1115:   info = TaoOptionDouble("-tao_nls_gamma4", "", "", ls->gamma4, &ls->gamma4, 0); CHKERRQ(info);
1116:   info = TaoOptionDouble("-tao_nls_theta", "", "", ls->theta, &ls->theta, 0); CHKERRQ(info);
1117:   info = TaoOptionDouble("-tao_nls_min_radius", "lower bound on initial radius", "", ls->min_radius, &ls->min_radius, 0); CHKERRQ(info);
1118:   info = TaoOptionDouble("-tao_nls_max_radius", "upper bound on radius", "", ls->max_radius, &ls->max_radius, 0); CHKERRQ(info);
1119:   info = TaoOptionDouble("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", ls->epsilon, &ls->epsilon, 0); CHKERRQ(info);
1120:   info = TaoLineSearchSetFromOptions(tao); CHKERRQ(info);
1121:   info = TaoOptionsTail(); CHKERRQ(info);
1122:   TaoFunctionReturn(0);
1123: }


1126: /*------------------------------------------------------------*/
1129: static int TaoView_NLS(TAO_SOLVER tao,void* solver)
1130: {
1131:   TAO_NLS *ls = (TAO_NLS *)solver;
1132:   int info;

1134:   TaoFunctionBegin;
1135:   if (NLS_PC_BFGS == ls->pc_type && ls->M) {
1136:     info = TaoPrintInt(tao, "  Rejected matrix updates: %d\n", ls->M->GetRejects()); CHKERRQ(info);
1137:   }
1138:   info = TaoPrintInt(tao, "  Newton steps: %d\n", ls->newt); CHKERRQ(info);
1139:   info = TaoPrintInt(tao, "  BFGS steps: %d\n", ls->bfgs); CHKERRQ(info);
1140:   info = TaoPrintInt(tao, "  Scaled gradient steps: %d\n", ls->sgrad); CHKERRQ(info);
1141:   info = TaoPrintInt(tao, "  Gradient steps: %d\n", ls->grad); CHKERRQ(info);

1143:   info = TaoPrintInt(tao, "  nls ksp atol: %d\n", ls->ksp_atol); CHKERRQ(info);
1144:   info = TaoPrintInt(tao, "  nls ksp rtol: %d\n", ls->ksp_rtol); CHKERRQ(info);
1145:   info = TaoPrintInt(tao, "  nls ksp ctol: %d\n", ls->ksp_ctol); CHKERRQ(info);
1146:   info = TaoPrintInt(tao, "  nls ksp negc: %d\n", ls->ksp_negc); CHKERRQ(info);
1147:   info = TaoPrintInt(tao, "  nls ksp dtol: %d\n", ls->ksp_dtol); CHKERRQ(info);
1148:   info = TaoPrintInt(tao, "  nls ksp iter: %d\n", ls->ksp_iter); CHKERRQ(info);
1149:   info = TaoPrintInt(tao, "  nls ksp othr: %d\n", ls->ksp_othr); CHKERRQ(info);
1150:   info = TaoLineSearchView(tao); CHKERRQ(info);
1151:   TaoFunctionReturn(0);
1152: }

1154: /* ---------------------------------------------------------- */
1158: int TaoCreate_NLS(TAO_SOLVER tao)
1159: {
1160:   TAO_NLS *ls;
1161:   int info;

1163:   TaoFunctionBegin;

1165:   info = TaoNew(TAO_NLS, &ls); CHKERRQ(info);
1166:   info = PetscLogObjectMemory(tao, sizeof(TAO_NLS)); CHKERRQ(info);

1168:   info = TaoSetTaoSolveRoutine(tao, TaoSolve_NLS, (void *)ls); CHKERRQ(info);
1169:   info = TaoSetTaoSetUpDownRoutines(tao, TaoSetUp_NLS, TaoDestroy_NLS); CHKERRQ(info);
1170:   info = TaoSetTaoOptionsRoutine(tao, TaoSetOptions_NLS); CHKERRQ(info);
1171:   info = TaoSetTaoViewRoutine(tao, TaoView_NLS); CHKERRQ(info);

1173:   info = TaoSetMaximumIterates(tao, 50); CHKERRQ(info);
1174:   info = TaoSetTolerances(tao, 1e-10, 1e-10, 0, 0); CHKERRQ(info);

1176:   info = TaoSetTrustRegionRadius(tao, 100.0); CHKERRQ(info);
1177:   info = TaoSetTrustRegionTolerance(tao, 1.0e-12); CHKERRQ(info);

1179:   ls->sval   = 0.0;
1180:   ls->imin   = 1.0e-4;
1181:   ls->imax   = 1.0e+2;
1182:   ls->imfac  = 1.0e-1;

1184:   ls->pmin   = 1.0e-12;
1185:   ls->pmax   = 1.0e+2;
1186:   ls->pgfac  = 1.0e+1;
1187:   ls->psfac  = 4.0e-1;
1188:   ls->pmgfac = 1.0e-1;
1189:   ls->pmsfac = 1.0e-1;

1191:   // Default values for trust-region radius update based on steplength
1192:   ls->nu1 = 0.25;
1193:   ls->nu2 = 0.50;
1194:   ls->nu3 = 1.00;
1195:   ls->nu4 = 1.25;

1197:   ls->omega1 = 0.25;
1198:   ls->omega2 = 0.50;
1199:   ls->omega3 = 1.00;
1200:   ls->omega4 = 2.00;
1201:   ls->omega5 = 4.00;

1203:   // Default values for trust-region radius update based on reduction
1204:   ls->eta1 = 1.0e-4;
1205:   ls->eta2 = 0.25;
1206:   ls->eta3 = 0.50;
1207:   ls->eta4 = 0.90;

1209:   ls->alpha1 = 0.25;
1210:   ls->alpha2 = 0.50;
1211:   ls->alpha3 = 1.00;
1212:   ls->alpha4 = 2.00;
1213:   ls->alpha5 = 4.00;

1215:   // Default values for trust-region radius update based on interpolation
1216:   ls->mu1 = 0.10;
1217:   ls->mu2 = 0.50;

1219:   ls->gamma1 = 0.25;
1220:   ls->gamma2 = 0.50;
1221:   ls->gamma3 = 2.00;
1222:   ls->gamma4 = 4.00;

1224:   ls->theta = 0.05;

1226:   // Default values for trust region initialization based on interpolation
1227:   ls->mu1_i = 0.35;
1228:   ls->mu2_i = 0.50;

1230:   ls->gamma1_i = 0.0625;
1231:   ls->gamma2_i = 0.5;
1232:   ls->gamma3_i = 2.0;
1233:   ls->gamma4_i = 5.0;
1234:   
1235:   ls->theta_i = 0.25;

1237:   // Remaining parameters
1238:   ls->min_radius = 1.0e-10;
1239:   ls->max_radius = 1.0e10;
1240:   ls->epsilon = 1.0e-6;

1242:   ls->ksp_type        = NLS_KSP_STCG;
1243:   ls->pc_type         = NLS_PC_BFGS;
1244:   ls->bfgs_scale_type = BFGS_SCALE_PHESS;
1245:   ls->init_type       = NLS_INIT_INTERPOLATION;
1246:   ls->update_type     = NLS_UPDATE_STEP;

1248:   info = TaoCreateMoreThuenteLineSearch(tao, 0, 0); CHKERRQ(info);
1249:   TaoFunctionReturn(0);
1250: }

1253: #endif