Actual source code: wp.c

      1: #define PETSCSNES_DLL

      3: /*MC
      4:      MATSNESMF_WP - Implements an alternative approach for computing the differencing parameter
      5:         h used with the finite difference based matrix-free Jacobian.  This code
      6:         implements the strategy of M. Pernice and H. Walker:

      8:       h = error_rel * sqrt(1 + ||U||) / ||a||

     10:       Notes:
     11:         1) || U || does not change between linear iterations so can be reused
     12:         2) In GMRES || a || == 1 and so does not need to ever be computed if you never
     13:            have a restart. Unfortunately a RESTART computes a matrix vector product
     14:            with ||a|| != 0 which breaks this

     16:       Reference:  M. Pernice and H. F. Walker, "NITSOL: A Newton Iterative
     17:       Solver for Nonlinear Systems", SIAM J. Sci. Stat. Comput.", 1998,
     18:       vol 19, pp. 302--318.

     20:    Options Database Keys:
     21: +   -snes_mf_compute_norma - compute the norm of a everytime see MatSNESMFWPSetComputeNormA()
     22: -   -snes_mf_compute_normu -Compute the norm of u everytime see MatSNESMFWPSetComputeNormU()

     25:    Level: intermediate

     27:    Notes: Requires no global collectives when used with GMRES

     29:    Formula used:
     30:      F'(u)*a = [F(u+h*a) - F(u)]/h where
     31:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
     32:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
     33:  where
     34:      error_rel = square root of relative error in function evaluation
     35:      umin = minimum iterate parameter

     37: .seealso: MATMFFD, MatCreateMF(), MatCreateSNESMF(), MATSNESMF_DEFAULT

     39: M*/

     41: /*
     42:     This include file defines the data structure  MatSNESMF that
     43:    includes information about the computation of h. It is shared by
     44:    all implementations that people provide.

     46:    See snesmfjdef.c for  a full set of comments on the routines below.
     47: */
     48:  #include src/mat/matimpl.h
     49:  #include src/snes/mf/snesmfj.h

     51: typedef struct {
     52:   PetscReal  normUfact;                   /* previous sqrt(1.0 + || U ||) */
     53:   PetscTruth computenorma,computenormU;
     54: } MatSNESMFWP;

     58: /*
     59:      MatSNESMFCompute_WP - Standard PETSc code for
     60:    computing h with matrix-free finite differences.

     62:   Input Parameters:
     63: +   ctx - the matrix free context
     64: .   U - the location at which you want the Jacobian
     65: -   a - the direction you want the derivative

     67:   Output Parameter:
     68: .   h - the scale computed

     70: */
     71: static PetscErrorCode MatSNESMFCompute_WP(MatSNESMFCtx ctx,Vec U,Vec a,PetscScalar *h,PetscTruth *zeroa)
     72: {
     73:   MatSNESMFWP    *hctx = (MatSNESMFWP*)ctx->hctx;
     74:   PetscReal      normU,norma = 1.0;

     79:   if (!(ctx->count % ctx->recomputeperiod)) {
     80:     if (hctx->computenorma && (hctx->computenormU || !ctx->ncurrenth)) {
     81:       VecNormBegin(U,NORM_2,&normU);
     82:       VecNormBegin(a,NORM_2,&norma);
     83:       VecNormEnd(U,NORM_2,&normU);
     84:       VecNormEnd(a,NORM_2,&norma);
     85:       hctx->normUfact = sqrt(1.0+normU);
     86:     } else if (hctx->computenormU || !ctx->ncurrenth) {
     87:       VecNorm(U,NORM_2,&normU);
     88:       hctx->normUfact = sqrt(1.0+normU);
     89:     } else if (hctx->computenorma) {
     90:       VecNorm(a,NORM_2,&norma);
     91:     }
     92:     if (norma == 0.0) {
     93:       *zeroa = PETSC_TRUE;
     94:       return(0);
     95:     }
     96:     *zeroa = PETSC_FALSE;
     97:     *h = ctx->error_rel*hctx->normUfact/norma;
     98:   } else {
     99:     *h = ctx->currenth;
    100:   }
    101:   return(0);
    102: }

    106: /*
    107:    MatSNESMFView_WP - Prints information about this particular
    108:      method for computing h. Note that this does not print the general
    109:      information about the matrix free, that is printed by the calling
    110:      routine.

    112:   Input Parameters:
    113: +   ctx - the matrix free context
    114: -   viewer - the PETSc viewer

    116: */
    117: static PetscErrorCode MatSNESMFView_WP(MatSNESMFCtx ctx,PetscViewer viewer)
    118: {
    119:   MatSNESMFWP    *hctx = (MatSNESMFWP *)ctx->hctx;
    121:   PetscTruth     iascii;

    124:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    125:   if (iascii) {
    126:     if (hctx->computenorma){PetscViewerASCIIPrintf(viewer,"    Computes normA\n");}
    127:     else                   { PetscViewerASCIIPrintf(viewer,"    Does not compute normA\n");}
    128:     if (hctx->computenormU){ PetscViewerASCIIPrintf(viewer,"    Computes normU\n");}
    129:     else                   { PetscViewerASCIIPrintf(viewer,"    Does not compute normU\n");}
    130:   } else {
    131:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES matrix-free WP",((PetscObject)viewer)->type_name);
    132:   }
    133:   return(0);
    134: }

    138: /*
    139:    MatSNESMFSetFromOptions_WP - Looks in the options database for
    140:      any options appropriate for this method

    142:   Input Parameter:
    143: .  ctx - the matrix free context

    145: */
    146: static PetscErrorCode MatSNESMFSetFromOptions_WP(MatSNESMFCtx ctx)
    147: {
    149:   MatSNESMFWP    *hctx = (MatSNESMFWP*)ctx->hctx;

    152:   PetscOptionsHead("Walker-Pernice options");
    153:     PetscOptionsTruth("-snes_mf_compute_norma","Compute the norm of a","MatSNESMFWPSetComputeNormA",
    154:                           hctx->computenorma,&hctx->computenorma,0);
    155:     PetscOptionsTruth("-snes_mf_compute_normu","Compute the norm of u","MatSNESMFWPSetComputeNormU",
    156:                           hctx->computenorma,&hctx->computenorma,0);
    157:   PetscOptionsTail();
    158:   return(0);
    159: }

    163: /*
    164:    MatSNESMFDestroy_WP - Frees the space allocated by
    165:        MatSNESMFCreate_WP().

    167:   Input Parameter:
    168: .  ctx - the matrix free context

    170:    Notes: does not free the ctx, that is handled by the calling routine

    172: */
    173: static PetscErrorCode MatSNESMFDestroy_WP(MatSNESMFCtx ctx)
    174: {
    177:   PetscFree(ctx->hctx);
    178:   return(0);
    179: }

    184: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFWPSetComputeNormA_P(Mat mat,PetscTruth flag)
    185: {
    186:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
    187:   MatSNESMFWP  *hctx;

    190:   hctx               = (MatSNESMFWP*)ctx->hctx;
    191:   hctx->computenorma = flag;
    192:   return(0);
    193: }

    198: /*@
    199:     MatSNESMFWPSetComputeNormA - Sets whether it computes the ||a|| used by the WP
    200:              PETSc routine for computing h. With GMRES since the ||a|| is always
    201:              one, you can save communication by setting this to false.

    203:   Input Parameters:
    204: +   A - the matrix created with MatCreateSNESMF()
    205: -   flag - PETSC_TRUE causes it to compute ||a||, PETSC_FALSE assumes it is 1.

    207:   Options Database Key:
    208: .   -snes_mf_compute_norma <true,false> - can be set false for GMRES to speed up code

    210:   Level: advanced

    212:   Notes:
    213:    See the manual page for MatCreateSNESMF() for a complete description of the
    214:    algorithm used to compute h.

    216: .seealso: MatSNESMFSetFunctionError(), MatCreateSNESMF()

    218: @*/
    219: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFWPSetComputeNormA(Mat A,PetscTruth flag)
    220: {
    221:   PetscErrorCode ierr,(*f)(Mat,PetscTruth);

    225:   PetscObjectQueryFunction((PetscObject)A,"MatSNESMFWPSetComputeNormA_C",(void (**)(void))&f);
    226:   if (f) {
    227:     (*f)(A,flag);
    228:   }
    229:   return(0);
    230: }

    235: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFWPSetComputeNormU_P(Mat mat,PetscTruth flag)
    236: {
    237:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
    238:   MatSNESMFWP  *hctx;

    241:   hctx               = (MatSNESMFWP*)ctx->hctx;
    242:   hctx->computenormU = flag;
    243:   return(0);
    244: }

    249: /*@
    250:     MatSNESMFWPSetComputeNormU - Sets whether it computes the ||U|| used by the WP
    251:              PETSc routine for computing h. With any Krylov solver this need only
    252:              be computed during the first iteration and kept for later.

    254:   Input Parameters:
    255: +   A - the matrix created with MatCreateSNESMF()
    256: -   flag - PETSC_TRUE causes it to compute ||U||, PETSC_FALSE uses the previous value

    258:   Options Database Key:
    259: .   -snes_mf_compute_normu <true,false> - true by default, false causes extra calculations

    261:   Level: advanced

    263:   Notes:
    264:    See the manual page for MatCreateSNESMF() for a complete description of the
    265:    algorithm used to compute h.

    267: .seealso: MatSNESMFSetFunctionError(), MatCreateSNESMF()

    269: @*/
    270: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFWPSetComputeNormU(Mat A,PetscTruth flag)
    271: {
    272:   PetscErrorCode ierr,(*f)(Mat,PetscTruth);

    276:   PetscObjectQueryFunction((PetscObject)A,"MatSNESMFWPSetComputeNormU_C",(void (**)(void))&f);
    277:   if (f) {
    278:     (*f)(A,flag);
    279:   }
    280:   return(0);
    281: }

    286: /*
    287:      MatSNESMFCreate_WP - Standard PETSc code for
    288:    computing h with matrix-free finite differences.

    290:    Input Parameter:
    291: .  ctx - the matrix free context created by MatSNESMFCreate()

    293: */
    294: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFCreate_WP(MatSNESMFCtx ctx)
    295: {
    297:   MatSNESMFWP    *hctx;

    301:   /* allocate my own private data structure */
    302:   PetscNew(MatSNESMFWP,&hctx);
    303:   ctx->hctx          = (void*)hctx;
    304:   hctx->computenormU = PETSC_FALSE;
    305:   hctx->computenorma = PETSC_TRUE;

    307:   /* set the functions I am providing */
    308:   ctx->ops->compute        = MatSNESMFCompute_WP;
    309:   ctx->ops->destroy        = MatSNESMFDestroy_WP;
    310:   ctx->ops->view           = MatSNESMFView_WP;
    311:   ctx->ops->setfromoptions = MatSNESMFSetFromOptions_WP;

    313:   PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFWPSetComputeNormA_C",
    314:                             "MatSNESMFWPSetComputeNormA_P",
    315:                              MatSNESMFWPSetComputeNormA_P);
    316:   PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFWPSetComputeNormU_C",
    317:                             "MatSNESMFWPSetComputeNormU_P",
    318:                              MatSNESMFWPSetComputeNormU_P);

    320:   return(0);
    321: }