Actual source code: snesmfjdef.c

  1: #define PETSCSNES_DLL

  3: /*
  4:   Implements the DS 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 include/private/matimpl.h
 36:  #include src/snes/mf/snesmfj.h

 38: /*
 39:       The  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: } MatSNESMF_DS;

 52: /*
 53:    MatSNESMFCompute_DS - 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_DS(MatSNESMFCtx ctx,Vec U,Vec a,PetscScalar *h,PetscTruth *zeroa)
 67: {
 68:   MatSNESMF_DS     *hctx = (MatSNESMF_DS*)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_DS - 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_DS(MatSNESMFCtx ctx,PetscViewer viewer)
125: {
126:   MatSNESMF_DS     *hctx = (MatSNESMF_DS *)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_DS - 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_DS(MatSNESMFCtx ctx)
156: {
157:   PetscErrorCode   ierr;
158:   MatSNESMF_DS     *hctx = (MatSNESMF_DS*)ctx->hctx;

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

169: /*
170:    MatSNESMFDestroy_DS - Frees the space allocated by 
171:    MatSNESMFCreate_DS(). 

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_DS(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 MatSNESMFDSSetUmin_Private(Mat mat,PetscReal umin)
196: {
197:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
198:   MatSNESMF_DS *hctx;

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

212: /*@
213:     MatSNESMFDSSetUmin - Sets the "umin" parameter used by the 
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  MatSNESMFDSSetUmin(Mat A,PetscReal umin)
231: {
232:   PetscErrorCode ierr,(*f)(Mat,PetscReal);

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

243: /*MC
244:      MATSNESMF_DS - the code for compute the "h" used in the finite difference
245:             matrix-free matrix vector product.  This code
246:         implements the strategy in Dennis and Schnabel, "Numerical Methods for Unconstrained
247:         Optimization and Nonlinear Equations".

249:    Options Database Keys:
250: .  -snes_mf_umin <umin> see MatSNESMFDSSetUmin()

252:    Level: intermediate

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

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

266: .seealso: MATMFFD, MatCreateMF(), MatCreateSNESMF(), MATSNESMF_WP, MatSNESMFDSSetUmin()

268: M*/
272: PetscErrorCode  MatSNESMFCreate_DS(MatSNESMFCtx ctx)
273: {
274:   MatSNESMF_DS     *hctx;
275:   PetscErrorCode   ierr;


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

285:   /* set the functions I am providing */
286:   ctx->ops->compute        = MatSNESMFCompute_DS;
287:   ctx->ops->destroy        = MatSNESMFDestroy_DS;
288:   ctx->ops->view           = MatSNESMFView_DS;
289:   ctx->ops->setfromoptions = MatSNESMFSetFromOptions_DS;

291:   PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFDSSetUmin_C",
292:                             "MatSNESMFDSSetUmin_Private",
293:                              MatSNESMFDSSetUmin_Private);
294:   return(0);
295: }