Actual source code: mffddef.c

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

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

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

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

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

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

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

 25: */

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

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

 47: static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx, Vec U, Vec a, PetscScalar *h, PetscBool *zeroa)
 48: {
 49:   MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
 50:   PetscReal   nrm, sum, umin = hctx->umin;
 51:   PetscScalar dot;

 53:   PetscFunctionBegin;
 54:   if (!(ctx->count % ctx->recomputeperiod)) {
 55:     /*
 56:      This algorithm requires 2 norms and 1 inner product. Rather than
 57:      use directly the VecNorm() and VecDot() routines (and thus have
 58:      three separate collective operations, we use the VecxxxBegin/End() routines
 59:     */
 60:     PetscCall(VecDotBegin(U, a, &dot));
 61:     PetscCall(VecNormBegin(a, NORM_1, &sum));
 62:     PetscCall(VecNormBegin(a, NORM_2, &nrm));
 63:     PetscCall(VecDotEnd(U, a, &dot));
 64:     PetscCall(VecNormEnd(a, NORM_1, &sum));
 65:     PetscCall(VecNormEnd(a, NORM_2, &nrm));

 67:     if (nrm == 0.0) {
 68:       *zeroa = PETSC_TRUE;
 69:       PetscFunctionReturn(PETSC_SUCCESS);
 70:     }
 71:     *zeroa = PETSC_FALSE;

 73:     /*
 74:       Safeguard for step sizes that are "too small"
 75:     */
 76:     if (PetscAbsScalar(dot) < umin * sum && PetscRealPart(dot) >= 0.0) dot = umin * sum;
 77:     else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin * sum) dot = -umin * sum;
 78:     *h = ctx->error_rel * dot / (nrm * nrm);
 79:     PetscCheck(!PetscIsInfOrNanScalar(*h), 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);
 80:   } else {
 81:     *h = ctx->currenth;
 82:   }
 83:   ctx->count++;
 84:   PetscFunctionReturn(PETSC_SUCCESS);
 85: }

 87: static PetscErrorCode MatMFFDView_DS(MatMFFD ctx, PetscViewer viewer)
 88: {
 89:   MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
 90:   PetscBool   iascii;

 92:   PetscFunctionBegin;
 93:   /*
 94:      Currently this only handles the ascii file viewers, others
 95:      could be added, but for this type of object other viewers
 96:      make less sense
 97:   */
 98:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 99:   if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, "    umin=%g (minimum iterate parameter)\n", (double)hctx->umin));
100:   PetscFunctionReturn(PETSC_SUCCESS);
101: }

103: static PetscErrorCode MatMFFDSetFromOptions_DS(MatMFFD ctx, PetscOptionItems *PetscOptionsObject)
104: {
105:   MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;

107:   PetscFunctionBegin;
108:   PetscOptionsHeadBegin(PetscOptionsObject, "Finite difference matrix-free parameters");
109:   PetscCall(PetscOptionsReal("-mat_mffd_umin", "umin", "MatMFFDDSSetUmin", hctx->umin, &hctx->umin, NULL));
110:   PetscOptionsHeadEnd();
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

114: static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx)
115: {
116:   PetscFunctionBegin;
117:   PetscCall(PetscFree(ctx->hctx));
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

121: /*
122:    The following two routines use the PetscObjectCompose() and PetscObjectQuery()
123:    mechanism to allow the user to change the Umin parameter used in this method.
124: */
125: static PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat, PetscReal umin)
126: {
127:   MatMFFD     ctx = NULL;
128:   MatMFFD_DS *hctx;

130:   PetscFunctionBegin;
131:   PetscCall(MatShellGetContext(mat, &ctx));
132:   PetscCheck(ctx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatMFFDDSSetUmin() attached to non-shell matrix");
133:   hctx       = (MatMFFD_DS *)ctx->hctx;
134:   hctx->umin = umin;
135:   PetscFunctionReturn(PETSC_SUCCESS);
136: }

138: /*@
139:   MatMFFDDSSetUmin - Sets the "umin" parameter used by the
140:   PETSc routine for computing the differencing parameter, h, which is used
141:   for matrix-free Jacobian-vector products for a `MATMFFD` matrix.

143:   Input Parameters:
144: + A    - the `MATMFFD` matrix
145: - umin - the parameter

147:   Level: advanced

149:   Note:
150:   See the manual page for `MatCreateSNESMF()` for a complete description of the
151:   algorithm used to compute h.

153: .seealso: `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
154: @*/
155: PetscErrorCode MatMFFDDSSetUmin(Mat A, PetscReal umin)
156: {
157:   PetscFunctionBegin;
159:   PetscTryMethod(A, "MatMFFDDSSetUmin_C", (Mat, PetscReal), (A, umin));
160:   PetscFunctionReturn(PETSC_SUCCESS);
161: }

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

166:    Options Database Keys:
167: .  -mat_mffd_umin <umin> - see `MatMFFDDSSetUmin()`

169:    Level: intermediate

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

176:    Formula used:
177:      F'(u)*a = [F(u+h*a) - F(u)]/h where
178:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
179:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
180:   where
181:      error_rel = square root of relative error in function evaluation
182:      umin = minimum iterate parameter

184:   Method taken from {cite}`dennis:83`

186: .seealso: `MATMFFD`, `MATMFFD_WP`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_WP`, `MatMFFDDSSetUmin()`
187: M*/
188: PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx)
189: {
190:   MatMFFD_DS *hctx;

192:   PetscFunctionBegin;
193:   /* allocate my own private data structure */
194:   PetscCall(PetscNew(&hctx));
195:   ctx->hctx = (void *)hctx;
196:   /* set a default for my parameter */
197:   hctx->umin = 1.e-6;

199:   /* set the functions I am providing */
200:   ctx->ops->compute        = MatMFFDCompute_DS;
201:   ctx->ops->destroy        = MatMFFDDestroy_DS;
202:   ctx->ops->view           = MatMFFDView_DS;
203:   ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS;

205:   PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDDSSetUmin_C", MatMFFDDSSetUmin_DS));
206:   PetscFunctionReturn(PETSC_SUCCESS);
207: }