Actual source code: tr.c

      1: #define PETSCSNES_DLL
      2:
     3:  #include src/snes/impls/tr/tr.h

      5: /*
      6:    This convergence test determines if the two norm of the
      7:    solution lies outside the trust region, if so it halts.
      8: */
     11: PetscErrorCode SNES_TR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
     12: {
     13:   SNES                snes = (SNES) ctx;
     14:   SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx;
     15:   SNES_TR             *neP = (SNES_TR*)snes->data;
     16:   Vec                 x;
     17:   PetscReal           nrm;
     18:   PetscErrorCode      ierr;

     21:   if (snes->ksp_ewconv) {
     22:     if (!kctx) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Eisenstat-Walker onvergence context not created");
     23:     if (!n) {SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp);}
     24:     kctx->lresid_last = rnorm;
     25:   }
     26:   KSPDefaultConverged(ksp,n,rnorm,reason,ctx);
     27:   if (*reason) {
     28:     PetscLogInfo((snes,"SNES_TR_KSPConverged_Private: regular convergence test KSP iterations=%D, rnorm=%g\n",n,rnorm));
     29:   }

     31:   /* Determine norm of solution */
     32:   KSPBuildSolution(ksp,0,&x);
     33:   VecNorm(x,NORM_2,&nrm);
     34:   if (nrm >= neP->delta) {
     35:     PetscLogInfo((snes,"SNES_TR_KSPConverged_Private: KSP iterations=%D, rnorm=%g\n",n,rnorm));
     36:     PetscLogInfo((snes,"SNES_TR_KSPConverged_Private: Ending linear iteration early, delta=%g, length=%g\n",neP->delta,nrm));
     37:     *reason = KSP_CONVERGED_STEP_LENGTH;
     38:   }
     39:   return(0);
     40: }

     42: /*
     43:    SNESSolve_TR - Implements Newton's Method with a very simple trust
     44:    region approach for solving systems of nonlinear equations.

     46:
     47: */
     50: static PetscErrorCode SNESSolve_TR(SNES snes)
     51: {
     52:   SNES_TR             *neP = (SNES_TR*)snes->data;
     53:   Vec                 X,F,Y,G,TMP,Ytmp;
     54:   PetscErrorCode      ierr;
     55:   PetscInt            maxits,i,lits;
     56:   MatStructure        flg = DIFFERENT_NONZERO_PATTERN;
     57:   PetscReal           rho,fnorm,gnorm,gpnorm,xnorm,delta,nrm,ynorm,norm1;
     58:   PetscScalar         mone = -1.0,cnorm;
     59:   KSP                 ksp;
     60:   SNESConvergedReason reason;
     61:   PetscTruth          conv,breakout = PETSC_FALSE;

     64:   maxits        = snes->max_its;        /* maximum number of iterations */
     65:   X                = snes->vec_sol;        /* solution vector */
     66:   F                = snes->vec_func;        /* residual vector */
     67:   Y                = snes->work[0];        /* work vectors */
     68:   G                = snes->work[1];
     69:   Ytmp          = snes->work[2];

     71:   PetscObjectTakeAccess(snes);
     72:   snes->iter = 0;
     73:   PetscObjectGrantAccess(snes);
     74:   VecNorm(X,NORM_2,&xnorm);         /* xnorm = || X || */

     76:   SNESComputeFunction(snes,X,F);          /* F(X) */
     77:   VecNorm(F,NORM_2,&fnorm);             /* fnorm <- || F || */
     78:   PetscObjectTakeAccess(snes);
     79:   snes->norm = fnorm;
     80:   PetscObjectGrantAccess(snes);
     81:   delta = neP->delta0*fnorm;
     82:   neP->delta = delta;
     83:   SNESLogConvHistory(snes,fnorm,0);
     84:   SNESMonitor(snes,0,fnorm);
     85:   SNESGetKSP(snes,&ksp);

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

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

     92:   /* Set the stopping criteria to use the More' trick. */
     93:   PetscOptionsHasName(PETSC_NULL,"-snes_tr_ksp_regular_convergence_test",&conv);
     94:   if (!conv) {
     95:     KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void*)snes);
     96:     PetscLogInfo((snes,"SNESSolve_TR: Using Krylov convergence test SNES_TR_KSPConverged_Private\n"));
     97:   }
     98:
     99:   for (i=0; i<maxits; i++) {
    100:     SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
    101:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);

    103:     /* Solve J Y = F, where J is Jacobian matrix */
    104:     KSPSolve(snes->ksp,F,Ytmp);
    105:     KSPGetIterationNumber(ksp,&lits);
    106:     snes->linear_its += lits;
    107:     PetscLogInfo((snes,"SNESSolve_TR: iter=%D, linear solve iterations=%D\n",snes->iter,lits));
    108:     VecNorm(Ytmp,NORM_2,&nrm);
    109:     norm1 = nrm;
    110:     while(1) {
    111:       VecCopy(Ytmp,Y);
    112:       nrm = norm1;

    114:       /* Scale Y if need be and predict new value of F norm */
    115:       if (nrm >= delta) {
    116:         nrm = delta/nrm;
    117:         gpnorm = (1.0 - nrm)*fnorm;
    118:         cnorm = nrm;
    119:         PetscLogInfo((snes,"SNESSolve_TR: Scaling direction by %g\n",nrm));
    120:         VecScale(Y,cnorm);
    121:         nrm = gpnorm;
    122:         ynorm = delta;
    123:       } else {
    124:         gpnorm = 0.0;
    125:         PetscLogInfo((snes,"SNESSolve_TR: Direction is in Trust Region\n"));
    126:         ynorm = nrm;
    127:       }
    128:       VecAYPX(Y,mone,X);            /* Y <- X - Y */
    129:       VecCopy(X,snes->vec_sol_update_always);
    130:       SNESComputeFunction(snes,Y,G); /*  F(X) */
    131:       VecNorm(G,NORM_2,&gnorm);      /* gnorm <- || g || */
    132:       if (fnorm == gpnorm) rho = 0.0;
    133:       else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);

    135:       /* Update size of trust region */
    136:       if      (rho < neP->mu)  delta *= neP->delta1;
    137:       else if (rho < neP->eta) delta *= neP->delta2;
    138:       else                     delta *= neP->delta3;
    139:       PetscLogInfo((snes,"SNESSolve_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm));
    140:       PetscLogInfo((snes,"SNESSolve_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta));
    141:       neP->delta = delta;
    142:       if (rho > neP->sigma) break;
    143:       PetscLogInfo((snes,"SNESSolve_TR: Trying again in smaller region\n"));
    144:       /* check to see if progress is hopeless */
    145:       neP->itflag = PETSC_FALSE;
    146:       (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);
    147:       if (reason) {
    148:         /* We're not progressing, so return with the current iterate */
    149:         SNESMonitor(snes,i+1,fnorm);
    150:         breakout = PETSC_TRUE;
    151:         break;
    152:       }
    153:       snes->numFailures++;
    154:     }
    155:     if (!breakout) {
    156:       fnorm = gnorm;
    157:       PetscObjectTakeAccess(snes);
    158:       snes->iter = i+1;
    159:       snes->norm = fnorm;
    160:       PetscObjectGrantAccess(snes);
    161:       TMP = F; F = G; snes->vec_func_always = F; G = TMP;
    162:       TMP = X; X = Y; snes->vec_sol_always  = X; Y = TMP;
    163:       VecNorm(X,NORM_2,&xnorm);                /* xnorm = || X || */
    164:       SNESLogConvHistory(snes,fnorm,lits);
    165:       SNESMonitor(snes,i+1,fnorm);

    167:       /* Test for convergence */
    168:       neP->itflag = PETSC_TRUE;
    169:       (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);
    170:       if (reason) {
    171:         break;
    172:       }
    173:     } else {
    174:       break;
    175:     }
    176:   }
    177:   /* Verify solution is in corect location */
    178:   if (X != snes->vec_sol) {
    179:     VecCopy(X,snes->vec_sol);
    180:   }
    181:   if (F != snes->vec_func) {
    182:     VecCopy(F,snes->vec_func);
    183:   }
    184:   snes->vec_sol_always  = snes->vec_sol;
    185:   snes->vec_func_always = snes->vec_func;
    186:   if (i == maxits) {
    187:     PetscLogInfo((snes,"SNESSolve_TR: Maximum number of iterations has been reached: %D\n",maxits));
    188:     reason = SNES_DIVERGED_MAX_IT;
    189:   }
    190:   PetscObjectTakeAccess(snes);
    191:   snes->reason = reason;
    192:   PetscObjectGrantAccess(snes);
    193:   return(0);
    194: }
    195: /*------------------------------------------------------------*/
    198: static PetscErrorCode SNESSetUp_TR(SNES snes)
    199: {

    203:   snes->nwork = 4;
    204:   VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);
    205:   PetscLogObjectParents(snes,snes->nwork,snes->work);
    206:   snes->vec_sol_update_always = snes->work[3];
    207:   return(0);
    208: }
    209: /*------------------------------------------------------------*/
    212: static PetscErrorCode SNESDestroy_TR(SNES snes)
    213: {

    217:   if (snes->nwork) {
    218:     VecDestroyVecs(snes->work,snes->nwork);
    219:   }
    220:   PetscFree(snes->data);
    221:   return(0);
    222: }
    223: /*------------------------------------------------------------*/

    227: static PetscErrorCode SNESSetFromOptions_TR(SNES snes)
    228: {
    229:   SNES_TR *ctx = (SNES_TR *)snes->data;

    233:   PetscOptionsHead("SNES trust region options for nonlinear equations");
    234:     PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,0);
    235:     PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,0);
    236:     PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,0);
    237:     PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);
    238:     PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);
    239:     PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);
    240:     PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);
    241:     PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,0);
    242:   PetscOptionsTail();
    243:   return(0);
    244: }

    248: static PetscErrorCode SNESView_TR(SNES snes,PetscViewer viewer)
    249: {
    250:   SNES_TR *tr = (SNES_TR *)snes->data;
    252:   PetscTruth iascii;

    255:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    256:   if (iascii) {
    257:     PetscViewerASCIIPrintf(viewer,"  mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma);
    258:     PetscViewerASCIIPrintf(viewer,"  delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);
    259:   } else {
    260:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name);
    261:   }
    262:   return(0);
    263: }

    265: /* ---------------------------------------------------------------- */
    268: /*@C
    269:    SNESConverged_TR - Monitors the convergence of the trust region
    270:    method SNESTR for solving systems of nonlinear equations (default).

    272:    Collective on SNES

    274:    Input Parameters:
    275: +  snes - the SNES context
    276: .  xnorm - 2-norm of current iterate
    277: .  pnorm - 2-norm of current step
    278: .  fnorm - 2-norm of function
    279: -  dummy - unused context

    281:    Output Parameter:
    282: .   reason - one of
    283: $  SNES_CONVERGED_FNORM_ABS       - (fnorm < abstol),
    284: $  SNES_CONVERGED_PNORM_RELATIVE  - (pnorm < xtol*xnorm),
    285: $  SNES_CONVERGED_FNORM_RELATIVE  - (fnorm < rtol*fnorm0),
    286: $  SNES_DIVERGED_FUNCTION_COUNT   - (nfct > maxf),
    287: $  SNES_DIVERGED_FNORM_NAN        - (fnorm == NaN),
    288: $  SNES_CONVERGED_TR_DELTA        - (delta < xnorm*deltatol),
    289: $  SNES_CONVERGED_ITERATING       - (otherwise)

    291:    where
    292: +    delta    - trust region paramenter
    293: .    deltatol - trust region size tolerance,
    294:                 set with SNESSetTrustRegionTolerance()
    295: .    maxf - maximum number of function evaluations,
    296:             set with SNESSetTolerances()
    297: .    nfct - number of function evaluations,
    298: .    abstol - absolute function norm tolerance,
    299:             set with SNESSetTolerances()
    300: -    xtol - relative function norm tolerance,
    301:             set with SNESSetTolerances()

    303:    Level: intermediate

    305: .keywords: SNES, nonlinear, default, converged, convergence

    307: .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged()
    308: @*/
    309: PetscErrorCode PETSCSNES_DLLEXPORT SNESConverged_TR(SNES snes,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
    310: {
    311:   SNES_TR *neP = (SNES_TR *)snes->data;

    315:   if (fnorm != fnorm) {
    316:     PetscLogInfo((snes,"SNESConverged_TR:Failed to converged, function norm is NaN\n"));
    317:     *reason = SNES_DIVERGED_FNORM_NAN;
    318:   } else if (neP->delta < xnorm * snes->deltatol) {
    319:     PetscLogInfo((snes,"SNESConverged_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol));
    320:     *reason = SNES_CONVERGED_TR_DELTA;
    321:   } else if (neP->itflag) {
    322:     SNESConverged_LS(snes,xnorm,pnorm,fnorm,reason,dummy);
    323:   } else if (snes->nfuncs >= snes->max_funcs) {
    324:     PetscLogInfo((snes,"SNESConverged_TR: Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs));
    325:     *reason = SNES_DIVERGED_FUNCTION_COUNT;
    326:   } else {
    327:     *reason = SNES_CONVERGED_ITERATING;
    328:   }
    329:   return(0);
    330: }
    331: /* ------------------------------------------------------------ */
    332: /*MC
    333:       SNESTR - Newton based nonlinear solver that uses a trust region

    335:    Options Database:
    336: +    -snes_trtol <tol> Trust region tolerance
    337: .    -snes_tr_mu <mu>
    338: .    -snes_tr_eta <eta>
    339: .    -snes_tr_sigma <sigma>
    340: .    -snes_tr_delta0 <delta0>
    341: .    -snes_tr_delta1 <delta1>
    342: .    -snes_tr_delta2 <delta2>
    343: -    -snes_tr_delta3 <delta3>

    345:    The basic algorithm is taken from "The Minpack Project", by More',
    346:    Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
    347:    of Mathematical Software", Wayne Cowell, editor.

    349:    This is intended as a model implementation, since it does not
    350:    necessarily have many of the bells and whistles of other
    351:    implementations.

    353:    Level: intermediate

    355: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESLS, SNESSetTrustRegionTolerance()

    357: M*/
    361: PetscErrorCode PETSCSNES_DLLEXPORT SNESCreate_TR(SNES snes)
    362: {
    363:   SNES_TR        *neP;

    367:   snes->setup                = SNESSetUp_TR;
    368:   snes->solve                = SNESSolve_TR;
    369:   snes->destroy                = SNESDestroy_TR;
    370:   snes->converged        = SNESConverged_TR;
    371:   snes->setfromoptions  = SNESSetFromOptions_TR;
    372:   snes->view            = SNESView_TR;
    373:   snes->nwork           = 0;
    374:
    375:   ierr                        = PetscNew(SNES_TR,&neP);
    376:   PetscLogObjectMemory(snes,sizeof(SNES_TR));
    377:   snes->data                = (void*)neP;
    378:   neP->mu                = 0.25;
    379:   neP->eta                = 0.75;
    380:   neP->delta                = 0.0;
    381:   neP->delta0                = 0.2;
    382:   neP->delta1                = 0.3;
    383:   neP->delta2                = 0.75;
    384:   neP->delta3                = 2.0;
    385:   neP->sigma                = 0.0001;
    386:   neP->itflag                = PETSC_FALSE;
    387:   neP->rnorm0                = 0;
    388:   neP->ttol                = 0;
    389:   return(0);
    390: }