Actual source code: wp.c

  1: /*$Id: wp.c,v 1.36 2001/08/06 21:17:13 bsmith Exp $*/
  2: /*
  3:   Implements an alternative approach for computing the differencing parameter
  4:   h used with the finite difference based matrix-free Jacobian.  This code
  5:   implements the strategy of M. Pernice and H. Walker:

  7:       h = error_rel * sqrt(1 + ||U||) / ||a||

  9:       Notes:
 10:         1) || U || does not change between linear iterations so can be reused
 11:         2) In GMRES || a || == 1 and so does not need to ever be computed if you never 
 12:            have a restart. Unfortunately a RESTART computes a matrix vector product 
 13:            with ||a|| != 0 which breaks this

 15:       Reference:  M. Pernice and H. F. Walker, "NITSOL: A Newton Iterative 
 16:       Solver for Nonlinear Systems", SIAM J. Sci. Stat. Comput.", 1998, 
 17:       vol 19, pp. 302--318.

 19:    See snesmfjdef.c for  a full set of comments on the routines below.
 20: */

 22: /*
 23:     This include file defines the data structure  MatSNESMF that 
 24:    includes information about the computation of h. It is shared by 
 25:    all implementations that people provide.
 26: */
 27:  #include src/mat/matimpl.h
 28:  #include src/snes/mf/snesmfj.h

 30: typedef struct {
 31:   PetscReal  normUfact;                   /* previous sqrt(1.0 + || U ||) */
 32:   PetscTruth computenorma,computenormU;
 33: } MatSNESMFWP;

 35: #undef __FUNCT__  
 37: /*
 38:      MatSNESMFCompute_WP - Standard PETSc code for 
 39:    computing h with matrix-free finite differences.

 41:   Input Parameters:
 42: +   ctx - the matrix free context
 43: .   U - the location at which you want the Jacobian
 44: -   a - the direction you want the derivative

 46:   Output Parameter:
 47: .   h - the scale computed

 49: */
 50: static int MatSNESMFCompute_WP(MatSNESMFCtx ctx,Vec U,Vec a,PetscScalar *h)
 51: {
 52:   MatSNESMFWP   *hctx = (MatSNESMFWP*)ctx->hctx;
 53:   PetscReal     normU,norma = 1.0;
 54:   int           ierr;


 58:   if (!(ctx->count % ctx->recomputeperiod)) {
 59:     if (hctx->computenorma && (hctx->computenormU || !ctx->ncurrenth)) {
 60:       VecNormBegin(U,NORM_2,&normU);
 61:       VecNormBegin(a,NORM_2,&norma);
 62:       VecNormEnd(U,NORM_2,&normU);
 63:       VecNormEnd(a,NORM_2,&norma);
 64:       hctx->normUfact = sqrt(1.0+normU);
 65:     } else if (hctx->computenormU || !ctx->ncurrenth) {
 66:       VecNorm(U,NORM_2,&normU);
 67:       hctx->normUfact = sqrt(1.0+normU);
 68:     } else if (hctx->computenorma) {
 69:       VecNorm(a,NORM_2,&norma);
 70:     }
 71:     *h = ctx->error_rel*hctx->normUfact/norma;
 72:   } else {
 73:     *h = ctx->currenth;
 74:   }
 75:   return(0);
 76: }

 78: #undef __FUNCT__  
 80: /*
 81:    MatSNESMFView_WP - Prints information about this particular 
 82:      method for computing h. Note that this does not print the general
 83:      information about the matrix free, that is printed by the calling
 84:      routine.

 86:   Input Parameters:
 87: +   ctx - the matrix free context
 88: -   viewer - the PETSc viewer

 90: */
 91: static int MatSNESMFView_WP(MatSNESMFCtx ctx,PetscViewer viewer)
 92: {
 93:   MatSNESMFWP *hctx = (MatSNESMFWP *)ctx->hctx;
 94:   int         ierr;
 95:   PetscTruth  isascii;

 98:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
 99:   if (isascii) {
100:     if (hctx->computenorma){PetscViewerASCIIPrintf(viewer,"    Computes normAn");}
101:     else                   { PetscViewerASCIIPrintf(viewer,"    Does not compute normAn");}
102:     if (hctx->computenormU){ PetscViewerASCIIPrintf(viewer,"    Computes normUn");}
103:     else                   { PetscViewerASCIIPrintf(viewer,"    Does not compute normUn");}
104:   } else {
105:     SETERRQ1(1,"Viewer type %s not supported for SNES matrix-free WP",((PetscObject)viewer)->type_name);
106:   }
107:   return(0);
108: }

110: #undef __FUNCT__  
112: /*
113:    MatSNESMFSetFromOptions_WP - Looks in the options database for 
114:      any options appropriate for this method

116:   Input Parameter:
117: .  ctx - the matrix free context

119: */
120: static int MatSNESMFSetFromOptions_WP(MatSNESMFCtx ctx)
121: {
122:   int         ierr;
123:   MatSNESMFWP *hctx = (MatSNESMFWP*)ctx->hctx;

126:   PetscOptionsHead("Walker-Pernice options");
127:     PetscOptionsLogical("-snes_mf_compute_norma","Compute the norm of a","MatSNESMFWPSetComputeNormA",
128:                           hctx->computenorma,&hctx->computenorma,0);
129:     PetscOptionsLogical("-snes_mf_compute_normu","Compute the norm of u","MatSNESMFWPSetComputeNormU",
130:                           hctx->computenorma,&hctx->computenorma,0);
131:   PetscOptionsTail();
132:   return(0);
133: }

135: #undef __FUNCT__  
137: /*
138:    MatSNESMFDestroy_WP - Frees the space allocated by 
139:        MatSNESMFCreate_WP(). 

141:   Input Parameter:
142: .  ctx - the matrix free context

144:    Notes: does not free the ctx, that is handled by the calling routine

146: */
147: static int MatSNESMFDestroy_WP(MatSNESMFCtx ctx)
148: {
151:   PetscFree(ctx->hctx);
152:   return(0);
153: }

155: EXTERN_C_BEGIN
156: #undef __FUNCT__  
158: int MatSNESMFWPSetComputeNormA_P(Mat mat,PetscTruth flag)
159: {
160:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
161:   MatSNESMFWP  *hctx;

164:   hctx               = (MatSNESMFWP*)ctx->hctx;
165:   hctx->computenorma = flag;
166:   return(0);
167: }
168: EXTERN_C_END

170: #undef __FUNCT__  
172: /*@
173:     MatSNESMFWPSetComputeNormA - Sets whether it computes the ||a|| used by the WP
174:              PETSc routine for computing h. With GMRES since the ||a|| is always
175:              one, you can save communication by setting this to false.

177:   Input Parameters:
178: +   A - the matrix created with MatCreateSNESMF()
179: -   flag - PETSC_TRUE causes it to compute ||a||, PETSC_FALSE assumes it is 1.

181:   Level: advanced

183:   Notes:
184:    See the manual page for MatCreateSNESMF() for a complete description of the
185:    algorithm used to compute h.

187: .seealso: MatSNESMFSetFunctionError(), MatCreateSNESMF()

189: @*/
190: int MatSNESMFWPSetComputeNormA(Mat A,PetscTruth flag)
191: {
192:   int ierr,(*f)(Mat,PetscTruth);

196:   PetscObjectQueryFunction((PetscObject)A,"MatSNESMFWPSetComputeNormA_C",(void (**)(void))&f);
197:   if (f) {
198:     (*f)(A,flag);
199:   }
200:   return(0);
201: }

203: EXTERN_C_BEGIN
204: #undef __FUNCT__  
206: int MatSNESMFWPSetComputeNormU_P(Mat mat,PetscTruth flag)
207: {
208:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
209:   MatSNESMFWP  *hctx;

212:   hctx               = (MatSNESMFWP*)ctx->hctx;
213:   hctx->computenormU = flag;
214:   return(0);
215: }
216: EXTERN_C_END

218: #undef __FUNCT__  
220: /*@
221:     MatSNESMFWPSetComputeNormU - Sets whether it computes the ||U|| used by the WP
222:              PETSc routine for computing h. With any Krylov solver this need only 
223:              be computed during the first iteration and kept for later.

225:   Input Parameters:
226: +   A - the matrix created with MatCreateSNESMF()
227: -   flag - PETSC_TRUE causes it to compute ||U||, PETSC_FALSE uses the previous value

229:   Level: advanced

231:   Notes:
232:    See the manual page for MatCreateSNESMF() for a complete description of the
233:    algorithm used to compute h.

235: .seealso: MatSNESMFSetFunctionError(), MatCreateSNESMF()

237: @*/
238: int MatSNESMFWPSetComputeNormU(Mat A,PetscTruth flag)
239: {
240:   int ierr,(*f)(Mat,PetscTruth);

244:   PetscObjectQueryFunction((PetscObject)A,"MatSNESMFWPSetComputeNormU_C",(void (**)(void))&f);
245:   if (f) {
246:     (*f)(A,flag);
247:   }
248:   return(0);
249: }

251: EXTERN_C_BEGIN
252: #undef __FUNCT__  
254: /*
255:      MatSNESMFCreate_WP - Standard PETSc code for 
256:    computing h with matrix-free finite differences.

258:    Input Parameter:
259: .  ctx - the matrix free context created by MatSNESMFCreate()

261: */
262: int MatSNESMFCreate_WP(MatSNESMFCtx ctx)
263: {
264:   int         ierr;
265:   MatSNESMFWP *hctx;


269:   /* allocate my own private data structure */
270:   ierr               = PetscNew(MatSNESMFWP,&hctx);
271:   ctx->hctx          = (void*)hctx;
272:   hctx->computenormU = PETSC_FALSE;
273:   hctx->computenorma = PETSC_TRUE;

275:   /* set the functions I am providing */
276:   ctx->ops->compute        = MatSNESMFCompute_WP;
277:   ctx->ops->destroy        = MatSNESMFDestroy_WP;
278:   ctx->ops->view           = MatSNESMFView_WP;
279:   ctx->ops->setfromoptions = MatSNESMFSetFromOptions_WP;

281:   PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFWPSetComputeNormA_C",
282:                             "MatSNESMFWPSetComputeNormA_P",
283:                              MatSNESMFWPSetComputeNormA_P);
284:   PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFWPSetComputeNormU_C",
285:                             "MatSNESMFWPSetComputeNormU_P",
286:                              MatSNESMFWPSetComputeNormU_P);

288:   return(0);
289: }
290: EXTERN_C_END