1: #define PETSCSNES_DLL

  3: /*MC
  4:      MATSNESMF_WP - Implements an alternative approach for computing the differencing parameter
  5:         h used with the finite difference based matrix-free Jacobian.  This code
  6:         implements the strategy of M. Pernice and H. Walker:

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

 10:       Notes:
 11:         1) || U || does not change between linear iterations so is reused
 12:         2) In GMRES || a || == 1 and so does not need to ever be computed except at restart
 13:            when it is recomputed.

 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:    Options Database Keys:
 20: .   -snes_mf_compute_normu -Compute the norm of u everytime see MatSNESMFWPSetComputeNormU()


 23:    Level: intermediate

 25:    Notes: Requires no global collectives when used with GMRES

 27:    Formula used:
 28:      F'(u)*a = [F(u+h*a) - F(u)]/h where

 30: .seealso: MATMFFD, MatCreateMF(), MatCreateSNESMF(), MATSNESMF_DS

 32: M*/

 34: /*
 35:     This include file defines the data structure  MatSNESMF that 
 36:    includes information about the computation of h. It is shared by 
 37:    all implementations that people provide.

 39:    See snesmfjdef.c for  a full set of comments on the routines below.
 40: */
 41:  #include include/private/matimpl.h
 42:  #include src/snes/mf/snesmfj.h

 44: typedef struct {
 45:   PetscReal  normUfact;                   /* previous sqrt(1.0 + || U ||) */
 46:   PetscTruth computenorma,computenormU;
 47: } MatSNESMFWP;

 51: /*
 52:      MatSNESMFCompute_WP - Standard PETSc code for 
 53:    computing h with matrix-free finite differences.

 55:   Input Parameters:
 56: +   ctx - the matrix free context
 57: .   U - the location at which you want the Jacobian
 58: -   a - the direction you want the derivative

 60:   Output Parameter:
 61: .   h - the scale computed

 63: */
 64: static PetscErrorCode MatSNESMFCompute_WP(MatSNESMFCtx ctx,Vec U,Vec a,PetscScalar *h,PetscTruth *zeroa)
 65: {
 66:   MatSNESMFWP    *hctx = (MatSNESMFWP*)ctx->hctx;
 67:   PetscReal      normU,norma = 1.0;


 72:   if (!(ctx->count % ctx->recomputeperiod)) {
 73:     if (hctx->computenorma && (hctx->computenormU || !ctx->ncurrenth)) {
 74:       VecNormBegin(U,NORM_2,&normU);
 75:       VecNormBegin(a,NORM_2,&norma);
 76:       VecNormEnd(U,NORM_2,&normU);
 77:       VecNormEnd(a,NORM_2,&norma);
 78:       hctx->normUfact = sqrt(1.0+normU);
 79:     } else if (hctx->computenormU || !ctx->ncurrenth) {
 80:       VecNorm(U,NORM_2,&normU);
 81:       hctx->normUfact = sqrt(1.0+normU);
 82:     } else if (hctx->computenorma) {
 83:       VecNorm(a,NORM_2,&norma);
 84:     }
 85:     if (norma == 0.0) {
 86:       *zeroa = PETSC_TRUE;
 87:       return(0);
 88:     }
 89:     *zeroa = PETSC_FALSE;
 90:     *h = ctx->error_rel*hctx->normUfact/norma;
 91:   } else {
 92:     *h = ctx->currenth;
 93:   }
 94:   return(0);
 95: }

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

105:   Input Parameters:
106: +   ctx - the matrix free context
107: -   viewer - the PETSc viewer

109: */
110: static PetscErrorCode MatSNESMFView_WP(MatSNESMFCtx ctx,PetscViewer viewer)
111: {
112:   MatSNESMFWP    *hctx = (MatSNESMFWP *)ctx->hctx;
114:   PetscTruth     iascii;

117:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
118:   if (iascii) {
119:     if (hctx->computenorma){PetscViewerASCIIPrintf(viewer,"    Computes normA\n");}
120:     else                   { PetscViewerASCIIPrintf(viewer,"    Does not compute normA\n");}
121:     if (hctx->computenormU){ PetscViewerASCIIPrintf(viewer,"    Computes normU\n");}
122:     else                   { PetscViewerASCIIPrintf(viewer,"    Does not compute normU\n");}
123:   } else {
124:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES matrix-free WP",((PetscObject)viewer)->type_name);
125:   }
126:   return(0);
127: }

131: /*
132:    MatSNESMFSetFromOptions_WP - Looks in the options database for 
133:      any options appropriate for this method

135:   Input Parameter:
136: .  ctx - the matrix free context

138: */
139: static PetscErrorCode MatSNESMFSetFromOptions_WP(MatSNESMFCtx ctx)
140: {
142:   MatSNESMFWP    *hctx = (MatSNESMFWP*)ctx->hctx;

145:   PetscOptionsHead("Walker-Pernice options");
146:     PetscOptionsTruth("-snes_mf_compute_normu","Compute the norm of u","MatSNESMFWPSetComputeNormU",
147:                           hctx->computenorma,&hctx->computenorma,0);
148:   PetscOptionsTail();
149:   return(0);
150: }

154: /*
155:    MatSNESMFDestroy_WP - Frees the space allocated by 
156:        MatSNESMFCreate_WP(). 

158:   Input Parameter:
159: .  ctx - the matrix free context

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

163: */
164: static PetscErrorCode MatSNESMFDestroy_WP(MatSNESMFCtx ctx)
165: {
168:   PetscFree(ctx->hctx);
169:   return(0);
170: }

175: PetscErrorCode  MatSNESMFWPSetComputeNormU_P(Mat mat,PetscTruth flag)
176: {
177:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
178:   MatSNESMFWP  *hctx;

181:   hctx               = (MatSNESMFWP*)ctx->hctx;
182:   hctx->computenormU = flag;
183:   return(0);
184: }

189: /*@
190:     MatSNESMFWPSetComputeNormU - Sets whether it computes the ||U|| used by the WP
191:              PETSc routine for computing h. With any Krylov solver this need only 
192:              be computed during the first iteration and kept for later.

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

198:   Options Database Key:
199: .   -snes_mf_compute_normu <true,false> - true by default, false can save calculations but you 
200:               must be sure that ||U|| has not changed in the mean time.

202:   Level: advanced

204:   Notes:
205:    See the manual page for MATSNESMF_WP for a complete description of the
206:    algorithm used to compute h.

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

210: @*/
211: PetscErrorCode  MatSNESMFWPSetComputeNormU(Mat A,PetscTruth flag)
212: {
213:   PetscErrorCode ierr,(*f)(Mat,PetscTruth);

217:   PetscObjectQueryFunction((PetscObject)A,"MatSNESMFWPSetComputeNormU_C",(void (**)(void))&f);
218:   if (f) {
219:     (*f)(A,flag);
220:   }
221:   return(0);
222: }

227: /*
228:      MatSNESMFCreate_WP - Standard PETSc code for 
229:    computing h with matrix-free finite differences.

231:    Input Parameter:
232: .  ctx - the matrix free context created by MatSNESMFCreate()

234: */
235: PetscErrorCode  MatSNESMFCreate_WP(MatSNESMFCtx ctx)
236: {
238:   MatSNESMFWP    *hctx;


242:   /* allocate my own private data structure */
243:   PetscNew(MatSNESMFWP,&hctx);
244:   ctx->hctx          = (void*)hctx;
245:   hctx->computenormU = PETSC_FALSE;
246:   hctx->computenorma = PETSC_TRUE;

248:   /* set the functions I am providing */
249:   ctx->ops->compute        = MatSNESMFCompute_WP;
250:   ctx->ops->destroy        = MatSNESMFDestroy_WP;
251:   ctx->ops->view           = MatSNESMFView_WP;
252:   ctx->ops->setfromoptions = MatSNESMFSetFromOptions_WP;

254:   PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatSNESMFWPSetComputeNormU_C",
255:                             "MatSNESMFWPSetComputeNormU_P",
256:                              MatSNESMFWPSetComputeNormU_P);

258:   return(0);
259: }