Actual source code: mffddef.c


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

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

  9:   Mandatory functions:
 10:   -------------------
 11:       MatMFFDCompute_  - for a given point and direction computes h

 13:       MatCreateMFFD _ - fills in the MatMFFD data structure
 14:                            for this particular implementation

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

 21:       MatMFFDSetFromOptions_ - checks the options database for options that
 22:                                apply to this method.

 24:       MatMFFDDestroy_ - frees any space allocated by the routines above

 26: */

 28: /*
 29:     This include file defines the data structure  MatMFFD that
 30:    includes information about the computation of h. It is shared by
 31:    all implementations that people provide
 32: */
 33: #include <petsc/private/matimpl.h>
 34: #include <../src/mat/impls/mffd/mffdimpl.h>

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

 48: /*
 49:    MatMFFDCompute_DS - Standard PETSc code for computing the
 50:    differencing parameter (h) for use with matrix-free finite differences.

 52:    Input Parameters:
 53: +  ctx - the matrix free context
 54: .  U - the location at which you want the Jacobian
 55: -  a - the direction you want the derivative

 57:    Output Parameter:
 58: .  h - the scale computed

 60: */
 61: static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx,Vec U,Vec a,PetscScalar *h,PetscBool  *zeroa)
 62: {
 63:   MatMFFD_DS     *hctx = (MatMFFD_DS*)ctx->hctx;
 64:   PetscReal      nrm,sum,umin = hctx->umin;
 65:   PetscScalar    dot;

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

 82:     if (nrm == 0.0) {
 83:       *zeroa = PETSC_TRUE;
 84:       return(0);
 85:     }
 86:     *zeroa = PETSC_FALSE;

 88:     /*
 89:       Safeguard for step sizes that are "too small"
 90:     */
 91:     if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
 92:     else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
 93:     *h = ctx->error_rel*dot/(nrm*nrm);
 94:     if (PetscIsInfOrNanScalar(*h)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Differencing parameter is not a number sum = %g dot = %g norm = %g",(double)sum,(double)PetscRealPart(dot),(double)nrm);
 95:   } else {
 96:     *h = ctx->currenth;
 97:   }
 98:   ctx->count++;
 99:   return(0);
100: }

102: /*
103:    MatMFFDView_DS - Prints information about this particular
104:    method for computing h. Note that this does not print the general
105:    information about the matrix-free method, as such info is printed
106:    by the calling routine.

108:    Input Parameters:
109: +  ctx - the matrix free context
110: -  viewer - the PETSc viewer
111: */
112: static PetscErrorCode MatMFFDView_DS(MatMFFD ctx,PetscViewer viewer)
113: {
114:   MatMFFD_DS     *hctx = (MatMFFD_DS*)ctx->hctx;
116:   PetscBool      iascii;

119:   /*
120:      Currently this only handles the ascii file viewers, others
121:      could be added, but for this type of object other viewers
122:      make less sense
123:   */
124:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
125:   if (iascii) {
126:     PetscViewerASCIIPrintf(viewer,"    umin=%g (minimum iterate parameter)\n",(double)hctx->umin);
127:   }
128:   return(0);
129: }

131: /*
132:    MatMFFDSetFromOptions_DS - Looks in the options database for
133:    any options appropriate for this method.

135:    Input Parameter:
136: .  ctx - the matrix free context

138: */
139: static PetscErrorCode MatMFFDSetFromOptions_DS(PetscOptionItems *PetscOptionsObject,MatMFFD ctx)
140: {
142:   MatMFFD_DS     *hctx = (MatMFFD_DS*)ctx->hctx;

145:   PetscOptionsHead(PetscOptionsObject,"Finite difference matrix free parameters");
146:   PetscOptionsReal("-mat_mffd_umin","umin","MatMFFDDSSetUmin",hctx->umin,&hctx->umin,NULL);
147:   PetscOptionsTail();
148:   return(0);
149: }

151: /*
152:    MatMFFDDestroy_DS - Frees the space allocated by
153:    MatCreateMFFD_DS().

155:    Input Parameter:
156: .  ctx - the matrix free context

158:    Notes:
159:    Does not free the ctx, that is handled by the calling routine
160: */
161: static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx)
162: {

166:   PetscFree(ctx->hctx);
167:   return(0);
168: }

170: /*
171:    The following two routines use the PetscObjectCompose() and PetscObjectQuery()
172:    mechanism to allow the user to change the Umin parameter used in this method.
173: */
174: PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat,PetscReal umin)
175: {
176:   MatMFFD        ctx=NULL;
177:   MatMFFD_DS     *hctx;

181:   MatShellGetContext(mat,&ctx);
182:   if (!ctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatMFFDDSSetUmin() attached to non-shell matrix");
183:   hctx       = (MatMFFD_DS*)ctx->hctx;
184:   hctx->umin = umin;
185:   return(0);
186: }

188: /*@
189:     MatMFFDDSSetUmin - Sets the "umin" parameter used by the
190:     PETSc routine for computing the differencing parameter, h, which is used
191:     for matrix-free Jacobian-vector products.

193:    Input Parameters:
194: +  A - the matrix created with MatCreateSNESMF()
195: -  umin - the parameter

197:    Level: advanced

199:    Notes:
200:    See the manual page for MatCreateSNESMF() for a complete description of the
201:    algorithm used to compute h.

203: .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF()

205: @*/
206: PetscErrorCode  MatMFFDDSSetUmin(Mat A,PetscReal umin)
207: {

212:   PetscTryMethod(A,"MatMFFDDSSetUmin_C",(Mat,PetscReal),(A,umin));
213:   return(0);
214: }

216: /*MC
217:      MATMFFD_DS - the code for compute the "h" used in the finite difference
218:             matrix-free matrix vector product.  This code
219:         implements the strategy in Dennis and Schnabel, "Numerical Methods for Unconstrained
220:         Optimization and Nonlinear Equations".

222:    Options Database Keys:
223: .  -mat_mffd_umin <umin> see MatMFFDDSSetUmin()

225:    Level: intermediate

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

232:    Formula used:
233:      F'(u)*a = [F(u+h*a) - F(u)]/h where
234:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
235:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
236:  where
237:      error_rel = square root of relative error in function evaluation
238:      umin = minimum iterate parameter

240: .seealso: MATMFFD, MatCreateMFFD(), MatCreateSNESMF(), MATMFFD_WP, MatMFFDDSSetUmin()

242: M*/
243: PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx)
244: {
245:   MatMFFD_DS     *hctx;

249:   /* allocate my own private data structure */
250:   PetscNewLog(ctx,&hctx);
251:   ctx->hctx = (void*)hctx;
252:   /* set a default for my parameter */
253:   hctx->umin = 1.e-6;

255:   /* set the functions I am providing */
256:   ctx->ops->compute        = MatMFFDCompute_DS;
257:   ctx->ops->destroy        = MatMFFDDestroy_DS;
258:   ctx->ops->view           = MatMFFDView_DS;
259:   ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS;

261:   PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDDSSetUmin_C",MatMFFDDSSetUmin_DS);
262:   return(0);
263: }