Actual source code: snesmfjdef.c

      1: #define PETSCSNES_DLL

      3: /*
      4:   Implements the default PETSc approach for computing the h
      5:   parameter used with the finite difference based matrix-free
      6:   Jacobian-vector products.

      8:   To make your own: clone this file and modify for your needs.

     10:   Mandatory functions:
     11:   -------------------
     12:       MatSNESMFCompute_  - for a given point and direction computes h

     14:       MatSNESMFCreate_ - fills in the MatSNESMF data structure
     15:                            for this particular implementation

     17:
     18:    Optional functions:
     19:    -------------------
     20:       MatSNESMFView_ - prints information about the parameters being used.
     21:                        This is called when SNESView() or -snes_view is used.

     23:       MatSNESMFSetFromOptions_ - checks the options database for options that
     24:                                apply to this method.

     26:       MatSNESMFDestroy_ - frees any space allocated by the routines above

     28: */

     30: /*
     31:     This include file defines the data structure  MatSNESMF that
     32:    includes information about the computation of h. It is shared by
     33:    all implementations that people provide
     34: */
     35:  #include src/mat/matimpl.h
     36:  #include src/snes/mf/snesmfj.h

     38: /*
     39:       The default method has one parameter that is used to
     40:    "cutoff" very small values. This is stored in a data structure
     41:    that is only visible to this file. If your method has no parameters
     42:    it can omit this, if it has several simply reorganize the data structure.
     43:    The data structure is "hung-off" the MatSNESMF data structure in
     44:    the void *hctx; field.
     45: */
     46: typedef struct {
     47:   PetscReal umin;          /* minimum allowable u'a value relative to |u|_1 */
     48: } MatSNESMFDefault;

     52: /*
     53:    MatSNESMFCompute_Default - Standard PETSc code for computing the
     54:    differencing paramter (h) for use with matrix-free finite differences.

     56:    Input Parameters:
     57: +  ctx - the matrix free context
     58: .  U - the location at which you want the Jacobian
     59: -  a - the direction you want the derivative

     61:
     62:    Output Parameter:
     63: .  h - the scale computed

     65: */
     66: static PetscErrorCode MatSNESMFCompute_Default(MatSNESMFCtx ctx,Vec U,Vec a,PetscScalar *h,PetscTruth *zeroa)
     67: {
     68:   MatSNESMFDefault *hctx = (MatSNESMFDefault*)ctx->hctx;
     69:   PetscReal        nrm,sum,umin = hctx->umin;
     70:   PetscScalar      dot;
     71:   PetscErrorCode   ierr;

     74:   if (!(ctx->count % ctx->recomputeperiod)) {
     75:     /*
     76:      This algorithm requires 2 norms and 1 inner product. Rather than
     77:      use directly the VecNorm() and VecDot() routines (and thus have
     78:      three separate collective operations, we use the VecxxxBegin/End() routines
     79:     */
     80:     VecDotBegin(U,a,&dot);
     81:     VecNormBegin(a,NORM_1,&sum);
     82:     VecNormBegin(a,NORM_2,&nrm);
     83:     VecDotEnd(U,a,&dot);
     84:     VecNormEnd(a,NORM_1,&sum);
     85:     VecNormEnd(a,NORM_2,&nrm);

     87:     if (nrm == 0.0) {
     88:       *zeroa = PETSC_TRUE;
     89:       return(0);
     90:     }
     91:     *zeroa = PETSC_FALSE;

     93:     /*
     94:       Safeguard for step sizes that are "too small"
     95:     */
     96: #if defined(PETSC_USE_COMPLEX)
     97:     if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
     98:     else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
     99: #else
    100:     if (dot < umin*sum && dot >= 0.0) dot = umin*sum;
    101:     else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum;
    102: #endif
    103:     *h = ctx->error_rel*dot/(nrm*nrm);
    104:   } else {
    105:     *h = ctx->currenth;
    106:   }
    107:   if (*h != *h) SETERRQ3(PETSC_ERR_PLIB,"Differencing parameter is not a number sum = %g dot = %g norm = %g",sum,PetscRealPart(dot),nrm);
    108:   ctx->count++;
    109:   return(0);
    110: }

    114: /*
    115:    MatSNESMFView_Default - Prints information about this particular
    116:    method for computing h. Note that this does not print the general
    117:    information about the matrix-free method, as such info is printed
    118:    by the calling routine.

    120:    Input Parameters:
    121: +  ctx - the matrix free context
    122: -  viewer - the PETSc viewer
    123: */
    124: static PetscErrorCode MatSNESMFView_Default(MatSNESMFCtx ctx,PetscViewer viewer)
    125: {
    126:   MatSNESMFDefault *hctx = (MatSNESMFDefault *)ctx->hctx;
    127:   PetscErrorCode   ierr;
    128:   PetscTruth       iascii;

    131:   /*
    132:      Currently this only handles the ascii file viewers, others
    133:      could be added, but for this type of object other viewers
    134:      make less sense
    135:   */
    136:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    137:   if (iascii) {
    138:     PetscViewerASCIIPrintf(viewer,"    umin=%g (minimum iterate parameter)\n",hctx->umin);
    139:   } else {
    140:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for this SNES matrix free matrix",((PetscObject)viewer)->type_name);
    141:   }
    142:   return(0);
    143: }

    147: /*
    148:    MatSNESMFSetFromOptions_Default - Looks in the options database for
    149:    any options appropriate for this method.

    151:    Input Parameter:
    152: .  ctx - the matrix free context

    154: */
    155: static PetscErrorCode MatSNESMFSetFromOptions_Default(MatSNESMFCtx ctx)
    156: {
    157:   PetscErrorCode   ierr;
    158:   MatSNESMFDefault *hctx = (MatSNESMFDefault*)ctx->hctx;

    161:   PetscOptionsHead("Default matrix free parameters");
    162:     PetscOptionsReal("-snes_mf_umin","umin","MatSNESMFDefaultSetUmin",hctx->umin,&hctx->umin,0);
    163:   PetscOptionsTail();
    164:   return(0);
    165: }

    169: /*
    170:    MatSNESMFDestroy_Default - Frees the space allocated by
    171:    MatSNESMFCreate_Default().

    173:    Input Parameter:
    174: .  ctx - the matrix free context

    176:    Notes:
    177:    Does not free the ctx, that is handled by the calling routine
    178: */
    179: static PetscErrorCode MatSNESMFDestroy_Default(MatSNESMFCtx ctx)
    180: {

    184:   PetscFree(ctx->hctx);
    185:   return(0);
    186: }

    191: /*
    192:    The following two routines use the PetscObjectCompose() and PetscObjectQuery()
    193:    mechanism to allow the user to change the Umin parameter used in this method.
    194: */
    195: PetscErrorCode MatSNESMFDefaultSetUmin_Private(Mat mat,PetscReal umin)
    196: {
    197:   MatSNESMFCtx     ctx = (MatSNESMFCtx)mat->data;
    198:   MatSNESMFDefault *hctx;

    201:   if (!ctx) {
    202:     SETERRQ(PETSC_ERR_ARG_WRONG,"MatSNESMFDefaultSetUmin() attached to non-shell matrix");
    203:   }
    204:   hctx = (MatSNESMFDefault*)ctx->hctx;
    205:   hctx->umin = umin;
    206:   return(0);
    207: }

    212: /*@
    213:     MatSNESMFDefaultSetUmin - Sets the "umin" parameter used by the default
    214:     PETSc routine for computing the differencing parameter, h, which is used
    215:     for matrix-free Jacobian-vector products.

    217:    Input Parameters:
    218: +  A - the matrix created with MatCreateSNESMF()
    219: -  umin - the parameter

    221:    Level: advanced

    223:    Notes:
    224:    See the manual page for MatCreateSNESMF() for a complete description of the
    225:    algorithm used to compute h.

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

    229: @*/
    230: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFDefaultSetUmin(Mat A,PetscReal umin)
    231: {
    232:   PetscErrorCode ierr,(*f)(Mat,PetscReal);

    236:   PetscObjectQueryFunction((PetscObject)A,"MatSNESMFDefaultSetUmin_C",(void (**)(void))&f);
    237:   if (f) {
    238:     (*f)(A,umin);
    239:   }
    240:   return(0);
    241: }

    243: /*MC
    244:      MATSNESMF_DEFAULT - the default code for compute the "h" used in the finite difference
    245:             matrix-free matrix vector product

    247:    Options Database Keys:
    248: .  -snes_mf_umin <umin> see MatSNESMFDefaultSetUmin()

    251:    Level: intermediate

    253:    Notes: Requires 2 norms and 1 inner product, but they are computed together
    254:        so only one parallel collective operation is needed. See MATSNESMF_WP for a method
    255:        (with GMRES) that requires NO collective operations.

    257:    Formula used:
    258:      F'(u)*a = [F(u+h*a) - F(u)]/h where
    259:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
    260:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
    261:  where
    262:      error_rel = square root of relative error in function evaluation
    263:      umin = minimum iterate parameter

    265: .seealso: MATMFFD, MatCreateMF(), MatCreateSNESMF(), MATSNESMF_WP, MatSNESMFDefaultSetUmin()

    267: M*/
    271: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFCreate_Default(MatSNESMFCtx ctx)
    272: {
    273:   MatSNESMFDefault *hctx;
    274:   PetscErrorCode   ierr;

    278:   /* allocate my own private data structure */
    279:   PetscNew(MatSNESMFDefault,&hctx);
    280:   ctx->hctx  = (void*)hctx;
    281:   /* set a default for my parameter */
    282:   hctx->umin = 1.e-6;

    284:   /* set the functions I am providing */
    285:   ctx->ops->compute        = MatSNESMFCompute_Default;
    286:   ctx->ops->destroy        = MatSNESMFDestroy_Default;
    287:   ctx->ops->view           = MatSNESMFView_Default;
    288:   ctx->ops->setfromoptions = MatSNESMFSetFromOptions_Default;

    290:   PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFDefaultSetUmin_C",
    291:                             "MatSNESMFDefaultSetUmin_Private",
    292:                              MatSNESMFDefaultSetUmin_Private);
    293:   return(0);
    294: }