Actual source code: snesmfj2.c

petsc-3.9.0 2018-04-07
Report Typos and Errors

  2:  #include <petsc/private/snesimpl.h>
  3: /* matimpl.h is needed only for logging of matrix operation */
  4:  #include <petsc/private/matimpl.h>

  6: extern PetscErrorCode SNESDiffParameterCreate_More(SNES,Vec,void**);
  7: extern PetscErrorCode SNESDiffParameterCompute_More(SNES,void*,Vec,Vec,PetscReal*,PetscReal*);
  8: extern PetscErrorCode SNESDiffParameterDestroy_More(void*);

 10: typedef struct {  /* default context for matrix-free SNES */
 11:   SNES         snes;             /* SNES context */
 12:   Vec          w;                /* work vector */
 13:   MatNullSpace sp;               /* null space context */
 14:   PetscReal    error_rel;        /* square root of relative error in computing function */
 15:   PetscReal    umin;             /* minimum allowable u'a value relative to |u|_1 */
 16:   PetscBool    jorge;            /* flag indicating use of Jorge's method for determining the differencing parameter */
 17:   PetscReal    h;                /* differencing parameter */
 18:   PetscBool    need_h;           /* flag indicating whether we must compute h */
 19:   PetscBool    need_err;         /* flag indicating whether we must currently compute error_rel */
 20:   PetscBool    compute_err;      /* flag indicating whether we must ever compute error_rel */
 21:   PetscInt     compute_err_iter; /* last iter where we've computer error_rel */
 22:   PetscInt     compute_err_freq; /* frequency of computing error_rel */
 23:   void         *data;            /* implementation-specific data */
 24: } MFCtx_Private;

 26: PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat)
 27: {
 29:   MFCtx_Private  *ctx;

 32:   MatShellGetContext(mat,(void**)&ctx);
 33:   VecDestroy(&ctx->w);
 34:   MatNullSpaceDestroy(&ctx->sp);
 35:   if (ctx->jorge || ctx->compute_err) {SNESDiffParameterDestroy_More(ctx->data);}
 36:   PetscFree(ctx);
 37:   return(0);
 38: }

 40: /*
 41:    SNESMatrixFreeView2_Private - Views matrix-free parameters.
 42:  */
 43: PetscErrorCode SNESMatrixFreeView2_Private(Mat J,PetscViewer viewer)
 44: {
 46:   MFCtx_Private  *ctx;
 47:   PetscBool      iascii;

 50:   MatShellGetContext(J,(void**)&ctx);
 51:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 52:   if (iascii) {
 53:     PetscViewerASCIIPrintf(viewer,"  SNES matrix-free approximation:\n");
 54:     if (ctx->jorge) {
 55:       PetscViewerASCIIPrintf(viewer,"    using Jorge's method of determining differencing parameter\n");
 56:     }
 57:     PetscViewerASCIIPrintf(viewer,"    err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
 58:     PetscViewerASCIIPrintf(viewer,"    umin=%g (minimum iterate parameter)\n",(double)ctx->umin);
 59:     if (ctx->compute_err) {
 60:       PetscViewerASCIIPrintf(viewer,"    freq_err=%D (frequency for computing err)\n",ctx->compute_err_freq);
 61:     }
 62:   }
 63:   return(0);
 64: }

 66: /*
 67:   SNESMatrixFreeMult2_Private - Default matrix-free form for Jacobian-vector
 68:   product, y = F'(u)*a:
 69:         y = (F(u + ha) - F(u)) /h,
 70:   where F = nonlinear function, as set by SNESSetFunction()
 71:         u = current iterate
 72:         h = difference interval
 73: */
 74: PetscErrorCode SNESMatrixFreeMult2_Private(Mat mat,Vec a,Vec y)
 75: {
 76:   MFCtx_Private  *ctx;
 77:   SNES           snes;
 78:   PetscReal      h,norm,sum,umin,noise;
 79:   PetscScalar    hs,dot;
 80:   Vec            w,U,F;
 81:   PetscErrorCode ierr,(*eval_fct)(SNES,Vec,Vec);
 82:   MPI_Comm       comm;
 83:   PetscInt       iter;

 86:   /* We log matrix-free matrix-vector products separately, so that we can
 87:      separate the performance monitoring from the cases that use conventional
 88:      storage.  We may eventually modify event logging to associate events
 89:      with particular objects, hence alleviating the more general problem. */
 90:   PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);

 92:   PetscObjectGetComm((PetscObject)mat,&comm);
 93:   MatShellGetContext(mat,(void**)&ctx);
 94:   snes = ctx->snes;
 95:   w    = ctx->w;
 96:   umin = ctx->umin;

 98:   SNESGetSolution(snes,&U);
 99:   eval_fct = SNESComputeFunction;
100:   SNESGetFunction(snes,&F,NULL,NULL);

102:   /* Determine a "good" step size, h */
103:   if (ctx->need_h) {

105:     /* Use Jorge's method to compute h */
106:     if (ctx->jorge) {
107:       SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);

109:       /* Use the Brown/Saad method to compute h */
110:     } else {
111:       /* Compute error if desired */
112:       SNESGetIterationNumber(snes,&iter);
113:       if ((ctx->need_err) || ((ctx->compute_err_freq) && (ctx->compute_err_iter != iter) && (!((iter-1)%ctx->compute_err_freq)))) {
114:         /* Use Jorge's method to compute noise */
115:         SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);

117:         ctx->error_rel = PetscSqrtReal(noise);

119:         PetscInfo3(snes,"Using Jorge's noise: noise=%g, sqrt(noise)=%g, h_more=%g\n",(double)noise,(double)ctx->error_rel,(double)h);

121:         ctx->compute_err_iter = iter;
122:         ctx->need_err         = PETSC_FALSE;
123:       }

125:       VecDotBegin(U,a,&dot);
126:       VecNormBegin(a,NORM_1,&sum);
127:       VecNormBegin(a,NORM_2,&norm);
128:       VecDotEnd(U,a,&dot);
129:       VecNormEnd(a,NORM_1,&sum);
130:       VecNormEnd(a,NORM_2,&norm);


133:       /* Safeguard for step sizes too small */
134:       if (sum == 0.0) {
135:         dot  = 1.0;
136:         norm = 1.0;
137:       } else if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
138:       else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
139:       h = PetscRealPart(ctx->error_rel*dot/(norm*norm));
140:     }
141:   } else h = ctx->h;

143:   if (!ctx->jorge || !ctx->need_h) {PetscInfo1(snes,"h = %g\n",(double)h);}

145:   /* Evaluate function at F(u + ha) */
146:   hs   = h;
147:   VecWAXPY(w,hs,a,U);
148:   eval_fct(snes,w,y);
149:   VecAXPY(y,-1.0,F);
150:   VecScale(y,1.0/hs);
151:   if (mat->nullsp) {MatNullSpaceRemove(mat->nullsp,y);}

153:   PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
154:   return(0);
155: }

157: /*@C
158:    SNESMatrixFreeCreate2 - Creates a matrix-free matrix
159:    context for use with a SNES solver.  This matrix can be used as
160:    the Jacobian argument for the routine SNESSetJacobian().

162:    Input Parameters:
163: .  snes - the SNES context
164: .  x - vector where SNES solution is to be stored.

166:    Output Parameter:
167: .  J - the matrix-free matrix

169:    Level: advanced

171:    Notes:
172:    The matrix-free matrix context merely contains the function pointers
173:    and work space for performing finite difference approximations of
174:    Jacobian-vector products, J(u)*a, via

176: $       J(u)*a = [J(u+h*a) - J(u)]/h,
177: $   where by default
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
182: $                    function evaluation
183: $        umin = minimum iterate parameter
184: $   Alternatively, the differencing parameter, h, can be set using
185: $   Jorge's nifty new strategy if one specifies the option
186: $          -snes_mf_jorge

188:    The user can set these parameters via MatMFFDSetFunctionError().
189:    See Users-Manual: ch_snes for details

191:    The user should call MatDestroy() when finished with the matrix-free
192:    matrix context.

194:    Options Database Keys:
195: $  -snes_mf_err <error_rel>
196: $  -snes_mf_unim <umin>
197: $  -snes_mf_compute_err
198: $  -snes_mf_freq_err <freq>
199: $  -snes_mf_jorge

201: .keywords: SNES, default, matrix-free, create, matrix

203: .seealso: MatDestroy(), MatMFFDSetFunctionError()
204: @*/
205: PetscErrorCode  SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J)
206: {
207:   MPI_Comm       comm;
208:   MFCtx_Private  *mfctx;
210:   PetscInt       n,nloc;
211:   PetscBool      flg;
212:   char           p[64];

215:   PetscNewLog(snes,&mfctx);
216:   mfctx->sp               = 0;
217:   mfctx->snes             = snes;
218:   mfctx->error_rel        = PETSC_SQRT_MACHINE_EPSILON;
219:   mfctx->umin             = 1.e-6;
220:   mfctx->h                = 0.0;
221:   mfctx->need_h           = PETSC_TRUE;
222:   mfctx->need_err         = PETSC_FALSE;
223:   mfctx->compute_err      = PETSC_FALSE;
224:   mfctx->compute_err_freq = 0;
225:   mfctx->compute_err_iter = -1;
226:   mfctx->compute_err      = PETSC_FALSE;
227:   mfctx->jorge            = PETSC_FALSE;

229:   PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_err",&mfctx->error_rel,NULL);
230:   PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_umin",&mfctx->umin,NULL);
231:   PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_jorge",&mfctx->jorge,NULL);
232:   PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_compute_err",&mfctx->compute_err,NULL);
233:   PetscOptionsGetInt(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_freq_err",&mfctx->compute_err_freq,&flg);
234:   if (flg) {
235:     if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0;
236:     mfctx->compute_err = PETSC_TRUE;
237:   }
238:   if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE;
239:   if (mfctx->jorge || mfctx->compute_err) {
240:     SNESDiffParameterCreate_More(snes,x,&mfctx->data);
241:   } else mfctx->data = 0;

243:   PetscOptionsHasName(((PetscObject)snes)->options,NULL,"-help",&flg);
244:   PetscStrncpy(p,"-",sizeof(p));
245:   if (((PetscObject)snes)->prefix) {PetscStrlcat(p,((PetscObject)snes)->prefix,sizeof(p));}
246:   if (flg) {
247:     PetscPrintf(PetscObjectComm((PetscObject)snes)," Matrix-free Options (via SNES):\n");
248:     PetscPrintf(PetscObjectComm((PetscObject)snes),"   %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n",p,(double)mfctx->error_rel);
249:     PetscPrintf(PetscObjectComm((PetscObject)snes),"   %ssnes_mf_umin <umin>: see users manual (default %g)\n",p,(double)mfctx->umin);
250:     PetscPrintf(PetscObjectComm((PetscObject)snes),"   %ssnes_mf_jorge: use Jorge More's method\n",p);
251:     PetscPrintf(PetscObjectComm((PetscObject)snes),"   %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);
252:     PetscPrintf(PetscObjectComm((PetscObject)snes),"   %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);
253:     PetscPrintf(PetscObjectComm((PetscObject)snes),"   %ssnes_mf_noise_file <file>: set file for printing noise info\n",p);
254:   }
255:   VecDuplicate(x,&mfctx->w);
256:   PetscObjectGetComm((PetscObject)x,&comm);
257:   VecGetSize(x,&n);
258:   VecGetLocalSize(x,&nloc);
259:   MatCreate(comm,J);
260:   MatSetSizes(*J,nloc,n,n,n);
261:   MatSetType(*J,MATSHELL);
262:   MatShellSetContext(*J,mfctx);
263:   MatShellSetOperation(*J,MATOP_MULT,(void (*)(void))SNESMatrixFreeMult2_Private);
264:   MatShellSetOperation(*J,MATOP_DESTROY,(void (*)(void))SNESMatrixFreeDestroy2_Private);
265:   MatShellSetOperation(*J,MATOP_VIEW,(void (*)(void))SNESMatrixFreeView2_Private);
266:   MatSetUp(*J);

268:   PetscLogObjectParent((PetscObject)*J,(PetscObject)mfctx->w);
269:   PetscLogObjectParent((PetscObject)snes,(PetscObject)*J);
270:   return(0);
271: }

273: /*@C
274:    SNESDefaultMatrixFreeSetParameters2 - Sets the parameters for the approximation of
275:    matrix-vector products using finite differences.

277: $       J(u)*a = [J(u+h*a) - J(u)]/h where

279:    either the user sets h directly here, or this parameter is computed via

281: $        h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
282: $          = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
283: $

285:    Input Parameters:
286: +  mat - the matrix
287: .  error_rel - relative error (should be set to the square root of
288:      the relative error in the function evaluations)
289: .  umin - minimum allowable u-value
290: -  h - differencing parameter

292:    Level: advanced

294:    Notes:
295:    If the user sets the parameter h directly, then this value will be used
296:    instead of the default computation indicated above.

298: .keywords: SNES, matrix-free, parameters

300: .seealso: MatCreateSNESMF()
301: @*/
302: PetscErrorCode  SNESDefaultMatrixFreeSetParameters2(Mat mat,PetscReal error,PetscReal umin,PetscReal h)
303: {
304:   MFCtx_Private  *ctx;

308:   MatShellGetContext(mat,(void**)&ctx);
309:   if (ctx) {
310:     if (error != PETSC_DEFAULT) ctx->error_rel = error;
311:     if (umin  != PETSC_DEFAULT) ctx->umin = umin;
312:     if (h     != PETSC_DEFAULT) {
313:       ctx->h      = h;
314:       ctx->need_h = PETSC_FALSE;
315:     }
316:   }
317:   return(0);
318: }

320: PetscErrorCode  SNESUnSetMatrixFreeParameter(SNES snes)
321: {
322:   MFCtx_Private  *ctx;
324:   Mat            mat;

327:   SNESGetJacobian(snes,&mat,NULL,NULL,NULL);
328:   MatShellGetContext(mat,(void**)&ctx);
329:   if (ctx) ctx->need_h = PETSC_TRUE;
330:   return(0);
331: }