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 != TAOLINESEARCH_SUCCESS &&
692:            ls_reason != TAOLINESEARCH_SUCCESS_USER &&
693:            stepType != NLS_GRADIENT) {      /* Linesearch failed */
694:       f = fold;
695:       VecCopy(nlsP->Xold, tao->solution); 
696:       VecCopy(nlsP->Gold, tao->gradient); 

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

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

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

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

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

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

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

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

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

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

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

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

929:       default:
930:         if (stepType == NLS_NEWTON) {

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

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

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

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

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

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

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

1042: /* ---------------------------------------------------------- */
1045: static PetscErrorCode TaoSetUp_NLS(TaoSolver tao)
1046: {
1047:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;


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

1059:   nlsP->Diag = 0;
1060:   nlsP->M = 0;

1062:   return(0);
1063: }

1065: /*------------------------------------------------------------*/
1068: static PetscErrorCode TaoDestroy_NLS(TaoSolver tao)
1069: {
1070:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

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

1087:   PetscFree(tao->data); 
1088:   tao->data = PETSC_NULL;

1090:   return(0);
1091: }

1093: /*------------------------------------------------------------*/
1096: static PetscErrorCode TaoSetFromOptions_NLS(TaoSolver tao)
1097: {
1098:   TAO_NLS *nlsP = (TAO_NLS *)tao->data;

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


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


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

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

1198: /* ---------------------------------------------------------- */
1202: PetscErrorCode TaoCreate_NLS(TaoSolver tao)
1203: {
1204:   TAO_NLS *nlsP;
1205:   const char *morethuente_type = TAOLINESEARCH_MT;

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

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

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

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

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

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

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

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

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

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

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

1268:   nlsP->theta = 0.05;

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

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

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

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

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

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

1299:   return(0);
1300: }





1309: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x) 
1310: {
1312:     Mat M;
1317:     PCShellGetContext(pc,(void**)&M); 
1318:     MatLMVMSolve(M, b, x); 
1319:     return(0);
1320: }