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: static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx, Vec U, Vec a, PetscScalar *h, PetscBool *zeroa)
 49: {
 50:   MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx;
 51:   PetscReal   nrm, sum, umin = hctx->umin;
 52:   PetscScalar dot;

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

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

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

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

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

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

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

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

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

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

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

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

148:    Level: advanced

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

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

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

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

170:    Level: intermediate

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

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

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

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

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

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

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