Actual source code: ls.c

  1: /*$Id: ls.c,v 1.172 2001/08/07 03:04:11 balay Exp $*/

 3:  #include src/snes/impls/ls/ls.h

  5: /*
  6:      Checks if J^T F = 0 which implies we've found a local minimum of the function,
  7:     but not a zero. In the case when one cannot compute J^T F we use the fact that
  8:     0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More 
  9:     for this trick.
 10: */
 11: #undef __FUNCT__  
 13: int SNESLSCheckLocalMin_Private(Mat A,Vec F,Vec W,PetscReal fnorm,PetscTruth *ismin)
 14: {
 15:   PetscReal  a1;
 16:   int        ierr;
 17:   PetscTruth hastranspose;

 20:   *ismin = PETSC_FALSE;
 21:   MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
 22:   if (hastranspose) {
 23:     /* Compute || J^T F|| */
 24:     MatMultTranspose(A,F,W);
 25:     VecNorm(W,NORM_2,&a1);
 26:     PetscLogInfo(0,"SNESSolve_LS: || J^T F|| %g near zero implies found a local minimumn",a1/fnorm);
 27:     if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
 28:   } else {
 29:     Vec       work;
 30:     PetscScalar    result;
 31:     PetscReal wnorm;

 33:     VecSetRandom(PETSC_NULL,W);
 34:     VecNorm(W,NORM_2,&wnorm);
 35:     VecDuplicate(W,&work);
 36:     MatMult(A,W,work);
 37:     VecDot(F,work,&result);
 38:     VecDestroy(work);
 39:     a1   = PetscAbsScalar(result)/(fnorm*wnorm);
 40:     PetscLogInfo(0,"SNESSolve_LS: (F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimumn",a1);
 41:     if (a1 < 1.e-4) *ismin = PETSC_TRUE;
 42:   }
 43:   return(0);
 44: }

 46: /*
 47:      Checks if J^T(F - J*X) = 0 
 48: */
 49: #undef __FUNCT__  
 51: int SNESLSCheckResidual_Private(Mat A,Vec F,Vec X,Vec W1,Vec W2)
 52: {
 53:   PetscReal     a1,a2;
 54:   int           ierr;
 55:   PetscTruth    hastranspose;
 56:   PetscScalar   mone = -1.0;

 59:   MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
 60:   if (hastranspose) {
 61:     MatMult(A,X,W1);
 62:     VecAXPY(&mone,F,W1);

 64:     /* Compute || J^T W|| */
 65:     MatMultTranspose(A,W1,W2);
 66:     VecNorm(W1,NORM_2,&a1);
 67:     VecNorm(W2,NORM_2,&a2);
 68:     if (a1 != 0) {
 69:       PetscLogInfo(0,"SNESSolve_LS: ||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhsn",a2/a1);
 70:     }
 71:   }
 72:   return(0);
 73: }

 75: /*  -------------------------------------------------------------------- 

 77:      This file implements a truncated Newton method with a line search,
 78:      for solving a system of nonlinear equations, using the SLES, Vec, 
 79:      and Mat interfaces for linear solvers, vectors, and matrices, 
 80:      respectively.

 82:      The following basic routines are required for each nonlinear solver:
 83:           SNESCreate_XXX()          - Creates a nonlinear solver context
 84:           SNESSetFromOptions_XXX()  - Sets runtime options
 85:           SNESSolve_XXX()           - Solves the nonlinear system
 86:           SNESDestroy_XXX()         - Destroys the nonlinear solver context
 87:      The suffix "_XXX" denotes a particular implementation, in this case
 88:      we use _LS (e.g., SNESCreate_LS, SNESSolve_LS) for solving
 89:      systems of nonlinear equations with a line search (LS) method.
 90:      These routines are actually called via the common user interface
 91:      routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and 
 92:      SNESDestroy(), so the application code interface remains identical 
 93:      for all nonlinear solvers.

 95:      Another key routine is:
 96:           SNESSetUp_XXX()           - Prepares for the use of a nonlinear solver
 97:      by setting data structures and options.   The interface routine SNESSetUp()
 98:      is not usually called directly by the user, but instead is called by
 99:      SNESSolve() if necessary.

101:      Additional basic routines are:
102:           SNESView_XXX()            - Prints details of runtime options that
103:                                       have actually been used.
104:      These are called by application codes via the interface routines
105:      SNESView().

107:      The various types of solvers (preconditioners, Krylov subspace methods,
108:      nonlinear solvers, timesteppers) are all organized similarly, so the
109:      above description applies to these categories also.  

111:     -------------------------------------------------------------------- */
112: /*
113:    SNESSolve_LS - Solves a nonlinear system with a truncated Newton
114:    method with a line search.

116:    Input Parameters:
117: .  snes - the SNES context

119:    Output Parameter:
120: .  outits - number of iterations until termination

122:    Application Interface Routine: SNESSolve()

124:    Notes:
125:    This implements essentially a truncated Newton method with a
126:    line search.  By default a cubic backtracking line search 
127:    is employed, as described in the text "Numerical Methods for
128:    Unconstrained Optimization and Nonlinear Equations" by Dennis 
129:    and Schnabel.
130: */
131: #undef __FUNCT__  
133: int SNESSolve_LS(SNES snes,int *outits)
134: {
135:   SNES_LS      *neP = (SNES_LS*)snes->data;
136:   int          maxits,i,ierr,lits,lsfail;
137:   MatStructure flg = DIFFERENT_NONZERO_PATTERN;
138:   PetscReal    fnorm,gnorm,xnorm,ynorm;
139:   Vec          Y,X,F,G,W,TMP;

142:   snes->reason  = SNES_CONVERGED_ITERATING;

144:   maxits        = snes->max_its;        /* maximum number of iterations */
145:   X                = snes->vec_sol;        /* solution vector */
146:   F                = snes->vec_func;        /* residual vector */
147:   Y                = snes->work[0];        /* work vectors */
148:   G                = snes->work[1];
149:   W                = snes->work[2];

151:   PetscObjectTakeAccess(snes);
152:   snes->iter = 0;
153:   PetscObjectGrantAccess(snes);
154:   SNESComputeFunction(snes,X,F);  /*  F(X)      */
155:   VecNorm(F,NORM_2,&fnorm);        /* fnorm <- ||F||  */
156:   PetscObjectTakeAccess(snes);
157:   snes->norm = fnorm;
158:   PetscObjectGrantAccess(snes);
159:   SNESLogConvHistory(snes,fnorm,0);
160:   SNESMonitor(snes,0,fnorm);

162:   if (fnorm < snes->atol) {*outits = 0; snes->reason = SNES_CONVERGED_FNORM_ABS; return(0);}

164:   /* set parameter for default relative tolerance convergence test */
165:   snes->ttol = fnorm*snes->rtol;

167:   for (i=0; i<maxits; i++) {

169:     /* Call general purpose update function */
170:     if (snes->update != PETSC_NULL) {
171:       (*snes->update)(snes, snes->iter);
172:     }

174:     /* Solve J Y = F, where J is Jacobian matrix */
175:     SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
176:     SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);
177:     SLESSolve(snes->sles,F,Y,&lits);

179:     if (PetscLogPrintInfo){
180:       SNESLSCheckResidual_Private(snes->jacobian,F,Y,G,W);
181:     }

183:     /* should check what happened to the linear solve? */
184:     snes->linear_its += lits;
185:     PetscLogInfo(snes,"SNESSolve_LS: iter=%d, linear solve iterations=%dn",snes->iter,lits);

187:     /* Compute a (scaled) negative update in the line search routine: 
188:          Y <- X - lambda*Y 
189:        and evaluate G(Y) = function(Y)) 
190:     */
191:     VecCopy(Y,snes->vec_sol_update_always);
192:     (*neP->LineSearch)(snes,neP->lsP,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lsfail);
193:     PetscLogInfo(snes,"SNESSolve_LS: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lsfail=%dn",fnorm,gnorm,ynorm,lsfail);

195:     TMP = F; F = G; snes->vec_func_always = F; G = TMP;
196:     TMP = X; X = Y; snes->vec_sol_always = X;  Y = TMP;
197:     fnorm = gnorm;

199:     PetscObjectTakeAccess(snes);
200:     snes->iter = i+1;
201:     snes->norm = fnorm;
202:     PetscObjectGrantAccess(snes);
203:     SNESLogConvHistory(snes,fnorm,lits);
204:     SNESMonitor(snes,i+1,fnorm);

206:     if (lsfail) {
207:       PetscTruth ismin;

209:       if (++snes->numFailures >= snes->maxFailures) {
210:         snes->reason = SNES_DIVERGED_LS_FAILURE;
211:         SNESLSCheckLocalMin_Private(snes->jacobian,F,W,fnorm,&ismin);
212:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
213:         break;
214:       }
215:     }

217:     /* Test for convergence */
218:     if (snes->converged) {
219:       VecNorm(X,NORM_2,&xnorm);        /* xnorm = || X || */
220:       (*snes->converged)(snes,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
221:       if (snes->reason) {
222:         break;
223:       }
224:     }
225:   }
226:   if (X != snes->vec_sol) {
227:     VecCopy(X,snes->vec_sol);
228:   }
229:   if (F != snes->vec_func) {
230:     VecCopy(F,snes->vec_func);
231:   }
232:   snes->vec_sol_always  = snes->vec_sol;
233:   snes->vec_func_always = snes->vec_func;
234:   if (i == maxits) {
235:     PetscLogInfo(snes,"SNESSolve_LS: Maximum number of iterations has been reached: %dn",maxits);
236:     i--;
237:     snes->reason = SNES_DIVERGED_MAX_IT;
238:   }
239:   PetscObjectTakeAccess(snes);
240:   PetscObjectGrantAccess(snes);
241:   *outits = i+1;
242:   return(0);
243: }
244: /* -------------------------------------------------------------------------- */
245: /*
246:    SNESSetUp_LS - Sets up the internal data structures for the later use
247:    of the SNESLS nonlinear solver.

249:    Input Parameter:
250: .  snes - the SNES context
251: .  x - the solution vector

253:    Application Interface Routine: SNESSetUp()

255:    Notes:
256:    For basic use of the SNES solvers, the user need not explicitly call
257:    SNESSetUp(), since these actions will automatically occur during
258:    the call to SNESSolve().
259:  */
260: #undef __FUNCT__  
262: int SNESSetUp_LS(SNES snes)
263: {

267:   snes->nwork = 4;
268:   VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);
269:   PetscLogObjectParents(snes,snes->nwork,snes->work);
270:   snes->vec_sol_update_always = snes->work[3];
271:   return(0);
272: }
273: /* -------------------------------------------------------------------------- */
274: /*
275:    SNESDestroy_LS - Destroys the private SNES_LS context that was created
276:    with SNESCreate_LS().

278:    Input Parameter:
279: .  snes - the SNES context

281:    Application Interface Routine: SNESDestroy()
282:  */
283: #undef __FUNCT__  
285: int SNESDestroy_LS(SNES snes)
286: {
287:   int  ierr;

290:   if (snes->nwork) {
291:     VecDestroyVecs(snes->work,snes->nwork);
292:   }
293:   PetscFree(snes->data);
294:   return(0);
295: }
296: /* -------------------------------------------------------------------------- */
297: #undef __FUNCT__  

300: /*@C
301:    SNESNoLineSearch - This routine is not a line search at all; 
302:    it simply uses the full Newton step.  Thus, this routine is intended 
303:    to serve as a template and is not recommended for general use.  

305:    Collective on SNES and Vec

307:    Input Parameters:
308: +  snes - nonlinear context
309: .  lsctx - optional context for line search (not used here)
310: .  x - current iterate
311: .  f - residual evaluated at x
312: .  y - search direction (contains new iterate on output)
313: .  w - work vector
314: -  fnorm - 2-norm of f

316:    Output Parameters:
317: +  g - residual evaluated at new iterate y
318: .  y - new iterate (contains search direction on input)
319: .  gnorm - 2-norm of g
320: .  ynorm - 2-norm of search length
321: -  flag - set to 0, indicating a successful line search

323:    Options Database Key:
324: .  -snes_ls basic - Activates SNESNoLineSearch()

326:    Level: advanced

328: .keywords: SNES, nonlinear, line search, cubic

330: .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 
331:           SNESSetLineSearch(), SNESNoLineSearchNoNorms()
332: @*/
333: int SNESNoLineSearch(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,int *flag)
334: {
335:   int         ierr;
336:   PetscScalar mone = -1.0;
337:   SNES_LS     *neP = (SNES_LS*)snes->data;
338:   PetscTruth  change_y = PETSC_FALSE;

341:   *flag = 0;
342:   PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
343:   VecNorm(y,NORM_2,ynorm);  /* ynorm = || y || */
344:   VecAYPX(&mone,x,y);            /* y <- y - x      */
345:   if (neP->CheckStep) {
346:    (*neP->CheckStep)(snes,neP->checkP,y,&change_y);
347:   }
348:   SNESComputeFunction(snes,y,g); /* Compute F(y)    */
349:   VecNorm(g,NORM_2,gnorm);  /* gnorm = || g || */
350:   PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
351:   return(0);
352: }
353: /* -------------------------------------------------------------------------- */
354: #undef __FUNCT__  

357: /*@C
358:    SNESNoLineSearchNoNorms - This routine is not a line search at 
359:    all; it simply uses the full Newton step. This version does not
360:    even compute the norm of the function or search direction; this
361:    is intended only when you know the full step is fine and are
362:    not checking for convergence of the nonlinear iteration (for
363:    example, you are running always for a fixed number of Newton steps).

365:    Collective on SNES and Vec

367:    Input Parameters:
368: +  snes - nonlinear context
369: .  lsctx - optional context for line search (not used here)
370: .  x - current iterate
371: .  f - residual evaluated at x
372: .  y - search direction (contains new iterate on output)
373: .  w - work vector
374: -  fnorm - 2-norm of f

376:    Output Parameters:
377: +  g - residual evaluated at new iterate y
378: .  gnorm - not changed
379: .  ynorm - not changed
380: -  flag - set to 0, indicating a successful line search

382:    Options Database Key:
383: .  -snes_ls basicnonorms - Activates SNESNoLineSearchNoNorms()

385:    Notes:
386:    SNESNoLineSearchNoNorms() must be used in conjunction with
387:    either the options
388: $     -snes_no_convergence_test -snes_max_it <its>
389:    or alternatively a user-defined custom test set via
390:    SNESSetConvergenceTest(); or a -snes_max_it of 1, 
391:    otherwise, the SNES solver will generate an error.

393:    During the final iteration this will not evaluate the function at
394:    the solution point. This is to save a function evaluation while
395:    using pseudo-timestepping.

397:    The residual norms printed by monitoring routines such as
398:    SNESDefaultMonitor() (as activated via -snes_monitor) will not be 
399:    correct, since they are not computed.

401:    Level: advanced

403: .keywords: SNES, nonlinear, line search, cubic

405: .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 
406:           SNESSetLineSearch(), SNESNoLineSearch()
407: @*/
408: int SNESNoLineSearchNoNorms(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,int *flag)
409: {
410:   int         ierr;
411:   PetscScalar mone = -1.0;
412:   SNES_LS     *neP = (SNES_LS*)snes->data;
413:   PetscTruth  change_y = PETSC_FALSE;

416:   *flag = 0;
417:   PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
418:   VecAYPX(&mone,x,y);            /* y <- y - x      */
419:   if (neP->CheckStep) {
420:    (*neP->CheckStep)(snes,neP->checkP,y,&change_y);
421:   }
422: 
423:   /* don't evaluate function the last time through */
424:   if (snes->iter < snes->max_its-1) {
425:     SNESComputeFunction(snes,y,g); /* Compute F(y)    */
426:   }
427:   PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
428:   return(0);
429: }
430: /* -------------------------------------------------------------------------- */
431: #undef __FUNCT__  
433: /*@C
434:    SNESCubicLineSearch - Performs a cubic line search (default line search method).

436:    Collective on SNES

438:    Input Parameters:
439: +  snes - nonlinear context
440: .  lsctx - optional context for line search (not used here)
441: .  x - current iterate
442: .  f - residual evaluated at x
443: .  y - search direction (contains new iterate on output)
444: .  w - work vector
445: -  fnorm - 2-norm of f

447:    Output Parameters:
448: +  g - residual evaluated at new iterate y
449: .  y - new iterate (contains search direction on input)
450: .  gnorm - 2-norm of g
451: .  ynorm - 2-norm of search length
452: -  flag - 0 if line search succeeds; -1 on failure.

454:    Options Database Key:
455: $  -snes_ls cubic - Activates SNESCubicLineSearch()

457:    Notes:
458:    This line search is taken from "Numerical Methods for Unconstrained 
459:    Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.

461:    Level: advanced

463: .keywords: SNES, nonlinear, line search, cubic

465: .seealso: SNESQuadraticLineSearch(), SNESNoLineSearch(), SNESSetLineSearch(), SNESNoLineSearchNoNorms()
466: @*/
467: int SNESCubicLineSearch(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,int *flag)
468: {
469:   /* 
470:      Note that for line search purposes we work with with the related
471:      minimization problem:
472:         min  z(x):  R^n -> R,
473:      where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2.
474:    */
475: 
476:   PetscReal   steptol,initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength;
477:   PetscReal   maxstep,minlambda,alpha,lambda,lambdatemp,lambdaneg;
478: #if defined(PETSC_USE_COMPLEX)
479:   PetscScalar cinitslope,clambda;
480: #endif
481:   int         ierr,count;
482:   SNES_LS     *neP = (SNES_LS*)snes->data;
483:   PetscScalar mone = -1.0,scale;
484:   PetscTruth  change_y = PETSC_FALSE;

487:   PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
488:   *flag   = 0;
489:   alpha   = neP->alpha;
490:   maxstep = neP->maxstep;
491:   steptol = neP->steptol;

493:   VecNorm(y,NORM_2,ynorm);
494:   if (*ynorm < snes->atol) {
495:     PetscLogInfo(snes,"SNESCubicLineSearch: Search direction and size are nearly 0n");
496:     *gnorm = fnorm;
497:     VecCopy(x,y);
498:     VecCopy(f,g);
499:     goto theend1;
500:   }
501:   if (*ynorm > maxstep) {        /* Step too big, so scale back */
502:     scale = maxstep/(*ynorm);
503: #if defined(PETSC_USE_COMPLEX)
504:     PetscLogInfo(snes,"SNESCubicLineSearch: Scaling step by %g old ynorm %gn",PetscRealPart(scale),*ynorm);
505: #else
506:     PetscLogInfo(snes,"SNESCubicLineSearch: Scaling step by %g old ynorm %gn",scale,*ynorm);
507: #endif
508:     VecScale(&scale,y);
509:     *ynorm = maxstep;
510:   }
511:   ierr      = VecMaxPointwiseDivide(y,x,&rellength);
512:   minlambda = steptol/rellength;
513:   MatMult(snes->jacobian,y,w);
514: #if defined(PETSC_USE_COMPLEX)
515:   VecDot(f,w,&cinitslope);
516:   initslope = PetscRealPart(cinitslope);
517: #else
518:   VecDot(f,w,&initslope);
519: #endif
520:   if (initslope > 0.0) initslope = -initslope;
521:   if (initslope == 0.0) initslope = -1.0;

523:   VecCopy(y,w);
524:   VecAYPX(&mone,x,w);
525:   SNESComputeFunction(snes,w,g);
526:   VecNorm(g,NORM_2,gnorm);
527:   if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */
528:     VecCopy(w,y);
529:     PetscLogInfo(snes,"SNESCubicLineSearch: Using full stepn");
530:     goto theend1;
531:   }

533:   /* Fit points with quadratic */
534:   lambda = 1.0;
535:   lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
536:   lambdaprev = lambda;
537:   gnormprev = *gnorm;
538:   if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
539:   if (lambdatemp <= .1*lambda) lambda = .1*lambda;
540:   else                         lambda = lambdatemp;
541:   ierr   = VecCopy(x,w);
542:   lambdaneg = -lambda;
543: #if defined(PETSC_USE_COMPLEX)
544:   clambda = lambdaneg; VecAXPY(&clambda,y,w);
545: #else
546:   VecAXPY(&lambdaneg,y,w);
547: #endif
548:   SNESComputeFunction(snes,w,g);
549:   VecNorm(g,NORM_2,gnorm);
550:   if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
551:     VecCopy(w,y);
552:     PetscLogInfo(snes,"SNESCubicLineSearch: Quadratically determined step, lambda=%18.16en",lambda);
553:     goto theend1;
554:   }

556:   /* Fit points with cubic */
557:   count = 1;
558:   while (PETSC_TRUE) {
559:     if (lambda <= minlambda) { /* bad luck; use full step */
560:       PetscLogInfo(snes,"SNESCubicLineSearch:Unable to find good step length! %d n",count);
561:       PetscLogInfo(snes,"SNESCubicLineSearch:fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16en",fnorm,*gnorm,*ynorm,lambda,initslope);
562:       VecCopy(x,y);
563:       *flag = -1; break;
564:     }
565:     t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope;
566:     t2 = .5*(gnormprev*gnormprev  - fnorm*fnorm) - lambdaprev*initslope;
567:     a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
568:     b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
569:     d  = b*b - 3*a*initslope;
570:     if (d < 0.0) d = 0.0;
571:     if (a == 0.0) {
572:       lambdatemp = -initslope/(2.0*b);
573:     } else {
574:       lambdatemp = (-b + sqrt(d))/(3.0*a);
575:     }
576:     lambdaprev = lambda;
577:     gnormprev  = *gnorm;
578:     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
579:     if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
580:     else                         lambda     = lambdatemp;
581:     VecCopy(x,w);
582:     lambdaneg = -lambda;
583: #if defined(PETSC_USE_COMPLEX)
584:     clambda = lambdaneg;
585:     VecAXPY(&clambda,y,w);
586: #else
587:     VecAXPY(&lambdaneg,y,w);
588: #endif
589:     SNESComputeFunction(snes,w,g);
590:     VecNorm(g,NORM_2,gnorm);
591:     if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */
592:       VecCopy(w,y);
593:       PetscLogInfo(snes,"SNESCubicLineSearch: Cubically determined step, lambda=%18.16en",lambda);
594:       goto theend1;
595:     } else {
596:       PetscLogInfo(snes,"SNESCubicLineSearch: Cubic step no good, shrinking lambda,  lambda=%18.16en",lambda);
597:     }
598:     count++;
599:   }
600:   theend1:
601:   /* Optional user-defined check for line search step validity */
602:   if (neP->CheckStep) {
603:     (*neP->CheckStep)(snes,neP->checkP,y,&change_y);
604:     if (change_y == PETSC_TRUE) { /* recompute the function if the step has changed */
605:       SNESComputeFunction(snes,y,g);
606:       VecNormBegin(y,NORM_2,ynorm);
607:       VecNormBegin(g,NORM_2,gnorm);
608:       VecNormEnd(y,NORM_2,ynorm);
609:       VecNormEnd(g,NORM_2,gnorm);
610:     }
611:   }
612:   PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
613:   return(0);
614: }
615: /* -------------------------------------------------------------------------- */
616: #undef __FUNCT__  
618: /*@C
619:    SNESQuadraticLineSearch - Performs a quadratic line search.

621:    Collective on SNES and Vec

623:    Input Parameters:
624: +  snes - the SNES context
625: .  lsctx - optional context for line search (not used here)
626: .  x - current iterate
627: .  f - residual evaluated at x
628: .  y - search direction (contains new iterate on output)
629: .  w - work vector
630: -  fnorm - 2-norm of f

632:    Output Parameters:
633: +  g - residual evaluated at new iterate y
634: .  y - new iterate (contains search direction on input)
635: .  gnorm - 2-norm of g
636: .  ynorm - 2-norm of search length
637: -  flag - 0 if line search succeeds; -1 on failure.

639:    Options Database Key:
640: .  -snes_ls quadratic - Activates SNESQuadraticLineSearch()

642:    Notes:
643:    Use SNESSetLineSearch() to set this routine within the SNESLS method.  

645:    Level: advanced

647: .keywords: SNES, nonlinear, quadratic, line search

649: .seealso: SNESCubicLineSearch(), SNESNoLineSearch(), SNESSetLineSearch(), SNESNoLineSearchNoNorms()
650: @*/
651: int SNESQuadraticLineSearch(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,int *flag)
652: {
653:   /* 
654:      Note that for line search purposes we work with with the related
655:      minimization problem:
656:         min  z(x):  R^n -> R,
657:      where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2.
658:    */
659:   PetscReal   steptol,initslope,maxstep,minlambda,alpha,lambda,lambdatemp,lambdaneg,rellength;
660: #if defined(PETSC_USE_COMPLEX)
661:   PetscScalar cinitslope,clambda;
662: #endif
663:   int         ierr,count;
664:   SNES_LS     *neP = (SNES_LS*)snes->data;
665:   PetscScalar mone = -1.0,scale;
666:   PetscTruth  change_y = PETSC_FALSE;

669:   PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
670:   *flag   = 0;
671:   alpha   = neP->alpha;
672:   maxstep = neP->maxstep;
673:   steptol = neP->steptol;

675:   VecNorm(y,NORM_2,ynorm);
676:   if (*ynorm < snes->atol) {
677:     PetscLogInfo(snes,"SNESQuadraticLineSearch: Search direction and size is 0n");
678:     *gnorm = fnorm;
679:     VecCopy(x,y);
680:     VecCopy(f,g);
681:     goto theend2;
682:   }
683:   if (*ynorm > maxstep) {        /* Step too big, so scale back */
684:     scale = maxstep/(*ynorm);
685:     VecScale(&scale,y);
686:     *ynorm = maxstep;
687:   }
688:   ierr      = VecMaxPointwiseDivide(y,x,&rellength);
689:   minlambda = steptol/rellength;
690:   MatMult(snes->jacobian,y,w);
691: #if defined(PETSC_USE_COMPLEX)
692:   VecDot(f,w,&cinitslope);
693:   initslope = PetscRealPart(cinitslope);
694: #else
695:   VecDot(f,w,&initslope);
696: #endif
697:   if (initslope > 0.0) initslope = -initslope;
698:   if (initslope == 0.0) initslope = -1.0;

700:   VecCopy(y,w);
701:   VecAYPX(&mone,x,w);
702:   SNESComputeFunction(snes,w,g);
703:   VecNorm(g,NORM_2,gnorm);
704:   if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */
705:     VecCopy(w,y);
706:     PetscLogInfo(snes,"SNESQuadraticLineSearch: Using full stepn");
707:     goto theend2;
708:   }

710:   /* Fit points with quadratic */
711:   lambda = 1.0;
712:   count = 1;
713:   while (PETSC_TRUE) {
714:     if (lambda <= minlambda) { /* bad luck; use full step */
715:       PetscLogInfo(snes,"SNESQuadraticLineSearch:Unable to find good step length! %d n",count);
716:       PetscLogInfo(snes,"SNESQuadraticLineSearch:fnorm=%g, gnorm=%g, ynorm=%g, lambda=%g, initial slope=%gn",fnorm,*gnorm,*ynorm,lambda,initslope);
717:       VecCopy(x,y);
718:       *flag = -1; break;
719:     }
720:     lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
721:     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
722:     if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
723:     else                         lambda     = lambdatemp;
724:     VecCopy(x,w);
725:     lambdaneg = -lambda;
726: #if defined(PETSC_USE_COMPLEX)
727:     clambda = lambdaneg; VecAXPY(&clambda,y,w);
728: #else
729:     VecAXPY(&lambdaneg,y,w);
730: #endif
731:     SNESComputeFunction(snes,w,g);
732:     VecNorm(g,NORM_2,gnorm);
733:     if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
734:       VecCopy(w,y);
735:       PetscLogInfo(snes,"SNESQuadraticLineSearch:Quadratically determined step, lambda=%gn",lambda);
736:       break;
737:     }
738:     count++;
739:   }
740:   theend2:
741:   /* Optional user-defined check for line search step validity */
742:   if (neP->CheckStep) {
743:     (*neP->CheckStep)(snes,neP->checkP,y,&change_y);
744:     if (change_y == PETSC_TRUE) { /* recompute the function if the step has changed */
745:       SNESComputeFunction(snes,y,g);
746:       VecNormBegin(y,NORM_2,ynorm);
747:       VecNormBegin(g,NORM_2,gnorm);
748:       VecNormEnd(y,NORM_2,ynorm);
749:       VecNormEnd(g,NORM_2,gnorm);
750:     }
751:   }
752:   PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
753:   return(0);
754: }
755: /* -------------------------------------------------------------------------- */
756: #undef __FUNCT__  
758: /*@C
759:    SNESSetLineSearch - Sets the line search routine to be used
760:    by the method SNESLS.

762:    Input Parameters:
763: +  snes - nonlinear context obtained from SNESCreate()
764: .  lsctx - optional user-defined context for use by line search 
765: -  func - pointer to int function

767:    Collective on SNES

769:    Available Routines:
770: +  SNESCubicLineSearch() - default line search
771: .  SNESQuadraticLineSearch() - quadratic line search
772: .  SNESNoLineSearch() - the full Newton step (actually not a line search)
773: -  SNESNoLineSearchNoNorms() - the full Newton step (calculating no norms; faster in parallel)

775:     Options Database Keys:
776: +   -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search
777: .   -snes_ls_alpha <alpha> - Sets alpha
778: .   -snes_ls_maxstep <max> - Sets maxstep
779: -   -snes_ls_steptol <steptol> - Sets steptol, this is the minimum step size that the line search code
780:                    will accept; min p[i]/x[i] < steptol. The -snes_stol <stol> is the minimum step length
781:                    the default convergence test will use and is based on 2-norm(p) < stol*2-norm(x)

783:    Calling sequence of func:
784: .vb
785:    func (SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,
786:          PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,*flag)
787: .ve

789:     Input parameters for func:
790: +   snes - nonlinear context
791: .   lsctx - optional user-defined context for line search
792: .   x - current iterate
793: .   f - residual evaluated at x
794: .   y - search direction (contains new iterate on output)
795: .   w - work vector
796: -   fnorm - 2-norm of f

798:     Output parameters for func:
799: +   g - residual evaluated at new iterate y
800: .   y - new iterate (contains search direction on input)
801: .   gnorm - 2-norm of g
802: .   ynorm - 2-norm of search length
803: -   flag - set to 0 if the line search succeeds; a nonzero integer 
804:            on failure.

806:     Level: advanced

808: .keywords: SNES, nonlinear, set, line search, routine

810: .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), SNESNoLineSearch(), SNESNoLineSearchNoNorms(), 
811:           SNESSetLineSearchCheck(), SNESSetLineSearchParams(), SNESGetLineSearchParams()
812: @*/
813: int SNESSetLineSearch(SNES snes,int (*func)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,int*),void *lsctx)
814: {
815:   int ierr,(*f)(SNES,int (*)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,int*),void*);

818:   PetscObjectQueryFunction((PetscObject)snes,"SNESSetLineSearch_C",(void (**)(void))&f);
819:   if (f) {
820:     (*f)(snes,func,lsctx);
821:   }
822:   return(0);
823: }
824: /* -------------------------------------------------------------------------- */
825: EXTERN_C_BEGIN
826: #undef __FUNCT__  
828: int SNESSetLineSearch_LS(SNES snes,int (*func)(SNES,void*,Vec,Vec,Vec,Vec,Vec,
829:                          PetscReal,PetscReal*,PetscReal*,int*),void *lsctx)
830: {
832:   ((SNES_LS *)(snes->data))->LineSearch = func;
833:   ((SNES_LS *)(snes->data))->lsP        = lsctx;
834:   return(0);
835: }
836: EXTERN_C_END
837: /* -------------------------------------------------------------------------- */
838: #undef __FUNCT__  
840: /*@C
841:    SNESSetLineSearchCheck - Sets a routine to check the validity of new iterate computed
842:    by the line search routine in the Newton-based method SNESLS.

844:    Input Parameters:
845: +  snes - nonlinear context obtained from SNESCreate()
846: .  func - pointer to int function
847: -  checkctx - optional user-defined context for use by step checking routine 

849:    Collective on SNES

851:    Calling sequence of func:
852: .vb
853:    int func (SNES snes, void *checkctx, Vec x, PetscTruth *flag)
854: .ve
855:    where func returns an error code of 0 on success and a nonzero
856:    on failure.

858:    Input parameters for func:
859: +  snes - nonlinear context
860: .  checkctx - optional user-defined context for use by step checking routine 
861: -  x - current candidate iterate

863:    Output parameters for func:
864: +  x - current iterate (possibly modified)
865: -  flag - flag indicating whether x has been modified (either
866:            PETSC_TRUE of PETSC_FALSE)

868:    Level: advanced

870:    Notes:
871:    SNESNoLineSearch() and SNESNoLineSearchNoNorms() accept the new
872:    iterate computed by the line search checking routine.  In particular,
873:    these routines (1) compute a candidate iterate u_{i+1}, (2) pass control 
874:    to the checking routine, and then (3) compute the corresponding nonlinear
875:    function f(u_{i+1}) with the (possibly altered) iterate u_{i+1}.

877:    SNESQuadraticLineSearch() and SNESCubicLineSearch() also accept the
878:    new iterate computed by the line search checking routine.  In particular,
879:    these routines (1) compute a candidate iterate u_{i+1} as well as a
880:    candidate nonlinear function f(u_{i+1}), (2) pass control to the checking 
881:    routine, and then (3) force a re-evaluation of f(u_{i+1}) if any changes 
882:    were made to the candidate iterate in the checking routine (as indicated 
883:    by flag=PETSC_TRUE).  The overhead of this function re-evaluation can be
884:    very costly, so use this feature with caution!

886: .keywords: SNES, nonlinear, set, line search check, step check, routine

888: .seealso: SNESSetLineSearch()
889: @*/
890: int SNESSetLineSearchCheck(SNES snes,int (*func)(SNES,void*,Vec,PetscTruth*),void *checkctx)
891: {
892:   int ierr,(*f)(SNES,int (*)(SNES,void*,Vec,PetscTruth*),void*);

895:   PetscObjectQueryFunction((PetscObject)snes,"SNESSetLineSearchCheck_C",(void (**)(void))&f);
896:   if (f) {
897:     (*f)(snes,func,checkctx);
898:   }
899:   return(0);
900: }
901: /* -------------------------------------------------------------------------- */
902: EXTERN_C_BEGIN
903: #undef __FUNCT__  
905: int SNESSetLineSearchCheck_LS(SNES snes,int (*func)(SNES,void*,Vec,PetscTruth*),void *checkctx)
906: {
908:   ((SNES_LS *)(snes->data))->CheckStep = func;
909:   ((SNES_LS *)(snes->data))->checkP    = checkctx;
910:   return(0);
911: }
912: EXTERN_C_END
913: /* -------------------------------------------------------------------------- */
914: /*
915:    SNESPrintHelp_LS - Prints all options for the SNES_LS method.

917:    Input Parameter:
918: .  snes - the SNES context

920:    Application Interface Routine: SNESPrintHelp()
921: */
922: #undef __FUNCT__  
924: static int SNESPrintHelp_LS(SNES snes,char *p)
925: {
926:   SNES_LS *ls = (SNES_LS *)snes->data;

929:   (*PetscHelpPrintf)(snes->comm," method SNES_LS (ls) for systems of nonlinear equations:n");
930:   (*PetscHelpPrintf)(snes->comm,"   %ssnes_ls [cubic,quadratic,basic,basicnonorms,...]n",p);
931:   (*PetscHelpPrintf)(snes->comm,"   %ssnes_ls_alpha <alpha> (default %g)n",p,ls->alpha);
932:   (*PetscHelpPrintf)(snes->comm,"   %ssnes_ls_maxstep <max> (default %g)n",p,ls->maxstep);
933:   (*PetscHelpPrintf)(snes->comm,"   %ssnes_ls_steptol <tol> (default %g)n",p,ls->steptol);
934:   return(0);
935: }

937: /*
938:    SNESView_LS - Prints info from the SNESLS data structure.

940:    Input Parameters:
941: .  SNES - the SNES context
942: .  viewer - visualization context

944:    Application Interface Routine: SNESView()
945: */
946: #undef __FUNCT__  
948: static int SNESView_LS(SNES snes,PetscViewer viewer)
949: {
950:   SNES_LS    *ls = (SNES_LS *)snes->data;
951:   char       *cstr;
952:   int        ierr;
953:   PetscTruth isascii;

956:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
957:   if (isascii) {
958:     if (ls->LineSearch == SNESNoLineSearch)             cstr = "SNESNoLineSearch";
959:     else if (ls->LineSearch == SNESQuadraticLineSearch) cstr = "SNESQuadraticLineSearch";
960:     else if (ls->LineSearch == SNESCubicLineSearch)     cstr = "SNESCubicLineSearch";
961:     else                                                cstr = "unknown";
962:     PetscViewerASCIIPrintf(viewer,"  line search variant: %sn",cstr);
963:     PetscViewerASCIIPrintf(viewer,"  alpha=%g, maxstep=%g, steptol=%gn",ls->alpha,ls->maxstep,ls->steptol);
964:   } else {
965:     SETERRQ1(1,"Viewer type %s not supported for SNES EQ LS",((PetscObject)viewer)->type_name);
966:   }
967:   return(0);
968: }
969: /* -------------------------------------------------------------------------- */
970: /*
971:    SNESSetFromOptions_LS - Sets various parameters for the SNESLS method.

973:    Input Parameter:
974: .  snes - the SNES context

976:    Application Interface Routine: SNESSetFromOptions()
977: */
978: #undef __FUNCT__  
980: static int SNESSetFromOptions_LS(SNES snes)
981: {
982:   SNES_LS    *ls = (SNES_LS *)snes->data;
983:   char       ver[16],*lses[] = {"basic","basicnonorms","quadratic","cubic"};
984:   int        ierr;
985:   PetscTruth flg;

988:   PetscOptionsHead("SNES Line search options");
989:     PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",ls->alpha,&ls->alpha,0);
990:     PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",ls->maxstep,&ls->maxstep,0);
991:     PetscOptionsReal("-snes_ls_steptol","Step must be greater than","None",ls->steptol,&ls->steptol,0);

993:     PetscOptionsEList("-snes_ls","Line search used","SNESSetLineSearch",lses,4,"cubic",ver,16,&flg);
994:     if (flg) {
995:       PetscTruth isbasic,isnonorms,isquad,iscubic;

997:       PetscStrcmp(ver,lses[0],&isbasic);
998:       PetscStrcmp(ver,lses[1],&isnonorms);
999:       PetscStrcmp(ver,lses[2],&isquad);
1000:       PetscStrcmp(ver,lses[3],&iscubic);

1002:       if (isbasic) {
1003:         SNESSetLineSearch(snes,SNESNoLineSearch,PETSC_NULL);
1004:       } else if (isnonorms) {
1005:         SNESSetLineSearch(snes,SNESNoLineSearchNoNorms,PETSC_NULL);
1006:       } else if (isquad) {
1007:         SNESSetLineSearch(snes,SNESQuadraticLineSearch,PETSC_NULL);
1008:       } else if (iscubic) {
1009:         SNESSetLineSearch(snes,SNESCubicLineSearch,PETSC_NULL);
1010:       }
1011:       else {SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Unknown line search");}
1012:     }
1013:   PetscOptionsTail();
1014:   return(0);
1015: }
1016: /* -------------------------------------------------------------------------- */
1017: /*
1018:    SNESCreate_LS - Creates a nonlinear solver context for the SNESLS method,
1019:    SNES_LS, and sets this as the private data within the generic nonlinear solver
1020:    context, SNES, that was created within SNESCreate().


1023:    Input Parameter:
1024: .  snes - the SNES context

1026:    Application Interface Routine: SNESCreate()
1027:  */
1028: EXTERN_C_BEGIN
1029: #undef __FUNCT__  
1031: int SNESCreate_LS(SNES snes)
1032: {
1033:   int     ierr;
1034:   SNES_LS *neP;

1037:   snes->setup                = SNESSetUp_LS;
1038:   snes->solve                = SNESSolve_LS;
1039:   snes->destroy                = SNESDestroy_LS;
1040:   snes->converged        = SNESConverged_LS;
1041:   snes->printhelp       = SNESPrintHelp_LS;
1042:   snes->setfromoptions  = SNESSetFromOptions_LS;
1043:   snes->view            = SNESView_LS;
1044:   snes->nwork           = 0;

1046:   ierr                  = PetscNew(SNES_LS,&neP);
1047:   PetscLogObjectMemory(snes,sizeof(SNES_LS));
1048:   snes->data            = (void*)neP;
1049:   neP->alpha                = 1.e-4;
1050:   neP->maxstep                = 1.e8;
1051:   neP->steptol                = 1.e-12;
1052:   neP->LineSearch       = SNESCubicLineSearch;
1053:   neP->lsP              = PETSC_NULL;
1054:   neP->CheckStep        = PETSC_NULL;
1055:   neP->checkP           = PETSC_NULL;

1057:   PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESSetLineSearch_C","SNESSetLineSearch_LS",SNESSetLineSearch_LS);
1058:   PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESSetLineSearchCheck_C","SNESSetLineSearchCheck_LS",SNESSetLineSearchCheck_LS);

1060:   return(0);
1061: }
1062: EXTERN_C_END