Actual source code: nls.c

  1: #include "taolinesearch.h"
  2: #include "src/matrix/lmvmmat.h"
  3: #include "nls.h"

  5: #include "petscksp.h"
  6: #include "petscpc.h"
  7: #include "private/kspimpl.h"
  8: #include "private/pcimpl.h"

 10: #define NLS_KSP_CG        0
 11: #define NLS_KSP_NASH        1
 12: #define NLS_KSP_STCG        2
 13: #define NLS_KSP_GLTR        3
 14: #define NLS_KSP_PETSC        4
 15: #define NLS_KSP_TYPES        5

 17: #define NLS_PC_NONE        0
 18: #define NLS_PC_AHESS        1
 19: #define NLS_PC_BFGS        2
 20: #define NLS_PC_PETSC        3
 21: #define NLS_PC_TYPES        4

 23: #define BFGS_SCALE_AHESS        0
 24: #define BFGS_SCALE_PHESS        1
 25: #define BFGS_SCALE_BFGS                2
 26: #define BFGS_SCALE_TYPES        3

 28: #define NLS_INIT_CONSTANT         0
 29: #define NLS_INIT_DIRECTION        1
 30: #define NLS_INIT_INTERPOLATION    2
 31: #define NLS_INIT_TYPES            3

 33: #define NLS_UPDATE_STEP           0
 34: #define NLS_UPDATE_REDUCTION      1
 35: #define NLS_UPDATE_INTERPOLATION  2
 36: #define NLS_UPDATE_TYPES          3

 38: static const char *NLS_KSP[64] = {
 39:   "cg", "nash", "stcg", "gltr", "petsc"
 40: };

 42: static const char *NLS_PC[64] = {
 43:   "none", "ahess", "bfgs", "petsc"
 44: };

 46: static const char *BFGS_SCALE[64] = {
 47:   "ahess", "phess", "bfgs"
 48: };

 50: static const char *NLS_INIT[64] = {
 51:   "constant", "direction", "interpolation"
 52: };

 54: static const char *NLS_UPDATE[64] = {
 55:   "step", "reduction", "interpolation"
 56: };

 58: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x) ;
 59: /* Routine for BFGS preconditioner 


 62:  Implements Newton's Method with a line search approach for solving
 63:  unconstrained minimization problems.  A More'-Thuente line search 
 64:  is used to guarantee that the bfgs preconditioner remains positive
 65:  definite.

 67:  The method can shift the Hessian matrix.  The shifting procedure is
 68:  adapted from the PATH algorithm for solving complementarity
 69:  problems.

 71:  The linear system solve should be done with a conjugate gradient
 72:  method, although any method can be used. */

 74: #define NLS_NEWTON                 0
 75: #define NLS_BFGS                 1
 76: #define NLS_SCALED_GRADIENT         2
 77: #define NLS_GRADIENT                 3

 81: static PetscErrorCode TaoSolve_NLS(TaoSolver tao)
 82: {
 84:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

 86:   PC pc;

 88:   KSPConvergedReason ksp_reason;
 89:   TaoLineSearchTerminationReason ls_reason;
 90:   TaoSolverTerminationReason reason;
 91:   
 92:   PetscReal fmin, ftrial, f_full, prered, actred, kappa, sigma;
 93:   PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
 94:   PetscReal f, fold, gdx, gnorm, pert;
 95:   PetscReal step = 1.0;

 97:   PetscReal delta;
 98:   PetscReal norm_d = 0.0, e_min;

100:   MatStructure matflag;

102:   PetscInt stepType;
103:   PetscInt iter = 0;
104:   PetscInt bfgsUpdates = 0;
105:   PetscInt n,N,kspits;
106:   PetscInt needH;
107:   
108:   PetscInt i_max = 5;
109:   PetscInt j_max = 1;
110:   PetscInt i, j;


114:   if (tao->XL || tao->XU || tao->ops->computebounds) {
115:     PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by nls algorithm\n"); 
116:   }

118:   /* Initialized variables */
119:   pert = nlsP->sval;

121:   nlsP->ksp_atol = 0;
122:   nlsP->ksp_rtol = 0;
123:   nlsP->ksp_dtol = 0;
124:   nlsP->ksp_ctol = 0;
125:   nlsP->ksp_negc = 0;
126:   nlsP->ksp_iter = 0;
127:   nlsP->ksp_othr = 0;



131:   /* Modify the linear solver to a trust region method if desired */

133:   switch(nlsP->ksp_type) {
134:   case NLS_KSP_CG:
135:     KSPSetType(tao->ksp, KSPCG); 
136:     if (tao->ksp->ops->setfromoptions) {
137:         (*tao->ksp->ops->setfromoptions)(tao->ksp);
138:     }
139:     break;

141:   case NLS_KSP_NASH:
142:     KSPSetType(tao->ksp, KSPNASH); 
143:     if (tao->ksp->ops->setfromoptions) {
144:       (*tao->ksp->ops->setfromoptions)(tao->ksp);
145:     }
146:     break;

148:   case NLS_KSP_STCG:
149:     KSPSetType(tao->ksp, KSPSTCG); 
150:     if (tao->ksp->ops->setfromoptions) {
151:       (*tao->ksp->ops->setfromoptions)(tao->ksp);
152:     }
153:     break;

155:   case NLS_KSP_GLTR:
156:     KSPSetType(tao->ksp, KSPGLTR); 
157:     if (tao->ksp->ops->setfromoptions) {
158:       (*tao->ksp->ops->setfromoptions)(tao->ksp);
159:     }
160:     break;

162:   default:
163:     /* Use the method set by the ksp_type */
164:     break;
165:   }

167:   /* Initialize trust-region radius when using nash, stcg, or gltr
168:      Will be reset during the first iteration */
169:   if (NLS_KSP_NASH == nlsP->ksp_type) {
170:       KSPNASHSetRadius(tao->ksp,nlsP->max_radius); 
171:   } else if (NLS_KSP_STCG == nlsP->ksp_type) {
172:       KSPSTCGSetRadius(tao->ksp,nlsP->max_radius); 
173:   } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
174:       KSPGLTRSetRadius(tao->ksp,nlsP->max_radius); 
175:   }
176:   
177:   
178:   if (NLS_KSP_NASH == nlsP->ksp_type ||
179:       NLS_KSP_STCG == nlsP->ksp_type || 
180:       NLS_KSP_GLTR == nlsP->ksp_type) {
181:     tao->trust = tao->trust0;

183:     if (tao->trust < 0.0) {
184:       SETERRQ(PETSC_COMM_SELF,1, "Initial radius negative");
185:     }

187:     /* Modify the radius if it is too large or small */
188:     tao->trust = PetscMax(tao->trust, nlsP->min_radius);
189:     tao->trust = PetscMin(tao->trust, nlsP->max_radius);
190:   }

192:   /* Get vectors we will need */

194:   if (NLS_PC_BFGS == nlsP->pc_type && !nlsP->M) {
195:     VecGetLocalSize(tao->solution,&n); 
196:     VecGetSize(tao->solution,&N); 
197:     MatCreateLMVM(((PetscObject)tao)->comm,n,N,&nlsP->M); 
198:     MatLMVMAllocateVectors(nlsP->M,tao->solution); 
199:   }

201:   /* Check convergence criteria */
202:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient); 
203:   VecNorm(tao->gradient,NORM_2,&gnorm); 
204:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
205:     SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
206:   }
207:   needH = 1;

209:   TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); 
210:   if (reason != TAO_CONTINUE_ITERATING) {
211:     return(0);
212:   }

214:   /* create vectors for the limited memory preconditioner */
215:   if ((NLS_PC_BFGS == nlsP->pc_type) && 
216:       (BFGS_SCALE_BFGS != nlsP->bfgs_scale_type)) {
217:     if (!nlsP->Diag) {
218:         VecDuplicate(tao->solution,&nlsP->Diag); 
219:     }
220:   }


223:   /* Modify the preconditioner to use the bfgs approximation */
224:   KSPGetPC(tao->ksp, &pc); 
225:   switch(nlsP->pc_type) {
226:   case NLS_PC_NONE:
227:     PCSetType(pc, PCNONE); 
228:     if (pc->ops->setfromoptions) {
229:       (*pc->ops->setfromoptions)(pc);
230:     }
231:     break;

233:   case NLS_PC_AHESS:
234:     PCSetType(pc, PCJACOBI); 
235:     if (pc->ops->setfromoptions) {
236:       (*pc->ops->setfromoptions)(pc);
237:     }
238:     PCJacobiSetUseAbs(pc); 
239:     break;

241:   case NLS_PC_BFGS:
242:     PCSetType(pc, PCSHELL); 
243:     if (pc->ops->setfromoptions) {
244:       (*pc->ops->setfromoptions)(pc);
245:     }
246:     PCShellSetName(pc, "bfgs"); 
247:     PCShellSetContext(pc, nlsP->M); 
248:     PCShellSetApply(pc, MatLMVMSolveShell); 
249:     break;

251:   default:
252:     /* Use the pc method set by pc_type */
253:     break;
254:   }

256:   /* Initialize trust-region radius.  The initialization is only performed 
257:      when we are using Nash, Steihaug-Toint or the Generalized Lanczos method. */
258:   if (NLS_KSP_NASH == nlsP->ksp_type ||
259:       NLS_KSP_STCG == nlsP->ksp_type || 
260:       NLS_KSP_GLTR == nlsP->ksp_type) {
261:     switch(nlsP->init_type) {
262:     case NLS_INIT_CONSTANT:
263:       /* Use the initial radius specified */
264:       break;

266:     case NLS_INIT_INTERPOLATION:
267:       /* Use the initial radius specified */
268:       max_radius = 0.0;
269:   
270:       for (j = 0; j < j_max; ++j) {
271:         fmin = f;
272:         sigma = 0.0;
273:   
274:         if (needH) {
275:             TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag); 
276:             needH = 0;
277:         }
278:   
279:         for (i = 0; i < i_max; ++i) {
280:           VecCopy(tao->solution,nlsP->W); 
281:           VecAXPY(nlsP->W,-tao->trust/gnorm,tao->gradient); 
282:           TaoComputeObjective(tao, nlsP->W, &ftrial); 
283:           if (PetscIsInfOrNanReal(ftrial)) {
284:             tau = nlsP->gamma1_i;
285:           }
286:           else {
287:             if (ftrial < fmin) {
288:               fmin = ftrial;
289:               sigma = -tao->trust / gnorm;
290:             }
291:             
292:             MatMult(tao->hessian, tao->gradient, nlsP->D); 
293:             VecDot(tao->gradient, nlsP->D, &prered); 
294:   
295:             prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
296:             actred = f - ftrial;
297:             if ((PetscAbsScalar(actred) <= nlsP->epsilon) && 
298:                 (PetscAbsScalar(prered) <= nlsP->epsilon)) {
299:               kappa = 1.0;
300:             }
301:             else {
302:               kappa = actred / prered;
303:             }
304:   
305:             tau_1 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust + (1.0 - nlsP->theta_i) * prered - actred);
306:             tau_2 = nlsP->theta_i * gnorm * tao->trust / (nlsP->theta_i * gnorm * tao->trust - (1.0 + nlsP->theta_i) * prered + actred);
307:             tau_min = PetscMin(tau_1, tau_2);
308:             tau_max = PetscMax(tau_1, tau_2);
309:   
310:             if (PetscAbsScalar(kappa - 1.0) <= nlsP->mu1_i) {
311:               /* Great agreement */
312:               max_radius = PetscMax(max_radius, tao->trust);
313:   
314:               if (tau_max < 1.0) {
315:                 tau = nlsP->gamma3_i;
316:               }
317:               else if (tau_max > nlsP->gamma4_i) {
318:                 tau = nlsP->gamma4_i;
319:               }
320:               else if (tau_1 >= 1.0 && tau_1 <= nlsP->gamma4_i && tau_2 < 1.0) {
321:                 tau = tau_1;
322:               }
323:               else if (tau_2 >= 1.0 && tau_2 <= nlsP->gamma4_i && tau_1 < 1.0) {
324:                 tau = tau_2;
325:               }
326:               else {
327:                 tau = tau_max;
328:               }
329:             }
330:             else if (PetscAbsScalar(kappa - 1.0) <= nlsP->mu2_i) {
331:               /* Good agreement */
332:               max_radius = PetscMax(max_radius, tao->trust);
333:   
334:               if (tau_max < nlsP->gamma2_i) {
335:                 tau = nlsP->gamma2_i;
336:               }
337:               else if (tau_max > nlsP->gamma3_i) {
338:                 tau = nlsP->gamma3_i;
339:               }
340:               else {
341:                 tau = tau_max;
342:               }
343:             }
344:             else {
345:               /* Not good agreement */
346:               if (tau_min > 1.0) {
347:                 tau = nlsP->gamma2_i;
348:               }
349:               else if (tau_max < nlsP->gamma1_i) {
350:                 tau = nlsP->gamma1_i;
351:               }
352:               else if ((tau_min < nlsP->gamma1_i) && (tau_max >= 1.0)) {
353:                 tau = nlsP->gamma1_i;
354:               }
355:               else if ((tau_1 >= nlsP->gamma1_i) && (tau_1 < 1.0) &&
356:                        ((tau_2 < nlsP->gamma1_i) || (tau_2 >= 1.0))) {
357:                 tau = tau_1;
358:               }
359:               else if ((tau_2 >= nlsP->gamma1_i) && (tau_2 < 1.0) &&
360:                        ((tau_1 < nlsP->gamma1_i) || (tau_2 >= 1.0))) {
361:                 tau = tau_2;
362:               }
363:               else {
364:                 tau = tau_max;
365:               }
366:             }
367:           }
368:           tao->trust = tau * tao->trust;
369:         }
370:   
371:         if (fmin < f) {
372:           f = fmin;
373:           VecAXPY(tao->solution,sigma,tao->gradient); 
374:           TaoComputeGradient(tao,tao->solution,tao->gradient); 

376:           VecNorm(tao->gradient,NORM_2,&gnorm); 
377:           if (PetscIsInfOrNanReal(gnorm)) {
378:             SETERRQ(PETSC_COMM_SELF,1, "User provided compute gradient generated Inf or NaN");
379:           }
380:           needH = 1;
381:   
382:           TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason); 
383:           if (reason != TAO_CONTINUE_ITERATING) {
384:             return(0);
385:           }
386:         }
387:       }
388:       tao->trust = PetscMax(tao->trust, max_radius);

390:       /* Modify the radius if it is too large or small */
391:       tao->trust = PetscMax(tao->trust, nlsP->min_radius);
392:       tao->trust = PetscMin(tao->trust, nlsP->max_radius);
393:       break;

395:     default:
396:       /* Norm of the first direction will initialize radius */
397:       tao->trust = 0.0;
398:       break;
399:     }
400:   } 

402:   /* Set initial scaling for the BFGS preconditioner
403:      This step is done after computing the initial trust-region radius
404:      since the function value may have decreased */
405:   if (NLS_PC_BFGS == nlsP->pc_type) {
406:     if (f != 0.0) {
407:       delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
408:     }
409:     else {
410:       delta = 2.0 / (gnorm*gnorm);
411:     }
412:     MatLMVMSetDelta(nlsP->M,delta); 
413:   }

415:   /* Set counter for gradient/reset steps*/
416:   nlsP->newt = 0;
417:   nlsP->bfgs = 0;
418:   nlsP->sgrad = 0;
419:   nlsP->grad = 0;

421:   /* Have not converged; continue with Newton method */
422:   while (reason == TAO_CONTINUE_ITERATING) {
423:     ++iter;

425:     /* Compute the Hessian */
426:     if (needH) {
427:         TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag); 
428:       needH = 0;
429:     }

431:     if ((NLS_PC_BFGS == nlsP->pc_type) && 
432:         (BFGS_SCALE_AHESS == nlsP->bfgs_scale_type)) {
433:       /* Obtain diagonal for the bfgs preconditioner  */
434:       MatGetDiagonal(tao->hessian, nlsP->Diag); 
435:       VecAbs(nlsP->Diag); 
436:       VecReciprocal(nlsP->Diag); 
437:       MatLMVMSetScale(nlsP->M,nlsP->Diag); 
438:     }
439:  
440:     /* Shift the Hessian matrix */
441:     if (pert > 0) {
442:       MatShift(tao->hessian, pert);
443:       if (tao->hessian != tao->hessian_pre) {
444:         MatShift(tao->hessian_pre, pert); 
445:       }
446:     }

448:     
449:     if (NLS_PC_BFGS == nlsP->pc_type) {
450:       if (BFGS_SCALE_PHESS == nlsP->bfgs_scale_type) {
451:         /* Obtain diagonal for the bfgs preconditioner  */
452:           MatGetDiagonal(tao->hessian, nlsP->Diag); 
453:           VecAbs(nlsP->Diag); 
454:           VecReciprocal(nlsP->Diag); 
455:           MatLMVMSetScale(nlsP->M,nlsP->Diag); 
456:       }
457:       /* Update the limited memory preconditioner */
458:       MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient); 
459:       ++bfgsUpdates;
460:     }

462:     /* Solve the Newton system of equations */
463:     KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre,matflag); 
464:     if (NLS_KSP_NASH == nlsP->ksp_type ||
465:         NLS_KSP_STCG == nlsP->ksp_type || 
466:         NLS_KSP_GLTR == nlsP->ksp_type) {

468:       if (NLS_KSP_NASH == nlsP->ksp_type) {
469:         KSPNASHSetRadius(tao->ksp,nlsP->max_radius); 
470:       } else if (NLS_KSP_STCG == nlsP->ksp_type) {
471:          KSPSTCGSetRadius(tao->ksp,nlsP->max_radius); 
472:       } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
473:         KSPGLTRSetRadius(tao->ksp,nlsP->max_radius); 
474:       }
475:         
476:       KSPSolve(tao->ksp, tao->gradient, nlsP->D); 
477:       KSPGetIterationNumber(tao->ksp,&kspits); 
478:       tao->ksp_its+=kspits;

480:       if (NLS_KSP_NASH == nlsP->ksp_type) {
481:         KSPNASHGetNormD(tao->ksp,&norm_d); 
482:       } else if (NLS_KSP_STCG == nlsP->ksp_type) {
483:          KSPSTCGGetNormD(tao->ksp,&norm_d); 
484:       } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
485:         KSPGLTRGetNormD(tao->ksp,&norm_d); 
486:       }

488:       if (0.0 == tao->trust) {
489:         /* Radius was uninitialized; use the norm of the direction */
490:         if (norm_d > 0.0) {
491:           tao->trust = norm_d;

493:           /* Modify the radius if it is too large or small */
494:           tao->trust = PetscMax(tao->trust, nlsP->min_radius);
495:           tao->trust = PetscMin(tao->trust, nlsP->max_radius);
496:         }
497:         else {
498:           /* The direction was bad; set radius to default value and re-solve 
499:              the trust-region subproblem to get a direction */
500:           tao->trust = tao->trust0;

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

506:           if (NLS_KSP_NASH == nlsP->ksp_type) {
507:             KSPNASHSetRadius(tao->ksp,nlsP->max_radius); 
508:           } else if (NLS_KSP_STCG == nlsP->ksp_type) {
509:             KSPSTCGSetRadius(tao->ksp,nlsP->max_radius); 
510:           } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
511:             KSPGLTRSetRadius(tao->ksp,nlsP->max_radius); 
512:           }
513:         
514:           KSPSolve(tao->ksp, tao->gradient, nlsP->D); 
515:           KSPGetIterationNumber(tao->ksp,&kspits); 
516:           tao->ksp_its+=kspits;
517:           if (NLS_KSP_NASH == nlsP->ksp_type) {
518:             KSPNASHGetNormD(tao->ksp,&norm_d); 
519:           } else if (NLS_KSP_STCG == nlsP->ksp_type) {
520:             KSPSTCGGetNormD(tao->ksp,&norm_d); 
521:           } else if (NLS_KSP_GLTR == nlsP->ksp_type) {
522:             KSPGLTRGetNormD(tao->ksp,&norm_d); 
523:           }

525:           if (norm_d == 0.0) {
526:             SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
527:           }
528:         }
529:       }
530:     }
531:     else {
532:       KSPSolve(tao->ksp, tao->gradient, nlsP->D); 
533:       KSPGetIterationNumber(tao->ksp, &kspits); 
534:       tao->ksp_its += kspits;
535:     }
536:     VecScale(nlsP->D, -1.0); 
537:     KSPGetConvergedReason(tao->ksp, &ksp_reason); 
538:     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
539:         (NLS_PC_BFGS == nlsP->pc_type) && (bfgsUpdates > 1)) {
540:       /* Preconditioner is numerically indefinite; reset the 
541:          approximate if using BFGS preconditioning. */

543:       if (f != 0.0) {
544:         delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
545:       }
546:       else {
547:         delta = 2.0 / (gnorm*gnorm);
548:       }
549:       MatLMVMSetDelta(nlsP->M,delta); 
550:       MatLMVMReset(nlsP->M); 
551:       MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient); 
552:       bfgsUpdates = 1;
553:     }

555:     if (KSP_CONVERGED_ATOL == ksp_reason) {
556:       ++nlsP->ksp_atol;
557:     }
558:     else if (KSP_CONVERGED_RTOL == ksp_reason) {
559:       ++nlsP->ksp_rtol;
560:     }
561:     else if (KSP_CONVERGED_CG_CONSTRAINED == ksp_reason) {
562:       ++nlsP->ksp_ctol;
563:     }
564:     else if (KSP_CONVERGED_CG_NEG_CURVE == ksp_reason) {
565:       ++nlsP->ksp_negc;
566:     }
567:     else if (KSP_DIVERGED_DTOL == ksp_reason) {
568:       ++nlsP->ksp_dtol;
569:     }
570:     else if (KSP_DIVERGED_ITS == ksp_reason) {
571:       ++nlsP->ksp_iter;
572:     }
573:     else {
574:       ++nlsP->ksp_othr;
575:     } 

577:     /* Check for success (descent direction) */
578:     VecDot(nlsP->D, tao->gradient, &gdx); 
579:     if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
580:       /* Newton step is not descent or direction produced Inf or NaN
581:          Update the perturbation for next time */
582:       if (pert <= 0.0) {
583:         /* Initialize the perturbation */
584:         pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
585:         if (NLS_KSP_GLTR == nlsP->ksp_type) {
586:           KSPGLTRGetMinEig(tao->ksp,&e_min); 
587:           pert = PetscMax(pert, -e_min);
588:         }
589:       }
590:       else {
591:         /* Increase the perturbation */
592:         pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
593:       }

595:       if (NLS_PC_BFGS != nlsP->pc_type) {
596:         /* We don't have the bfgs matrix around and updated
597:            Must use gradient direction in this case */
598:         VecCopy(tao->gradient, nlsP->D); 
599:         VecScale(nlsP->D, -1.0); 
600:         ++nlsP->grad;
601:         stepType = NLS_GRADIENT;
602:       }
603:       else {
604:         /* Attempt to use the BFGS direction */
605:         MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D); 
606:         VecScale(nlsP->D, -1.0); 

608:         /* Check for success (descent direction) */
609:         VecDot(tao->gradient, nlsP->D, &gdx); 
610:         if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
611:           /* BFGS direction is not descent or direction produced not a number
612:              We can assert bfgsUpdates > 1 in this case because
613:              the first solve produces the scaled gradient direction,
614:              which is guaranteed to be descent */

616:           /* Use steepest descent direction (scaled) */

618:           if (f != 0.0) {
619:             delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
620:           }
621:           else {
622:             delta = 2.0 / (gnorm*gnorm);
623:           }
624:           MatLMVMSetDelta(nlsP->M, delta); 
625:           MatLMVMReset(nlsP->M); 
626:           MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient); 
627:           MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D); 
628:           VecScale(nlsP->D, -1.0); 
629:   
630:           bfgsUpdates = 1;
631:           ++nlsP->sgrad;
632:           stepType = NLS_SCALED_GRADIENT;
633:         }
634:         else {
635:           if (1 == bfgsUpdates) {
636:             /* The first BFGS direction is always the scaled gradient */
637:             ++nlsP->sgrad;
638:             stepType = NLS_SCALED_GRADIENT;
639:           }
640:           else {
641:             ++nlsP->bfgs;
642:             stepType = NLS_BFGS;
643:           }
644:         }
645:       }
646:     }
647:     else {
648:       /* Computed Newton step is descent */
649:       switch (ksp_reason) {
650:       case KSP_DIVERGED_NAN:
651:       case KSP_DIVERGED_BREAKDOWN:
652:       case KSP_DIVERGED_INDEFINITE_MAT:
653:       case KSP_DIVERGED_INDEFINITE_PC:
654:       case KSP_CONVERGED_CG_NEG_CURVE:
655:         /* Matrix or preconditioner is indefinite; increase perturbation */
656:         if (pert <= 0.0) {
657:           /* Initialize the perturbation */
658:           pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
659:           if (NLS_KSP_GLTR == nlsP->ksp_type) {
660:             KSPGLTRGetMinEig(tao->ksp, &e_min); 
661:             pert = PetscMax(pert, -e_min);
662:           }
663:         }
664:         else {
665:           /* Increase the perturbation */
666:           pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
667:         }
668:         break;

670:       default:
671:         /* Newton step computation is good; decrease perturbation */
672:         pert = PetscMin(nlsP->psfac * pert, nlsP->pmsfac * gnorm);
673:         if (pert < nlsP->pmin) {
674:           pert = 0.0;
675:         }
676:         break; 
677:       }

679:       ++nlsP->newt;
680:       stepType = NLS_NEWTON;
681:     }

683:     /* Perform the linesearch */
684:     fold = f;
685:     VecCopy(tao->solution, nlsP->Xold); 
686:     VecCopy(tao->gradient, nlsP->Gold); 

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

691:     while (ls_reason < 0  && stepType != NLS_GRADIENT) {
692:       /* Linesearch failed */
693:       f = fold;
694:       VecCopy(nlsP->Xold, tao->solution); 
695:       VecCopy(nlsP->Gold, tao->gradient); 

697:       switch(stepType) {
698:       case NLS_NEWTON:
699:         /* Failed to obtain acceptable iterate with Newton 1step
700:            Update the perturbation for next time */
701:         if (pert <= 0.0) {
702:           /* Initialize the perturbation */
703:           pert = PetscMin(nlsP->imax, PetscMax(nlsP->imin, nlsP->imfac * gnorm));
704:           if (NLS_KSP_GLTR == nlsP->ksp_type) {
705:             KSPGLTRGetMinEig(tao->ksp,&e_min); 
706:             pert = PetscMax(pert, -e_min);
707:           }
708:         }
709:         else {
710:           /* Increase the perturbation */
711:           pert = PetscMin(nlsP->pmax, PetscMax(nlsP->pgfac * pert, nlsP->pmgfac * gnorm));
712:         }

714:         if (NLS_PC_BFGS != nlsP->pc_type) {
715:           /* We don't have the bfgs matrix around and being updated
716:              Must use gradient direction in this case */
717:           VecCopy(tao->gradient, nlsP->D); 
718:           ++nlsP->grad;
719:           stepType = NLS_GRADIENT;
720:         }
721:         else {
722:           /* Attempt to use the BFGS direction */
723:           MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D); 
724:           /* Check for success (descent direction) */
725:           VecDot(tao->solution, nlsP->D, &gdx); 
726:           if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
727:             /* BFGS direction is not descent or direction produced not a number
728:                We can assert bfgsUpdates > 1 in this case
729:                Use steepest descent direction (scaled) */
730:     
731:             if (f != 0.0) {
732:               delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
733:             }
734:             else {
735:               delta = 2.0 / (gnorm*gnorm);
736:             }
737:             MatLMVMSetDelta(nlsP->M, delta); 
738:             MatLMVMReset(nlsP->M); 
739:             MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient); 
740:             MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D); 
741:   
742:             bfgsUpdates = 1;
743:             ++nlsP->sgrad;
744:             stepType = NLS_SCALED_GRADIENT;
745:           }
746:           else {
747:             if (1 == bfgsUpdates) {
748:               /* The first BFGS direction is always the scaled gradient */
749:               ++nlsP->sgrad;
750:               stepType = NLS_SCALED_GRADIENT;
751:             }
752:             else {
753:               ++nlsP->bfgs;
754:               stepType = NLS_BFGS;
755:             }
756:           }
757:         }
758:         break;

760:       case NLS_BFGS:
761:         /* Can only enter if pc_type == NLS_PC_BFGS
762:            Failed to obtain acceptable iterate with BFGS step
763:            Attempt to use the scaled gradient direction */

765:         if (f != 0.0) {
766:           delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
767:         }
768:         else {
769:           delta = 2.0 / (gnorm*gnorm);
770:         }
771:         MatLMVMSetDelta(nlsP->M, delta); 
772:         MatLMVMReset(nlsP->M); 
773:         MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient); 
774:         MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D); 

776:         bfgsUpdates = 1;
777:         ++nlsP->sgrad;
778:         stepType = NLS_SCALED_GRADIENT;
779:         break;

781:       case NLS_SCALED_GRADIENT:
782:         /* Can only enter if pc_type == NLS_PC_BFGS
783:            The scaled gradient step did not produce a new iterate;
784:            attemp to use the gradient direction.
785:            Need to make sure we are not using a different diagonal scaling */
786:         
787:         MatLMVMSetScale(nlsP->M,0); 
788:         MatLMVMSetDelta(nlsP->M,1.0); 
789:         MatLMVMReset(nlsP->M); 
790:         MatLMVMUpdate(nlsP->M, tao->solution, tao->gradient); 
791:         MatLMVMSolve(nlsP->M, tao->gradient, nlsP->D); 
792:         
793:         bfgsUpdates = 1;
794:         ++nlsP->grad;
795:         stepType = NLS_GRADIENT;
796:         break;
797:       }
798:       VecScale(nlsP->D, -1.0); 

800:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, nlsP->D, &step, &ls_reason); 
801:       TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full); 
802:       TaoAddLineSearchCounts(tao); 
803:     }

805:     if (ls_reason < 0) {
806:       /* Failed to find an improving point */
807:       f = fold;
808:       VecCopy(nlsP->Xold, tao->solution); 
809:       VecCopy(nlsP->Gold, tao->gradient); 
810:       step = 0.0;
811:       reason = TAO_DIVERGED_LS_FAILURE;
812:       tao->reason = TAO_DIVERGED_LS_FAILURE;
813:       break;
814:     }

816:     /* Update trust region radius */
817:     if (NLS_KSP_NASH == nlsP->ksp_type ||
818:         NLS_KSP_STCG == nlsP->ksp_type || 
819:         NLS_KSP_GLTR == nlsP->ksp_type) {
820:       switch(nlsP->update_type) {
821:       case NLS_UPDATE_STEP:
822:         if (stepType == NLS_NEWTON) {
823:           if (step < nlsP->nu1) {
824:             /* Very bad step taken; reduce radius */
825:             tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
826:           }
827:           else if (step < nlsP->nu2) {
828:             /* Reasonably bad step taken; reduce radius */
829:             tao->trust = nlsP->omega2 * PetscMin(norm_d, tao->trust);
830:           }
831:           else if (step < nlsP->nu3) {
832:             /*  Reasonable step was taken; leave radius alone */
833:             if (nlsP->omega3 < 1.0) {
834:               tao->trust = nlsP->omega3 * PetscMin(norm_d, tao->trust);
835:             }
836:             else if (nlsP->omega3 > 1.0) {
837:               tao->trust = PetscMax(nlsP->omega3 * norm_d, tao->trust);  
838:             }
839:           }
840:           else if (step < nlsP->nu4) {
841:             /*  Full step taken; increase the radius */
842:             tao->trust = PetscMax(nlsP->omega4 * norm_d, tao->trust);  
843:           }
844:           else {
845:             /*  More than full step taken; increase the radius */
846:             tao->trust = PetscMax(nlsP->omega5 * norm_d, tao->trust);  
847:           }
848:         }
849:         else {
850:           /*  Newton step was not good; reduce the radius */
851:           tao->trust = nlsP->omega1 * PetscMin(norm_d, tao->trust);
852:         }
853:         break;

855:       case NLS_UPDATE_REDUCTION:
856:         if (stepType == NLS_NEWTON) {
857:           /*  Get predicted reduction */

859:           if (NLS_KSP_STCG == nlsP->ksp_type) {
860:               KSPSTCGGetObjFcn(tao->ksp,&prered);
861:           } else if (NLS_KSP_NASH == nlsP->ksp_type)  {
862:               KSPNASHGetObjFcn(tao->ksp,&prered);
863:           } else {
864:               KSPGLTRGetObjFcn(tao->ksp,&prered); 
865:           }
866:           
867:               
868:           

870:           if (prered >= 0.0) {
871:             /*  The predicted reduction has the wrong sign.  This cannot */
872:             /*  happen in infinite precision arithmetic.  Step should */
873:             /*  be rejected! */
874:             tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
875:           }
876:           else {
877:             if (PetscIsInfOrNanReal(f_full)) {
878:               tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
879:             }
880:             else {
881:               /*  Compute and actual reduction */
882:               actred = fold - f_full;
883:               prered = -prered;
884:               if ((PetscAbsScalar(actred) <= nlsP->epsilon) && 
885:                   (PetscAbsScalar(prered) <= nlsP->epsilon)) {
886:                 kappa = 1.0;
887:               }
888:               else {
889:                 kappa = actred / prered;
890:               }
891:   
892:               /*  Accept of reject the step and update radius */
893:               if (kappa < nlsP->eta1) {
894:                 /*  Very bad step */
895:                 tao->trust = nlsP->alpha1 * PetscMin(tao->trust, norm_d);
896:               }
897:               else if (kappa < nlsP->eta2) {
898:                 /*  Marginal bad step */
899:                 tao->trust = nlsP->alpha2 * PetscMin(tao->trust, norm_d);
900:               }
901:               else if (kappa < nlsP->eta3) {
902:                 /*  Reasonable step */
903:                 if (nlsP->alpha3 < 1.0) {
904:                   tao->trust = nlsP->alpha3 * PetscMin(norm_d, tao->trust);
905:                 }
906:                 else if (nlsP->alpha3 > 1.0) {
907:                   tao->trust = PetscMax(nlsP->alpha3 * norm_d, tao->trust);  
908:                 }
909:               }
910:               else if (kappa < nlsP->eta4) {
911:                 /*  Good step */
912:                 tao->trust = PetscMax(nlsP->alpha4 * norm_d, tao->trust);
913:               }
914:               else {
915:                 /*  Very good step */
916:                 tao->trust = PetscMax(nlsP->alpha5 * norm_d, tao->trust);
917:               }
918:             }
919:           }
920:         }
921:         else {
922:           /*  Newton step was not good; reduce the radius */
923:           tao->trust = nlsP->alpha1 * PetscMin(norm_d, tao->trust);
924:         }
925:         break;

927:       default:
928:         if (stepType == NLS_NEWTON) {

930:           if (NLS_KSP_STCG == nlsP->ksp_type) {
931:               KSPSTCGGetObjFcn(tao->ksp,&prered);
932:           } else if (NLS_KSP_NASH == nlsP->ksp_type)  {
933:               KSPNASHGetObjFcn(tao->ksp,&prered);
934:           } else {
935:               KSPGLTRGetObjFcn(tao->ksp,&prered); 
936:           }
937:           if (prered >= 0.0) {
938:             /*  The predicted reduction has the wrong sign.  This cannot */
939:             /*  happen in infinite precision arithmetic.  Step should */
940:             /*  be rejected! */
941:             tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
942:           }
943:           else {
944:             if (PetscIsInfOrNanReal(f_full)) {
945:               tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
946:             }
947:             else {
948:               actred = fold - f_full;
949:               prered = -prered;
950:               if ((PetscAbsScalar(actred) <= nlsP->epsilon) && 
951:                   (PetscAbsScalar(prered) <= nlsP->epsilon)) {
952:                 kappa = 1.0;
953:               }
954:               else {
955:                 kappa = actred / prered;
956:               }

958:               tau_1 = nlsP->theta * gdx / (nlsP->theta * gdx - (1.0 - nlsP->theta) * prered + actred);
959:               tau_2 = nlsP->theta * gdx / (nlsP->theta * gdx + (1.0 + nlsP->theta) * prered - actred);
960:               tau_min = PetscMin(tau_1, tau_2);
961:               tau_max = PetscMax(tau_1, tau_2);

963:               if (kappa >= 1.0 - nlsP->mu1) {
964:                 /*  Great agreement */
965:                 if (tau_max < 1.0) {
966:                   tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
967:                 }
968:                 else if (tau_max > nlsP->gamma4) {
969:                   tao->trust = PetscMax(tao->trust, nlsP->gamma4 * norm_d);
970:                 }
971:                 else {
972:                   tao->trust = PetscMax(tao->trust, tau_max * norm_d);
973:                 }
974:               }
975:               else if (kappa >= 1.0 - nlsP->mu2) {
976:                 /*  Good agreement */

978:                 if (tau_max < nlsP->gamma2) {
979:                   tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
980:                 }
981:                 else if (tau_max > nlsP->gamma3) {
982:                   tao->trust = PetscMax(tao->trust, nlsP->gamma3 * norm_d);
983:                 }
984:                 else if (tau_max < 1.0) {
985:                   tao->trust = tau_max * PetscMin(tao->trust, norm_d);
986:                 }
987:                 else {
988:                   tao->trust = PetscMax(tao->trust, tau_max * norm_d);
989:                 }
990:               }
991:               else {
992:                 /*  Not good agreement */
993:                 if (tau_min > 1.0) {
994:                   tao->trust = nlsP->gamma2 * PetscMin(tao->trust, norm_d);
995:                 }
996:                 else if (tau_max < nlsP->gamma1) {
997:                   tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
998:                 }
999:                 else if ((tau_min < nlsP->gamma1) && (tau_max >= 1.0)) {
1000:                   tao->trust = nlsP->gamma1 * PetscMin(tao->trust, norm_d);
1001:                 }
1002:                 else if ((tau_1 >= nlsP->gamma1) && (tau_1 < 1.0) &&
1003:                          ((tau_2 < nlsP->gamma1) || (tau_2 >= 1.0))) {
1004:                   tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
1005:                 }
1006:                 else if ((tau_2 >= nlsP->gamma1) && (tau_2 < 1.0) &&
1007:                          ((tau_1 < nlsP->gamma1) || (tau_2 >= 1.0))) {
1008:                   tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
1009:                 }
1010:                 else {
1011:                   tao->trust = tau_max * PetscMin(tao->trust, norm_d);
1012:                 }
1013:               }
1014:             } 
1015:           }
1016:         }
1017:         else {
1018:           /*  Newton step was not good; reduce the radius */
1019:           tao->trust = nlsP->gamma1 * PetscMin(norm_d, tao->trust);
1020:         }
1021:         break;
1022:       }

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

1028:     /*  Check for termination */
1029:     VecNorm(tao->gradient, NORM_2, &gnorm); 
1030:     if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
1031:       SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number");
1032:     }
1033:     needH = 1;

1035:     TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason); 
1036:   }
1037:   return(0);
1038: }

1040: /* ---------------------------------------------------------- */
1043: static PetscErrorCode TaoSetUp_NLS(TaoSolver tao)
1044: {
1045:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;


1050:   if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);   }
1051:   if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection);   }
1052:   if (!nlsP->W) {VecDuplicate(tao->solution,&nlsP->W);   }
1053:   if (!nlsP->D) {VecDuplicate(tao->solution,&nlsP->D);   }
1054:   if (!nlsP->Xold) {VecDuplicate(tao->solution,&nlsP->Xold);   }
1055:   if (!nlsP->Gold) {VecDuplicate(tao->solution,&nlsP->Gold);   }

1057:   nlsP->Diag = 0;
1058:   nlsP->M = 0;

1060:   return(0);
1061: }

1063: /*------------------------------------------------------------*/
1066: static PetscErrorCode TaoDestroy_NLS(TaoSolver tao)
1067: {
1068:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

1072:   if (tao->setupcalled) {
1073:     VecDestroy(&nlsP->D); 
1074:     VecDestroy(&nlsP->W); 
1075:     VecDestroy(&nlsP->Xold); 
1076:     VecDestroy(&nlsP->Gold); 
1077:   }
1078:   if (nlsP->Diag) {
1079:     VecDestroy(&nlsP->Diag); 
1080:   }
1081:   if (nlsP->M) {
1082:     MatDestroy(&nlsP->M); 
1083:   }

1085:   PetscFree(tao->data); 
1086:   tao->data = PETSC_NULL;

1088:   return(0);
1089: }

1091: /*------------------------------------------------------------*/
1094: static PetscErrorCode TaoSetFromOptions_NLS(TaoSolver tao)
1095: {
1096:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

1100:   PetscOptionsHead("Newton line search method for unconstrained optimization"); 
1101:   PetscOptionsEList("-tao_nls_ksp_type", "ksp type", "", NLS_KSP, NLS_KSP_TYPES, NLS_KSP[nlsP->ksp_type], &nlsP->ksp_type, 0); 
1102:   PetscOptionsEList("-tao_nls_pc_type", "pc type", "", NLS_PC, NLS_PC_TYPES, NLS_PC[nlsP->pc_type], &nlsP->pc_type, 0); 
1103:   PetscOptionsEList("-tao_nls_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[nlsP->bfgs_scale_type], &nlsP->bfgs_scale_type, 0); 
1104:   PetscOptionsEList("-tao_nls_init_type", "radius initialization type", "", NLS_INIT, NLS_INIT_TYPES, NLS_INIT[nlsP->init_type], &nlsP->init_type, 0); 
1105:   PetscOptionsEList("-tao_nls_update_type", "radius update type", "", NLS_UPDATE, NLS_UPDATE_TYPES, NLS_UPDATE[nlsP->update_type], &nlsP->update_type, 0); 
1106:  PetscOptionsReal("-tao_nls_sval", "perturbation starting value", "", nlsP->sval, &nlsP->sval, 0); 
1107:   PetscOptionsReal("-tao_nls_imin", "minimum initial perturbation", "", nlsP->imin, &nlsP->imin, 0); 
1108:   PetscOptionsReal("-tao_nls_imax", "maximum initial perturbation", "", nlsP->imax, &nlsP->imax, 0); 
1109:   PetscOptionsReal("-tao_nls_imfac", "initial merit factor", "", nlsP->imfac, &nlsP->imfac, 0); 
1110:   PetscOptionsReal("-tao_nls_pmin", "minimum perturbation", "", nlsP->pmin, &nlsP->pmin, 0); 
1111:   PetscOptionsReal("-tao_nls_pmax", "maximum perturbation", "", nlsP->pmax, &nlsP->pmax, 0); 
1112:   PetscOptionsReal("-tao_nls_pgfac", "growth factor", "", nlsP->pgfac, &nlsP->pgfac, 0); 
1113:   PetscOptionsReal("-tao_nls_psfac", "shrink factor", "", nlsP->psfac, &nlsP->psfac, 0); 
1114:   PetscOptionsReal("-tao_nls_pmgfac", "merit growth factor", "", nlsP->pmgfac, &nlsP->pmgfac, 0); 
1115:   PetscOptionsReal("-tao_nls_pmsfac", "merit shrink factor", "", nlsP->pmsfac, &nlsP->pmsfac, 0); 
1116:   PetscOptionsReal("-tao_nls_eta1", "poor steplength; reduce radius", "", nlsP->eta1, &nlsP->eta1, 0); 
1117:   PetscOptionsReal("-tao_nls_eta2", "reasonable steplength; leave radius alone", "", nlsP->eta2, &nlsP->eta2, 0); 
1118:   PetscOptionsReal("-tao_nls_eta3", "good steplength; increase radius", "", nlsP->eta3, &nlsP->eta3, 0); 
1119:   PetscOptionsReal("-tao_nls_eta4", "excellent steplength; greatly increase radius", "", nlsP->eta4, &nlsP->eta4, 0); 
1120:   PetscOptionsReal("-tao_nls_alpha1", "", "", nlsP->alpha1, &nlsP->alpha1, 0); 
1121:   PetscOptionsReal("-tao_nls_alpha2", "", "", nlsP->alpha2, &nlsP->alpha2, 0); 
1122:   PetscOptionsReal("-tao_nls_alpha3", "", "", nlsP->alpha3, &nlsP->alpha3, 0); 
1123:   PetscOptionsReal("-tao_nls_alpha4", "", "", nlsP->alpha4, &nlsP->alpha4, 0); 
1124:   PetscOptionsReal("-tao_nls_alpha5", "", "", nlsP->alpha5, &nlsP->alpha5, 0); 
1125:   PetscOptionsReal("-tao_nls_nu1", "poor steplength; reduce radius", "", nlsP->nu1, &nlsP->nu1, 0); 
1126:   PetscOptionsReal("-tao_nls_nu2", "reasonable steplength; leave radius alone", "", nlsP->nu2, &nlsP->nu2, 0); 
1127:   PetscOptionsReal("-tao_nls_nu3", "good steplength; increase radius", "", nlsP->nu3, &nlsP->nu3, 0); 
1128:   PetscOptionsReal("-tao_nls_nu4", "excellent steplength; greatly increase radius", "", nlsP->nu4, &nlsP->nu4, 0); 
1129:   PetscOptionsReal("-tao_nls_omega1", "", "", nlsP->omega1, &nlsP->omega1, 0); 
1130:   PetscOptionsReal("-tao_nls_omega2", "", "", nlsP->omega2, &nlsP->omega2, 0); 
1131:   PetscOptionsReal("-tao_nls_omega3", "", "", nlsP->omega3, &nlsP->omega3, 0); 
1132:   PetscOptionsReal("-tao_nls_omega4", "", "", nlsP->omega4, &nlsP->omega4, 0); 
1133:   PetscOptionsReal("-tao_nls_omega5", "", "", nlsP->omega5, &nlsP->omega5, 0); 
1134:   PetscOptionsReal("-tao_nls_mu1_i", "", "", nlsP->mu1_i, &nlsP->mu1_i, 0); 
1135:   PetscOptionsReal("-tao_nls_mu2_i", "", "", nlsP->mu2_i, &nlsP->mu2_i, 0); 
1136:   PetscOptionsReal("-tao_nls_gamma1_i", "", "", nlsP->gamma1_i, &nlsP->gamma1_i, 0); 
1137:   PetscOptionsReal("-tao_nls_gamma2_i", "", "", nlsP->gamma2_i, &nlsP->gamma2_i, 0); 
1138:   PetscOptionsReal("-tao_nls_gamma3_i", "", "", nlsP->gamma3_i, &nlsP->gamma3_i, 0); 
1139:   PetscOptionsReal("-tao_nls_gamma4_i", "", "", nlsP->gamma4_i, &nlsP->gamma4_i, 0); 
1140:   PetscOptionsReal("-tao_nls_theta_i", "", "", nlsP->theta_i, &nlsP->theta_i, 0); 
1141:   PetscOptionsReal("-tao_nls_mu1", "", "", nlsP->mu1, &nlsP->mu1, 0); 
1142:   PetscOptionsReal("-tao_nls_mu2", "", "", nlsP->mu2, &nlsP->mu2, 0); 
1143:   PetscOptionsReal("-tao_nls_gamma1", "", "", nlsP->gamma1, &nlsP->gamma1, 0); 
1144:   PetscOptionsReal("-tao_nls_gamma2", "", "", nlsP->gamma2, &nlsP->gamma2, 0); 
1145:   PetscOptionsReal("-tao_nls_gamma3", "", "", nlsP->gamma3, &nlsP->gamma3, 0); 
1146:   PetscOptionsReal("-tao_nls_gamma4", "", "", nlsP->gamma4, &nlsP->gamma4, 0); 
1147:   PetscOptionsReal("-tao_nls_theta", "", "", nlsP->theta, &nlsP->theta, 0); 
1148:   PetscOptionsReal("-tao_nls_min_radius", "lower bound on initial radius", "", nlsP->min_radius, &nlsP->min_radius, 0); 
1149:   PetscOptionsReal("-tao_nls_max_radius", "upper bound on radius", "", nlsP->max_radius, &nlsP->max_radius, 0); 
1150:   PetscOptionsReal("-tao_nls_epsilon", "tolerance used when computing actual and predicted reduction", "", nlsP->epsilon, &nlsP->epsilon, 0); 
1151:   PetscOptionsTail(); 
1152:   TaoLineSearchSetFromOptions(tao->linesearch);
1153:   KSPSetFromOptions(tao->ksp); 
1154:   return(0);
1155: }


1158: /*------------------------------------------------------------*/
1161: static PetscErrorCode TaoView_NLS(TaoSolver tao, PetscViewer viewer)
1162: {
1163:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;
1164:   PetscInt nrejects;
1165:   PetscBool isascii;


1170:   PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1171:   if (isascii) {
1172:     PetscViewerASCIIPushTab(viewer); 
1173:     if (NLS_PC_BFGS == nlsP->pc_type && nlsP->M) {
1174:       MatLMVMGetRejects(nlsP->M,&nrejects); 
1175:       PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n",nrejects); 
1176:     }
1177:     PetscViewerASCIIPrintf(viewer, "Newton steps: %D\n", nlsP->newt); 
1178:     PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", nlsP->bfgs); 
1179:     PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", nlsP->sgrad); 
1180:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", nlsP->grad); 

1182:     PetscViewerASCIIPrintf(viewer, "nls ksp atol: %D\n", nlsP->ksp_atol); 
1183:     PetscViewerASCIIPrintf(viewer, "nls ksp rtol: %D\n", nlsP->ksp_rtol); 
1184:     PetscViewerASCIIPrintf(viewer, "nls ksp ctol: %D\n", nlsP->ksp_ctol); 
1185:     PetscViewerASCIIPrintf(viewer, "nls ksp negc: %D\n", nlsP->ksp_negc); 
1186:     PetscViewerASCIIPrintf(viewer, "nls ksp dtol: %D\n", nlsP->ksp_dtol); 
1187:     PetscViewerASCIIPrintf(viewer, "nls ksp iter: %D\n", nlsP->ksp_iter); 
1188:     PetscViewerASCIIPrintf(viewer, "nls ksp othr: %D\n", nlsP->ksp_othr); 
1189:     PetscViewerASCIIPopTab(viewer); 
1190:   } else {
1191:     SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO NLS",((PetscObject)viewer)->type_name);
1192:   }
1193:   return(0);
1194: }

1196: /* ---------------------------------------------------------- */
1200: PetscErrorCode TaoCreate_NLS(TaoSolver tao)
1201: {
1202:   TAO_NLS *nlsP;
1203:   const char *morethuente_type = TAOLINESEARCH_MT;

1207:   PetscNewLog(tao,TAO_NLS,&nlsP); 

1209:   tao->ops->setup = TaoSetUp_NLS;
1210:   tao->ops->solve = TaoSolve_NLS;
1211:   tao->ops->view = TaoView_NLS;
1212:   tao->ops->setfromoptions = TaoSetFromOptions_NLS;
1213:   tao->ops->destroy = TaoDestroy_NLS;

1215:   tao->max_it = 50;
1216:   tao->fatol = 1e-10;
1217:   tao->frtol = 1e-10;
1218:   tao->data = (void*)nlsP;
1219:   tao->trust0 = 100.0;

1221:   nlsP->sval   = 0.0;
1222:   nlsP->imin   = 1.0e-4;
1223:   nlsP->imax   = 1.0e+2;
1224:   nlsP->imfac  = 1.0e-1;

1226:   nlsP->pmin   = 1.0e-12;
1227:   nlsP->pmax   = 1.0e+2;
1228:   nlsP->pgfac  = 1.0e+1;
1229:   nlsP->psfac  = 4.0e-1;
1230:   nlsP->pmgfac = 1.0e-1;
1231:   nlsP->pmsfac = 1.0e-1;

1233:   /*  Default values for trust-region radius update based on steplength */
1234:   nlsP->nu1 = 0.25;
1235:   nlsP->nu2 = 0.50;
1236:   nlsP->nu3 = 1.00;
1237:   nlsP->nu4 = 1.25;

1239:   nlsP->omega1 = 0.25;
1240:   nlsP->omega2 = 0.50;
1241:   nlsP->omega3 = 1.00;
1242:   nlsP->omega4 = 2.00;
1243:   nlsP->omega5 = 4.00;

1245:   /*  Default values for trust-region radius update based on reduction */
1246:   nlsP->eta1 = 1.0e-4;
1247:   nlsP->eta2 = 0.25;
1248:   nlsP->eta3 = 0.50;
1249:   nlsP->eta4 = 0.90;

1251:   nlsP->alpha1 = 0.25;
1252:   nlsP->alpha2 = 0.50;
1253:   nlsP->alpha3 = 1.00;
1254:   nlsP->alpha4 = 2.00;
1255:   nlsP->alpha5 = 4.00;

1257:   /*  Default values for trust-region radius update based on interpolation */
1258:   nlsP->mu1 = 0.10;
1259:   nlsP->mu2 = 0.50;

1261:   nlsP->gamma1 = 0.25;
1262:   nlsP->gamma2 = 0.50;
1263:   nlsP->gamma3 = 2.00;
1264:   nlsP->gamma4 = 4.00;

1266:   nlsP->theta = 0.05;

1268:   /*  Default values for trust region initialization based on interpolation */
1269:   nlsP->mu1_i = 0.35;
1270:   nlsP->mu2_i = 0.50;

1272:   nlsP->gamma1_i = 0.0625;
1273:   nlsP->gamma2_i = 0.5;
1274:   nlsP->gamma3_i = 2.0;
1275:   nlsP->gamma4_i = 5.0;
1276:   
1277:   nlsP->theta_i = 0.25;

1279:   /*  Remaining parameters */
1280:   nlsP->min_radius = 1.0e-10;
1281:   nlsP->max_radius = 1.0e10;
1282:   nlsP->epsilon = 1.0e-6;

1284:   nlsP->ksp_type        = NLS_KSP_STCG;
1285:   nlsP->pc_type         = NLS_PC_BFGS;
1286:   nlsP->bfgs_scale_type = BFGS_SCALE_PHESS;
1287:   nlsP->init_type       = NLS_INIT_INTERPOLATION;
1288:   nlsP->update_type     = NLS_UPDATE_STEP;

1290:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch); 
1291:   TaoLineSearchSetType(tao->linesearch,morethuente_type); 
1292:   TaoLineSearchUseTaoSolverRoutines(tao->linesearch,tao); 

1294:   /*  Set linear solver to default for symmetric matrices */
1295:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp); 

1297:   return(0);
1298: }





1307: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x) 
1308: {
1310:     Mat M;
1315:     PCShellGetContext(pc,(void**)&M); 
1316:     MatLMVMSolve(M, b, x); 
1317:     return(0);
1318: }