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;

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

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

 86:     /*
 87:       Safeguard for step sizes that are "too small"
 88:     */
 89:     if (PetscAbsScalar(dot) < umin * sum && PetscRealPart(dot) >= 0.0) dot = umin * sum;
 90:     else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin * sum) dot = -umin * sum;
 91:     *h = ctx->error_rel * dot / (nrm * nrm);
 93:   } else {
 94:     *h = ctx->currenth;
 95:   }
 96:   ctx->count++;
 97:   return 0;
 98: }

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

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

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

125: /*
126:    MatMFFDSetFromOptions_DS - Looks in the options database for
127:    any options appropriate for this method.

129:    Input Parameter:
130: .  ctx - the matrix free context

132: */
133: static PetscErrorCode MatMFFDSetFromOptions_DS(MatMFFD ctx, PetscOptionItems *PetscOptionsObject)
134: {
135:   MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;

137:   PetscOptionsHeadBegin(PetscOptionsObject, "Finite difference matrix free parameters");
138:   PetscOptionsReal("-mat_mffd_umin", "umin", "MatMFFDDSSetUmin", hctx->umin, &hctx->umin, NULL);
139:   PetscOptionsHeadEnd();
140:   return 0;
141: }

143: /*
144:    MatMFFDDestroy_DS - Frees the space allocated by
145:    MatCreateMFFD_DS().

147:    Input Parameter:
148: .  ctx - the matrix free context

150:    Note:
151:    Does not free the ctx, that is handled by the calling routine
152: */
153: static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx)
154: {
155:   PetscFree(ctx->hctx);
156:   return 0;
157: }

159: /*
160:    The following two routines use the PetscObjectCompose() and PetscObjectQuery()
161:    mechanism to allow the user to change the Umin parameter used in this method.
162: */
163: PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat, PetscReal umin)
164: {
165:   MatMFFD     ctx = NULL;
166:   MatMFFD_DS *hctx;

168:   MatShellGetContext(mat, &ctx);
170:   hctx       = (MatMFFD_DS *)ctx->hctx;
171:   hctx->umin = umin;
172:   return 0;
173: }

175: /*@
176:     MatMFFDDSSetUmin - Sets the "umin" parameter used by the
177:     PETSc routine for computing the differencing parameter, h, which is used
178:     for matrix-free Jacobian-vector products for a `MATMFFD` matrix.

180:    Input Parameters:
181: +  A - the `MATMFFD` matrix
182: -  umin - the parameter

184:    Level: advanced

186:    Note:
187:    See the manual page for `MatCreateSNESMF()` for a complete description of the
188:    algorithm used to compute h.

190: .seealso: `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
191: @*/
192: PetscErrorCode MatMFFDDSSetUmin(Mat A, PetscReal umin)
193: {
195:   PetscTryMethod(A, "MatMFFDDSSetUmin_C", (Mat, PetscReal), (A, umin));
196:   return 0;
197: }

199: /*MC
200:      MATMFFD_DS - algorithm for compute the "h" used in the finite difference matrix-free matrix vector product, `MatMult()`.

202:    Options Database Keys:
203: .  -mat_mffd_umin <umin> - see `MatMFFDDSSetUmin()`

205:    Level: intermediate

207:    Notes:
208:     Requires 2 norms and 1 inner product, but they are computed together
209:        so only one parallel collective operation is needed. See `MATMFFD_WP` for a method
210:        (with `KSPGMRES`) that requires NO collective operations.

212:    Formula used:
213:      F'(u)*a = [F(u+h*a) - F(u)]/h where
214:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
215:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
216:  where
217:      error_rel = square root of relative error in function evaluation
218:      umin = minimum iterate parameter

220:   References:
221: .  * -  Dennis and Schnabel, "Numerical Methods for Unconstrained Optimization and Nonlinear Equations"

223: .seealso: `MATMFFD`, `MATMFFD_WP`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_WP`, `MatMFFDDSSetUmin()`
224: M*/
225: PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx)
226: {
227:   MatMFFD_DS *hctx;

229:   /* allocate my own private data structure */
230:   PetscNew(&hctx);
231:   ctx->hctx = (void *)hctx;
232:   /* set a default for my parameter */
233:   hctx->umin = 1.e-6;

235:   /* set the functions I am providing */
236:   ctx->ops->compute        = MatMFFDCompute_DS;
237:   ctx->ops->destroy        = MatMFFDDestroy_DS;
238:   ctx->ops->view           = MatMFFDView_DS;
239:   ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS;

241:   PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDDSSetUmin_C", MatMFFDDSSetUmin_DS);
242:   return 0;
243: }