Actual source code: ls.c

      1: #define PETSCSNES_DLL

     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: */
     13: PetscErrorCode SNESLSCheckLocalMin_Private(Mat A,Vec F,Vec W,PetscReal fnorm,PetscTruth *ismin)
     14: {
     15:   PetscReal      a1;
     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 minimum\n",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(W,PETSC_NULL);
     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 minimum\n",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: */
     51: PetscErrorCode SNESLSCheckResidual_Private(Mat A,Vec F,Vec X,Vec W1,Vec W2)
     52: {
     53:   PetscReal      a1,a2;
     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(W1,mone,F);

     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) {
     69:       PetscLogInfo((0,"SNESSolve_LS: ||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n",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 KSP, 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: */
    133: PetscErrorCode SNESSolve_LS(SNES snes)
    134: {
    135:   SNES_LS        *neP = (SNES_LS*)snes->data;
    137:   PetscInt       maxits,i,lits;
    138:   PetscTruth     lssucceed;
    139:   MatStructure   flg = DIFFERENT_NONZERO_PATTERN;
    140:   PetscReal      fnorm,gnorm,xnorm,ynorm;
    141:   Vec            Y,X,F,G,W,TMP;
    142:   KSP            ksp;

    145:   SNESGetKSP(snes,&ksp);
    146:   snes->reason  = SNES_CONVERGED_ITERATING;

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

    155:   PetscObjectTakeAccess(snes);
    156:   snes->iter = 0;
    157:   PetscObjectGrantAccess(snes);
    158:   SNESComputeFunction(snes,X,F);
    159:   VecNorm(F,NORM_2,&fnorm);        /* fnorm <- ||F||  */
    160:   if (fnorm != fnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
    161:   PetscObjectTakeAccess(snes);
    162:   snes->norm = fnorm;
    163:   PetscObjectGrantAccess(snes);
    164:   SNESLogConvHistory(snes,fnorm,0);
    165:   SNESMonitor(snes,0,fnorm);

    167:   if (fnorm < snes->abstol) {snes->reason = SNES_CONVERGED_FNORM_ABS; return(0);}

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

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

    174:     /* Call general purpose update function */
    175:     if (snes->update) {
    176:       (*snes->update)(snes, snes->iter);
    177:     }

    179:     /* Solve J Y = F, where J is Jacobian matrix */
    180:     SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
    181:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);
    182:     KSPSolve(snes->ksp,F,Y);
    183:     KSPGetIterationNumber(ksp,&lits);

    185:     if (neP->precheckstep) {
    186:       PetscTruth changed_y = PETSC_FALSE;
    187:       (*neP->precheckstep)(snes,X,Y,neP->precheck,&changed_y);
    188:     }

    190:     if (PetscLogPrintInfo){
    191:       SNESLSCheckResidual_Private(snes->jacobian,F,Y,G,W);
    192:     }

    194:     /* should check what happened to the linear solve? */
    195:     snes->linear_its += lits;
    196:     PetscLogInfo((snes,"SNESSolve_LS: iter=%D, linear solve iterations=%D\n",snes->iter,lits));

    198:     /* Compute a (scaled) negative update in the line search routine:
    199:          Y <- X - lambda*Y
    200:        and evaluate G(Y) = function(Y))
    201:     */
    202:     VecCopy(Y,snes->vec_sol_update_always);
    203:     (*neP->LineSearch)(snes,neP->lsP,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lssucceed);
    204:     PetscLogInfo((snes,"SNESSolve_LS: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed));
    205:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;

    207:     TMP = F; F = G; snes->vec_func_always = F; G = TMP;
    208:     TMP = X; X = W; snes->vec_sol_always = X;  W = TMP;
    209:     fnorm = gnorm;

    211:     PetscObjectTakeAccess(snes);
    212:     snes->iter = i+1;
    213:     snes->norm = fnorm;
    214:     PetscObjectGrantAccess(snes);
    215:     SNESLogConvHistory(snes,fnorm,lits);
    216:     SNESMonitor(snes,i+1,fnorm);

    218:     if (!lssucceed) {
    219:       PetscTruth ismin;

    221:       if (++snes->numFailures >= snes->maxFailures) {
    222:         snes->reason = SNES_DIVERGED_LS_FAILURE;
    223:         SNESLSCheckLocalMin_Private(snes->jacobian,F,W,fnorm,&ismin);
    224:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
    225:         break;
    226:       }
    227:     }

    229:     /* Test for convergence */
    230:     if (snes->converged) {
    231:       VecNorm(X,NORM_2,&xnorm);        /* xnorm = || X || */
    232:       (*snes->converged)(snes,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
    233:       if (snes->reason) {
    234:         break;
    235:       }
    236:     }
    237:   }
    238:   if (X != snes->vec_sol) {
    239:     VecCopy(X,snes->vec_sol);
    240:   }
    241:   if (F != snes->vec_func) {
    242:     VecCopy(F,snes->vec_func);
    243:   }
    244:   snes->vec_sol_always  = snes->vec_sol;
    245:   snes->vec_func_always = snes->vec_func;
    246:   if (i == maxits) {
    247:     PetscLogInfo((snes,"SNESSolve_LS: Maximum number of iterations has been reached: %D\n",maxits));
    248:     snes->reason = SNES_DIVERGED_MAX_IT;
    249:   }
    250:   return(0);
    251: }
    252: /* -------------------------------------------------------------------------- */
    253: /*
    254:    SNESSetUp_LS - Sets up the internal data structures for the later use
    255:    of the SNESLS nonlinear solver.

    257:    Input Parameter:
    258: .  snes - the SNES context
    259: .  x - the solution vector

    261:    Application Interface Routine: SNESSetUp()

    263:    Notes:
    264:    For basic use of the SNES solvers, the user need not explicitly call
    265:    SNESSetUp(), since these actions will automatically occur during
    266:    the call to SNESSolve().
    267:  */
    270: PetscErrorCode SNESSetUp_LS(SNES snes)
    271: {

    275:   if (!snes->work) {
    276:     snes->nwork = 4;
    277:     VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);
    278:     PetscLogObjectParents(snes,snes->nwork,snes->work);
    279:     snes->vec_sol_update_always = snes->work[3];
    280:   }
    281:   return(0);
    282: }
    283: /* -------------------------------------------------------------------------- */
    284: /*
    285:    SNESDestroy_LS - Destroys the private SNES_LS context that was created
    286:    with SNESCreate_LS().

    288:    Input Parameter:
    289: .  snes - the SNES context

    291:    Application Interface Routine: SNESDestroy()
    292:  */
    295: PetscErrorCode SNESDestroy_LS(SNES snes)
    296: {

    300:   if (snes->nwork) {
    301:     VecDestroyVecs(snes->work,snes->nwork);
    302:   }
    303:   PetscFree(snes->data);
    304:   return(0);
    305: }
    306: /* -------------------------------------------------------------------------- */

    310: /*@C
    311:    SNESLineSearchNo - This routine is not a line search at all;
    312:    it simply uses the full Newton step.  Thus, this routine is intended
    313:    to serve as a template and is not recommended for general use.

    315:    Collective on SNES and Vec

    317:    Input Parameters:
    318: +  snes - nonlinear context
    319: .  lsctx - optional context for line search (not used here)
    320: .  x - current iterate
    321: .  f - residual evaluated at x
    322: .  y - search direction
    323: .  w - work vector
    324: -  fnorm - 2-norm of f

    326:    Output Parameters:
    327: +  g - residual evaluated at new iterate y
    328: .  w - new iterate
    329: .  gnorm - 2-norm of g
    330: .  ynorm - 2-norm of search length
    331: -  flag - PETSC_TRUE on success, PETSC_FALSE on failure

    333:    Options Database Key:
    334: .  -snes_ls basic - Activates SNESLineSearchNo()

    336:    Level: advanced

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

    340: .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(),
    341:           SNESLineSearchSet(), SNESLineSearchNoNorms()
    342: @*/
    343: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchNo(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
    344: {
    346:   PetscScalar    mone = -1.0;
    347:   SNES_LS        *neP = (SNES_LS*)snes->data;
    348:   PetscTruth     changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;

    351:   *flag = PETSC_TRUE;
    352:   PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
    353:   VecNorm(y,NORM_2,ynorm);         /* ynorm = || y || */
    354:   VecWAXPY(w,mone,y,x);            /* w <- x - y   */
    355:   if (neP->postcheckstep) {
    356:    (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);
    357:   }
    358:   if (changed_y) {
    359:     VecWAXPY(w,mone,y,x);            /* w <- x - y   */
    360:   }
    361:   SNESComputeFunction(snes,w,g);
    362:   if (PetscExceptionValue(ierr)) {
    363:     PetscErrorCode pPetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr);
    364:   }
    365:

    367:   VecNorm(g,NORM_2,gnorm);  /* gnorm = || g || */
    368:   if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
    369:   PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
    370:   return(0);
    371: }
    372: /* -------------------------------------------------------------------------- */

    376: /*@C
    377:    SNESLineSearchNoNorms - This routine is not a line search at
    378:    all; it simply uses the full Newton step. This version does not
    379:    even compute the norm of the function or search direction; this
    380:    is intended only when you know the full step is fine and are
    381:    not checking for convergence of the nonlinear iteration (for
    382:    example, you are running always for a fixed number of Newton steps).

    384:    Collective on SNES and Vec

    386:    Input Parameters:
    387: +  snes - nonlinear context
    388: .  lsctx - optional context for line search (not used here)
    389: .  x - current iterate
    390: .  f - residual evaluated at x
    391: .  y - search direction
    392: .  w - work vector
    393: -  fnorm - 2-norm of f

    395:    Output Parameters:
    396: +  g - residual evaluated at new iterate y
    397: .  w - new iterate
    398: .  gnorm - not changed
    399: .  ynorm - not changed
    400: -  flag - set to PETSC_TRUE indicating a successful line search

    402:    Options Database Key:
    403: .  -snes_ls basicnonorms - Activates SNESLineSearchNoNorms()

    405:    Notes:
    406:    SNESLineSearchNoNorms() must be used in conjunction with
    407:    either the options
    408: $     -snes_no_convergence_test -snes_max_it <its>
    409:    or alternatively a user-defined custom test set via
    410:    SNESSetConvergenceTest(); or a -snes_max_it of 1,
    411:    otherwise, the SNES solver will generate an error.

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

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

    421:    Level: advanced

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

    425: .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(),
    426:           SNESLineSearchSet(), SNESLineSearchNo()
    427: @*/
    428: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchNoNorms(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
    429: {
    431:   PetscScalar    mone = -1.0;
    432:   SNES_LS        *neP = (SNES_LS*)snes->data;
    433:   PetscTruth     changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;

    436:   *flag = PETSC_TRUE;
    437:   PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
    438:   VecWAXPY(w,mone,y,x);            /* w <- x - y      */
    439:   if (neP->postcheckstep) {
    440:    (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);
    441:   }
    442:   if (changed_y) {
    443:     VecWAXPY(w,mone,y,x);            /* w <- x - y   */
    444:   }
    445:
    446:   /* don't evaluate function the last time through */
    447:   if (snes->iter < snes->max_its-1) {
    448:     SNESComputeFunction(snes,w,g);
    449:     if (PetscExceptionValue(ierr)) {
    450:       PetscErrorCode pPetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr);
    451:     }
    452:
    453:   }
    454:   PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
    455:   return(0);
    456: }
    457: /* -------------------------------------------------------------------------- */
    460: /*@C
    461:    SNESLineSearchCubic - Performs a cubic line search (default line search method).

    463:    Collective on SNES

    465:    Input Parameters:
    466: +  snes - nonlinear context
    467: .  lsctx - optional context for line search (not used here)
    468: .  x - current iterate
    469: .  f - residual evaluated at x
    470: .  y - search direction
    471: .  w - work vector
    472: -  fnorm - 2-norm of f

    474:    Output Parameters:
    475: +  g - residual evaluated at new iterate y
    476: .  w - new iterate
    477: .  gnorm - 2-norm of g
    478: .  ynorm - 2-norm of search length
    479: -  flag - PETSC_TRUE if line search succeeds; PETSC_FALSE on failure.

    481:    Options Database Key:
    482: $  -snes_ls cubic - Activates SNESLineSearchCubic()

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

    488:    Level: advanced

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

    492: .seealso: SNESLineSearchQuadratic(), SNESLineSearchNo(), SNESLineSearchSet(), SNESLineSearchNoNorms()
    493: @*/
    494: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchCubic(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
    495: {
    496:   /*
    497:      Note that for line search purposes we work with with the related
    498:      minimization problem:
    499:         min  z(x):  R^n -> R,
    500:      where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2.
    501:    */
    502:
    503:   PetscReal      steptol,initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength;
    504:   PetscReal      maxstep,minlambda,alpha,lambda,lambdatemp,lambdaneg;
    505: #if defined(PETSC_USE_COMPLEX)
    506:   PetscScalar    cinitslope,clambda;
    507: #endif
    509:   PetscInt       count;
    510:   SNES_LS        *neP = (SNES_LS*)snes->data;
    511:   PetscScalar    mone = -1.0,scale;
    512:   PetscTruth     changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;

    515:   PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
    516:   *flag   = PETSC_TRUE;
    517:   alpha   = neP->alpha;
    518:   maxstep = neP->maxstep;
    519:   steptol = neP->steptol;

    521:   VecNorm(y,NORM_2,ynorm);
    522:   if (!*ynorm) {
    523:     PetscLogInfo((snes,"SNESLineSearchCubic: Search direction and size is 0\n"));
    524:     *gnorm = fnorm;
    525:     VecCopy(x,w);
    526:     VecCopy(f,g);
    527:     *flag  = PETSC_FALSE;
    528:     goto theend1;
    529:   }
    530:   if (*ynorm > maxstep) {        /* Step too big, so scale back */
    531:     scale = maxstep/(*ynorm);
    532: #if defined(PETSC_USE_COMPLEX)
    533:     PetscLogInfo((snes,"SNESLineSearchCubic: Scaling step by %g old ynorm %g\n",PetscRealPart(scale),*ynorm));
    534: #else
    535:     PetscLogInfo((snes,"SNESLineSearchCubic: Scaling step by %g old ynorm %g\n",scale,*ynorm));
    536: #endif
    537:     VecScale(y,scale);
    538:     *ynorm = maxstep;
    539:   }
    540:   VecMaxPointwiseDivide(y,x,&rellength);
    541:   minlambda = steptol/rellength;
    542:   MatMult(snes->jacobian,y,w);
    543: #if defined(PETSC_USE_COMPLEX)
    544:   VecDot(f,w,&cinitslope);
    545:   initslope = PetscRealPart(cinitslope);
    546: #else
    547:   VecDot(f,w,&initslope);
    548: #endif
    549:   if (initslope > 0.0)  initslope = -initslope;
    550:   if (initslope == 0.0) initslope = -1.0;

    552:   VecCopy(y,w);
    553:   VecAYPX(w,mone,x);
    554:   if (snes->nfuncs >= snes->max_funcs) {
    555:     PetscLogInfo((snes,"SNESLineSearchCubic:Exceeded maximum function evaluations, while checking full step length!\n"));
    556:     *flag = PETSC_FALSE;
    557:     snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
    558:     goto theend1;
    559:   }
    560:   SNESComputeFunction(snes,w,g);
    561:   if (PetscExceptionValue(ierr)) {
    562:     PetscErrorCode pPetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr);
    563:   }
    564:
    565:   VecNorm(g,NORM_2,gnorm);
    566:   if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
    567:   if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */
    568:     PetscLogInfo((snes,"SNESLineSearchCubic: Using full step\n"));
    569:     goto theend1;
    570:   }

    572:   /* Fit points with quadratic */
    573:   lambda     = 1.0;
    574:   lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
    575:   lambdaprev = lambda;
    576:   gnormprev  = *gnorm;
    577:   if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
    578:   if (lambdatemp <= .1*lambda) lambda = .1*lambda;
    579:   else                         lambda = lambdatemp;
    580:   VecCopy(x,w);
    581:   lambdaneg = -lambda;
    582: #if defined(PETSC_USE_COMPLEX)
    583:   clambda   = lambdaneg; VecAXPY(w,clambda,y);
    584: #else
    585:   VecAXPY(w,lambdaneg,y);
    586: #endif
    587:   if (snes->nfuncs >= snes->max_funcs) {
    588:     PetscLogInfo((snes,"SNESLineSearchCubic:Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n"));
    589:     *flag = PETSC_FALSE;
    590:     snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
    591:     goto theend1;
    592:   }
    593:   SNESComputeFunction(snes,w,g);
    594:   if (PetscExceptionValue(ierr)) {
    595:     PetscErrorCode pPetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr);
    596:   }
    597:
    598:   VecNorm(g,NORM_2,gnorm);
    599:   if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
    600:   if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
    601:     PetscLogInfo((snes,"SNESLineSearchCubic: Quadratically determined step, lambda=%18.16e\n",lambda));
    602:     goto theend1;
    603:   }

    605:   /* Fit points with cubic */
    606:   count = 1;
    607:   while (count < 10000) {
    608:     if (lambda <= minlambda) { /* bad luck; use full step */
    609:       PetscLogInfo((snes,"SNESLineSearchCubic:Unable to find good step length! %D \n",count));
    610:       PetscLogInfo((snes,"SNESLineSearchCubic:fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope));
    611:       *flag = PETSC_FALSE;
    612:       break;
    613:     }
    614:     t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope;
    615:     t2 = .5*(gnormprev*gnormprev  - fnorm*fnorm) - lambdaprev*initslope;
    616:     a  = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
    617:     b  = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev);
    618:     d  = b*b - 3*a*initslope;
    619:     if (d < 0.0) d = 0.0;
    620:     if (a == 0.0) {
    621:       lambdatemp = -initslope/(2.0*b);
    622:     } else {
    623:       lambdatemp = (-b + sqrt(d))/(3.0*a);
    624:     }
    625:     lambdaprev = lambda;
    626:     gnormprev  = *gnorm;
    627:     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
    628:     if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
    629:     else                         lambda     = lambdatemp;
    630:     VecCopy(x,w);
    631:     lambdaneg = -lambda;
    632: #if defined(PETSC_USE_COMPLEX)
    633:     clambda   = lambdaneg;
    634:     VecAXPY(w,clambda,y);
    635: #else
    636:     VecAXPY(w,lambdaneg,y);
    637: #endif
    638:     if (snes->nfuncs >= snes->max_funcs) {
    639:       PetscLogInfo((snes,"SNESLineSearchCubic:Exceeded maximum function evaluations, while looking for good step length! %D \n",count));
    640:       PetscLogInfo((snes,"SNESLineSearchCubic:fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope));
    641:       *flag = PETSC_FALSE;
    642:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
    643:       break;
    644:     }
    645:     SNESComputeFunction(snes,w,g);
    646:     if (PetscExceptionValue(ierr)) {
    647:       PetscErrorCode pPetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr);
    648:     }
    649:
    650:     VecNorm(g,NORM_2,gnorm);
    651:     if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
    652:     if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* is reduction enough? */
    653:       PetscLogInfo((snes,"SNESLineSearchCubic: Cubically determined step, lambda=%18.16e\n",lambda));
    654:       break;
    655:     } else {
    656:       PetscLogInfo((snes,"SNESLineSearchCubic: Cubic step no good, shrinking lambda,  lambda=%18.16e\n",lambda));
    657:     }
    658:     count++;
    659:   }
    660:   if (count >= 10000) {
    661:     SETERRQ(PETSC_ERR_LIB, "Lambda was decreased more than 10,000 times, so something is probably wrong with the function evaluation");
    662:   }
    663:   theend1:
    664:   /* Optional user-defined check for line search step validity */
    665:   if (neP->postcheckstep && *flag) {
    666:     (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);
    667:     if (changed_y) {
    668:       VecWAXPY(w,mone,y,x);
    669:     }
    670:     if (changed_y || changed_w) { /* recompute the function if the step has changed */
    671:       SNESComputeFunction(snes,w,g);
    672:       if (PetscExceptionValue(ierr)) {
    673:         PetscErrorCode pPetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr);
    674:       }
    675:
    676:       VecNormBegin(g,NORM_2,gnorm);
    677:       if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
    678:       VecNormBegin(w,NORM_2,ynorm);
    679:       VecNormEnd(g,NORM_2,gnorm);
    680:       VecNormEnd(w,NORM_2,ynorm);
    681:     }
    682:   }
    683:   PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
    684:   return(0);
    685: }
    686: /* -------------------------------------------------------------------------- */
    689: /*@C
    690:    SNESLineSearchQuadratic - Performs a quadratic line search.

    692:    Collective on SNES and Vec

    694:    Input Parameters:
    695: +  snes - the SNES context
    696: .  lsctx - optional context for line search (not used here)
    697: .  x - current iterate
    698: .  f - residual evaluated at x
    699: .  y - search direction
    700: .  w - work vector
    701: -  fnorm - 2-norm of f

    703:    Output Parameters:
    704: +  g - residual evaluated at new iterate y
    705: .  w - new iterate
    706: .  gnorm - 2-norm of g
    707: .  ynorm - 2-norm of search length
    708: -  flag - PETSC_TRUE if line search succeeds; PETSC_FALSE on failure.

    710:    Options Database Key:
    711: .  -snes_ls quadratic - Activates SNESLineSearchQuadratic()

    713:    Notes:
    714:    Use SNESLineSearchSet() to set this routine within the SNESLS method.

    716:    Level: advanced

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

    720: .seealso: SNESLineSearchCubic(), SNESLineSearchNo(), SNESLineSearchSet(), SNESLineSearchNoNorms()
    721: @*/
    722: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchQuadratic(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
    723: {
    724:   /*
    725:      Note that for line search purposes we work with with the related
    726:      minimization problem:
    727:         min  z(x):  R^n -> R,
    728:      where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2.
    729:    */
    730:   PetscReal      steptol,initslope,maxstep,minlambda,alpha,lambda,lambdatemp,lambdaneg,rellength;
    731: #if defined(PETSC_USE_COMPLEX)
    732:   PetscScalar    cinitslope,clambda;
    733: #endif
    735:   PetscInt       count;
    736:   SNES_LS        *neP = (SNES_LS*)snes->data;
    737:   PetscScalar    mone = -1.0,scale;
    738:   PetscTruth     changed_w = PETSC_FALSE,changed_y = PETSC_FALSE;

    741:   PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);
    742:   *flag   = PETSC_TRUE;
    743:   alpha   = neP->alpha;
    744:   maxstep = neP->maxstep;
    745:   steptol = neP->steptol;

    747:   VecNorm(y,NORM_2,ynorm);
    748:   if (*ynorm == 0.0) {
    749:     PetscLogInfo((snes,"SNESLineSearchQuadratic: Search direction and size is 0\n"));
    750:     *gnorm = fnorm;
    751:     VecCopy(x,y);
    752:     VecCopy(f,g);
    753:     *flag  = PETSC_FALSE;
    754:     goto theend2;
    755:   }
    756:   if (*ynorm > maxstep) {        /* Step too big, so scale back */
    757:     scale  = maxstep/(*ynorm);
    758:     VecScale(y,scale);
    759:     *ynorm = maxstep;
    760:   }
    761:   VecMaxPointwiseDivide(y,x,&rellength);
    762:   minlambda = steptol/rellength;
    763:   MatMult(snes->jacobian,y,w);
    764: #if defined(PETSC_USE_COMPLEX)
    765:   VecDot(f,w,&cinitslope);
    766:   initslope = PetscRealPart(cinitslope);
    767: #else
    768:   VecDot(f,w,&initslope);
    769: #endif
    770:   if (initslope > 0.0)  initslope = -initslope;
    771:   if (initslope == 0.0) initslope = -1.0;

    773:   VecCopy(y,w);
    774:   VecAYPX(w,mone,x);
    775:   if (snes->nfuncs >= snes->max_funcs) {
    776:     PetscLogInfo((snes,"SNESLineSearchQuadratic:Exceeded maximum function evaluations, while checking full step length!\n"));
    777:     VecCopy(w,y);
    778:     *flag = PETSC_FALSE;
    779:     snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
    780:     goto theend2;
    781:   }
    782:   SNESComputeFunction(snes,w,g);
    783:   if (PetscExceptionValue(ierr)) {
    784:     PetscErrorCode pPetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr);
    785:   }
    786:
    787:   VecNorm(g,NORM_2,gnorm);
    788:   if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
    789:   if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + alpha*initslope) { /* Sufficient reduction */
    790:     VecCopy(w,y);
    791:     PetscLogInfo((snes,"SNESLineSearchQuadratic: Using full step\n"));
    792:     goto theend2;
    793:   }

    795:   /* Fit points with quadratic */
    796:   lambda = 1.0;
    797:   count = 1;
    798:   while (PETSC_TRUE) {
    799:     if (lambda <= minlambda) { /* bad luck; use full step */
    800:       PetscLogInfo((snes,"SNESLineSearchQuadratic:Unable to find good step length! %D \n",count));
    801:       PetscLogInfo((snes,"SNESLineSearchQuadratic:fnorm=%g, gnorm=%g, ynorm=%g, lambda=%g, initial slope=%g\n",fnorm,*gnorm,*ynorm,lambda,initslope));
    802:       VecCopy(x,y);
    803:       *flag = PETSC_FALSE;
    804:       break;
    805:     }
    806:     lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope);
    807:     if (lambdatemp > .5*lambda)  lambdatemp = .5*lambda;
    808:     if (lambdatemp <= .1*lambda) lambda     = .1*lambda;
    809:     else                         lambda     = lambdatemp;
    810:     VecCopy(x,w);
    811:     lambdaneg = -lambda;
    812: #if defined(PETSC_USE_COMPLEX)
    813:     clambda   = lambdaneg; VecAXPY(w,clambda,y);
    814: #else
    815:     VecAXPY(w,lambdaneg,y);
    816: #endif
    817:     if (snes->nfuncs >= snes->max_funcs) {
    818:       PetscLogInfo((snes,"SNESLineSearchQuadratic:Exceeded maximum function evaluations, while looking for good step length! %D \n",count));
    819:       PetscLogInfo((snes,"SNESLineSearchQuadratic:fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope));
    820:       VecCopy(w,y);
    821:       *flag = PETSC_FALSE;
    822:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
    823:       break;
    824:     }
    825:     SNESComputeFunction(snes,w,g);
    826:     if (PetscExceptionValue(ierr)) {
    827:       PetscErrorCode pPetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr);
    828:     }
    829:
    830:     VecNorm(g,NORM_2,gnorm);
    831:     if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
    832:     if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*alpha*initslope) { /* sufficient reduction */
    833:       VecCopy(w,y);
    834:       PetscLogInfo((snes,"SNESLineSearchQuadratic:Quadratically determined step, lambda=%g\n",lambda));
    835:       break;
    836:     }
    837:     count++;
    838:   }
    839:   theend2:
    840:   /* Optional user-defined check for line search step validity */
    841:   if (neP->postcheckstep) {
    842:     (*neP->postcheckstep)(snes,x,y,w,neP->postcheck,&changed_y,&changed_w);
    843:     if (changed_y) {
    844:       VecWAXPY(w,mone,y,x);
    845:     }
    846:     if (changed_y || changed_w) { /* recompute the function if the step has changed */
    847:       SNESComputeFunction(snes,w,g);
    848:       if (PetscExceptionValue(ierr)) {
    849:         PetscErrorCode pPetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(pierr);
    850:       }
    851:
    852:       VecNormBegin(g,NORM_2,gnorm);
    853:       VecNormBegin(w,NORM_2,ynorm);
    854:       VecNormEnd(g,NORM_2,gnorm);
    855:       VecNormEnd(w,NORM_2,ynorm);
    856:       if (*gnorm != *gnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
    857:     }
    858:   }
    859:   PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);
    860:   return(0);
    861: }

    863: /* -------------------------------------------------------------------------- */
    866: /*@C
    867:    SNESLineSearchSet - Sets the line search routine to be used
    868:    by the method SNESLS.

    870:    Input Parameters:
    871: +  snes - nonlinear context obtained from SNESCreate()
    872: .  lsctx - optional user-defined context for use by line search
    873: -  func - pointer to int function

    875:    Collective on SNES

    877:    Available Routines:
    878: +  SNESLineSearchCubic() - default line search
    879: .  SNESLineSearchQuadratic() - quadratic line search
    880: .  SNESLineSearchNo() - the full Newton step (actually not a line search)
    881: -  SNESLineSearchNoNorms() - the full Newton step (calculating no norms; faster in parallel)

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

    891:    Calling sequence of func:
    892: .vb
    893:    func (SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,
    894:          PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag)
    895: .ve

    897:     Input parameters for func:
    898: +   snes - nonlinear context
    899: .   lsctx - optional user-defined context for line search
    900: .   x - current iterate
    901: .   f - residual evaluated at x
    902: .   y - search direction
    903: .   w - work vector
    904: -   fnorm - 2-norm of f

    906:     Output parameters for func:
    907: +   g - residual evaluated at new iterate y
    908: .   w - new iterate
    909: .   gnorm - 2-norm of g
    910: .   ynorm - 2-norm of search length
    911: -   flag - set to PETSC_TRUE if the line search succeeds; PETSC_FALSE on failure.

    913:     Level: advanced

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

    917: .seealso: SNESLineSearchCubic(), SNESLineSearchQuadratic(), SNESLineSearchNo(), SNESLineSearchNoNorms(),
    918:           SNESLineSearchSetPostCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams(), SNESLineSearchSetPreCheck()
    919: @*/
    920: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSet(SNES snes,PetscErrorCode (*func)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,PetscTruth*),void *lsctx)
    921: {
    922:   PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal*,PetscReal*,PetscTruth*),void*);

    925:   PetscObjectQueryFunction((PetscObject)snes,"SNESLineSearchSet_C",(void (**)(void))&f);
    926:   if (f) {
    927:     (*f)(snes,func,lsctx);
    928:   }
    929:   return(0);
    930: }

    933: /* -------------------------------------------------------------------------- */
    937: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSet_LS(SNES snes,FCN2 func,void *lsctx)
    938: {
    940:   ((SNES_LS *)(snes->data))->LineSearch = func;
    941:   ((SNES_LS *)(snes->data))->lsP        = lsctx;
    942:   return(0);
    943: }
    945: /* -------------------------------------------------------------------------- */
    948: /*@C
    949:    SNESLineSearchSetPostCheck - Sets a routine to check the validity of new iterate computed
    950:    by the line search routine in the Newton-based method SNESLS.

    952:    Input Parameters:
    953: +  snes - nonlinear context obtained from SNESCreate()
    954: .  func - pointer to function
    955: -  checkctx - optional user-defined context for use by step checking routine

    957:    Collective on SNES

    959:    Calling sequence of func:
    960: .vb
    961:    int func (SNES snes, Vec x,Vec y,Vec w,void *checkctx, PetscTruth *changed_y,PetscTruth *changed_w)
    962: .ve
    963:    where func returns an error code of 0 on success and a nonzero
    964:    on failure.

    966:    Input parameters for func:
    967: +  snes - nonlinear context
    968: .  checkctx - optional user-defined context for use by step checking routine
    969: .  x - previous iterate
    970: .  y - new search direction and length
    971: -  w - current candidate iterate

    973:    Output parameters for func:
    974: +  y - search direction (possibly changed)
    975: .  w - current iterate (possibly modified)
    976: .  changed_y - indicates search direction was changed by this routine
    977: -  changed_w - indicates current iterate was changed by this routine

    979:    Level: advanced

    981:    Notes: All line searches accept the new iterate computed by the line search checking routine.

    983:    Only one of changed_y and changed_w can  be PETSC_TRUE

    985:    On input w = x + y

    987:    SNESLineSearchNo() and SNESLineSearchNoNorms() (1) compute a candidate iterate u_{i+1}, (2) pass control
    988:    to the checking routine, and then (3) compute the corresponding nonlinear
    989:    function f(u_{i+1}) with the (possibly altered) iterate u_{i+1}.

    991:    SNESLineSearchQuadratic() and SNESLineSearchCubic() (1) compute a candidate iterate u_{i+1} as well as a
    992:    candidate nonlinear function f(u_{i+1}), (2) pass control to the checking
    993:    routine, and then (3) force a re-evaluation of f(u_{i+1}) if any changes
    994:    were made to the candidate iterate in the checking routine (as indicated
    995:    by flag=PETSC_TRUE).  The overhead of this extra function re-evaluation can be
    996:    very costly, so use this feature with caution!

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

    1000: .seealso: SNESLineSearchSet(), SNESLineSearchSetPreCheck()
    1001: @*/
    1002: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetPostCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,Vec,void*,PetscTruth*,PetscTruth*),void *checkctx)
    1003: {
    1004:   PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec,Vec,void*,PetscTruth*,PetscTruth*),void*);

    1007:   PetscObjectQueryFunction((PetscObject)snes,"SNESLineSearchSetPostCheck_C",(void (**)(void))&f);
    1008:   if (f) {
    1009:     (*f)(snes,func,checkctx);
    1010:   }
    1011:   return(0);
    1012: }
    1015: /*@C
    1016:    SNESLineSearchSetPreCheck - Sets a routine to check the validity of a new direction given by the linear solve

    1018:    Input Parameters:
    1019: +  snes - nonlinear context obtained from SNESCreate()
    1020: .  func - pointer to function
    1021: -  checkctx - optional user-defined context for use by step checking routine

    1023:    Collective on SNES

    1025:    Calling sequence of func:
    1026: .vb
    1027:    int func (SNES snes, Vec x,Vec y,,void *checkctx, PetscTruth *changed_y)
    1028: .ve
    1029:    where func returns an error code of 0 on success and a nonzero
    1030:    on failure.

    1032:    Input parameters for func:
    1033: +  snes - nonlinear context
    1034: .  checkctx - optional user-defined context for use by step checking routine
    1035: .  x - previous iterate
    1036: -  y - new search direction and length

    1038:    Output parameters for func:
    1039: +  y - search direction (possibly changed)
    1040: -  changed_y - indicates search direction was changed by this routine

    1042:    Level: advanced

    1044:    Notes: All line searches accept the new iterate computed by the line search checking routine.

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

    1048: .seealso: SNESLineSearchSet(), SNESLineSearchSetPostCheck()
    1049: @*/
    1050: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetPreCheck(SNES snes,PetscErrorCode (*func)(SNES,Vec,Vec,void*,PetscTruth*),void *checkctx)
    1051: {
    1052:   PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec,void*,PetscTruth*),void*);

    1055:   PetscObjectQueryFunction((PetscObject)snes,"SNESLineSearchSetPreCheck_C",(void (**)(void))&f);
    1056:   if (f) {
    1057:     (*f)(snes,func,checkctx);
    1058:   }
    1059:   return(0);
    1060: }

    1062: /* -------------------------------------------------------------------------- */
    1068: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetPostCheck_LS(SNES snes,FCN1 func,void *checkctx)
    1069: {
    1071:   ((SNES_LS *)(snes->data))->postcheckstep = func;
    1072:   ((SNES_LS *)(snes->data))->postcheck     = checkctx;
    1073:   return(0);
    1074: }

    1080: PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetPreCheck_LS(SNES snes,FCN3 func,void *checkctx)
    1081: {
    1083:   ((SNES_LS *)(snes->data))->precheckstep = func;
    1084:   ((SNES_LS *)(snes->data))->precheck     = checkctx;
    1085:   return(0);
    1086: }
    1088: /* -------------------------------------------------------------------------- */
    1089: /*
    1090:    SNESPrintHelp_LS - Prints all options for the SNES_LS method.

    1092:    Input Parameter:
    1093: .  snes - the SNES context

    1095:    Application Interface Routine: SNESPrintHelp()
    1096: */
    1099: static PetscErrorCode SNESPrintHelp_LS(SNES snes,char *p)
    1100: {
    1101:   SNES_LS *ls = (SNES_LS *)snes->data;

    1104:   (*PetscHelpPrintf)(snes->comm," method SNES_LS (ls) for systems of nonlinear equations:\n");
    1105:   (*PetscHelpPrintf)(snes->comm,"   %ssnes_ls [cubic,quadratic,basic,basicnonorms,...]\n",p);
    1106:   (*PetscHelpPrintf)(snes->comm,"   %ssnes_ls_alpha <alpha> (default %g)\n",p,ls->alpha);
    1107:   (*PetscHelpPrintf)(snes->comm,"   %ssnes_ls_maxstep <max> (default %g)\n",p,ls->maxstep);
    1108:   (*PetscHelpPrintf)(snes->comm,"   %ssnes_ls_steptol <tol> (default %g)\n",p,ls->steptol);
    1109:   return(0);
    1110: }

    1112: /*
    1113:    SNESView_LS - Prints info from the SNESLS data structure.

    1115:    Input Parameters:
    1116: .  SNES - the SNES context
    1117: .  viewer - visualization context

    1119:    Application Interface Routine: SNESView()
    1120: */
    1123: static PetscErrorCode SNESView_LS(SNES snes,PetscViewer viewer)
    1124: {
    1125:   SNES_LS        *ls = (SNES_LS *)snes->data;
    1126:   const char     *cstr;
    1128:   PetscTruth     iascii;

    1131:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    1132:   if (iascii) {
    1133:     if (ls->LineSearch == SNESLineSearchNo)             cstr = "SNESLineSearchNo";
    1134:     else if (ls->LineSearch == SNESLineSearchQuadratic) cstr = "SNESLineSearchQuadratic";
    1135:     else if (ls->LineSearch == SNESLineSearchCubic)     cstr = "SNESLineSearchCubic";
    1136:     else                                                cstr = "unknown";
    1137:     PetscViewerASCIIPrintf(viewer,"  line search variant: %s\n",cstr);
    1138:     PetscViewerASCIIPrintf(viewer,"  alpha=%g, maxstep=%g, steptol=%g\n",ls->alpha,ls->maxstep,ls->steptol);
    1139:   } else {
    1140:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ LS",((PetscObject)viewer)->type_name);
    1141:   }
    1142:   return(0);
    1143: }
    1144: /* -------------------------------------------------------------------------- */
    1145: /*
    1146:    SNESSetFromOptions_LS - Sets various parameters for the SNESLS method.

    1148:    Input Parameter:
    1149: .  snes - the SNES context

    1151:    Application Interface Routine: SNESSetFromOptions()
    1152: */
    1155: static PetscErrorCode SNESSetFromOptions_LS(SNES snes)
    1156: {
    1157:   SNES_LS        *ls = (SNES_LS *)snes->data;
    1158:   const char     *lses[] = {"basic","basicnonorms","quadratic","cubic"};
    1160:   PetscInt       indx;
    1161:   PetscTruth     flg;

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

    1169:     PetscOptionsEList("-snes_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);
    1170:     if (flg) {
    1171:       switch (indx) {
    1172:       case 0:
    1173:         SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);
    1174:         break;
    1175:       case 1:
    1176:         SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);
    1177:         break;
    1178:       case 2:
    1179:         SNESLineSearchSet(snes,SNESLineSearchQuadratic,PETSC_NULL);
    1180:         break;
    1181:       case 3:
    1182:         SNESLineSearchSet(snes,SNESLineSearchCubic,PETSC_NULL);
    1183:         break;
    1184:       }
    1185:     }
    1186:   PetscOptionsTail();
    1187:   return(0);
    1188: }
    1189: /* -------------------------------------------------------------------------- */
    1190: /*MC
    1191:       SNESLS - Newton based nonlinear solver that uses a line search

    1193:    Options Database:
    1194: +   -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search
    1195: .   -snes_ls_alpha <alpha> - Sets alpha
    1196: .   -snes_ls_maxstep <max> - Sets maxstep
    1197: -   -snes_ls_steptol <steptol> - Sets steptol, this is the minimum step size that the line search code
    1198:                    will accept; min p[i]/x[i] < steptol. The -snes_stol <stol> is the minimum step length
    1199:                    the default convergence test will use and is based on 2-norm(p) < stol*2-norm(x)

    1201:     Notes: This is the default nonlinear solver in SNES

    1203:    Level: beginner

    1205: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(),
    1206:            SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(),
    1207:           SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck()

    1209: M*/
    1213: PetscErrorCode PETSCSNES_DLLEXPORT SNESCreate_LS(SNES snes)
    1214: {
    1216:   SNES_LS        *neP;

    1219:   snes->setup                = SNESSetUp_LS;
    1220:   snes->solve                = SNESSolve_LS;
    1221:   snes->destroy                = SNESDestroy_LS;
    1222:   snes->converged        = SNESConverged_LS;
    1223:   snes->printhelp       = SNESPrintHelp_LS;
    1224:   snes->setfromoptions  = SNESSetFromOptions_LS;
    1225:   snes->view            = SNESView_LS;
    1226:   snes->nwork           = 0;

    1228:   PetscNew(SNES_LS,&neP);
    1229:   PetscLogObjectMemory(snes,sizeof(SNES_LS));
    1230:   snes->data            = (void*)neP;
    1231:   neP->alpha                = 1.e-4;
    1232:   neP->maxstep                = 1.e8;
    1233:   neP->steptol                = 1.e-12;
    1234:   neP->LineSearch       = SNESLineSearchCubic;
    1235:   neP->lsP              = PETSC_NULL;
    1236:   neP->postcheckstep    = PETSC_NULL;
    1237:   neP->postcheck        = PETSC_NULL;
    1238:   neP->precheckstep     = PETSC_NULL;
    1239:   neP->precheck         = PETSC_NULL;

    1241:   PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_LS",SNESLineSearchSet_LS);
    1242:   PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetPostCheck_C","SNESLineSearchSetPostCheck_LS",SNESLineSearchSetPostCheck_LS);
    1243:   PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetPreCheck_C","SNESLineSearchSetPreCheck_LS",SNESLineSearchSetPreCheck_LS);

    1245:   return(0);
    1246: }