Actual source code: snesmfj2.c

      1: #define PETSCSNES_DLL

     3:  #include src/snes/snesimpl.h

      5: EXTERN PetscErrorCode DiffParameterCreate_More(SNES,Vec,void**);
      6: EXTERN PetscErrorCode DiffParameterCompute_More(SNES,void*,Vec,Vec,PetscReal*,PetscReal*);
      7: EXTERN PetscErrorCode DiffParameterDestroy_More(void*);

      9: typedef struct {  /* default context for matrix-free SNES */
     10:   SNES         snes;             /* SNES context */
     11:   Vec          w;                /* work vector */
     12:   MatNullSpace sp;               /* null space context */
     13:   PetscReal    error_rel;        /* square root of relative error in computing function */
     14:   PetscReal    umin;             /* minimum allowable u'a value relative to |u|_1 */
     15:   PetscTruth   jorge;            /* flag indicating use of Jorge's method for determining
     16:                                    the differencing parameter */
     17:   PetscReal    h;                /* differencing parameter */
     18:   PetscTruth   need_h;           /* flag indicating whether we must compute h */
     19:   PetscTruth   need_err;         /* flag indicating whether we must currently compute error_rel */
     20:   PetscTruth   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;

     28: PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat)
     29: {
     31:   MFCtx_Private  *ctx;

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

     44: /*
     45:    SNESMatrixFreeView2_Private - Views matrix-free parameters.
     46:  */
     47: PetscErrorCode SNESMatrixFreeView2_Private(Mat J,PetscViewer viewer)
     48: {
     50:   MFCtx_Private  *ctx;
     51:   PetscTruth     iascii;

     54:   MatShellGetContext(J,(void **)&ctx);
     55:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     56:   if (iascii) {
     57:      PetscViewerASCIIPrintf(viewer,"  SNES matrix-free approximation:\n");
     58:      if (ctx->jorge) {
     59:        PetscViewerASCIIPrintf(viewer,"    using Jorge's method of determining differencing parameter\n");
     60:      }
     61:      PetscViewerASCIIPrintf(viewer,"    err=%g (relative error in function evaluation)\n",ctx->error_rel);
     62:      PetscViewerASCIIPrintf(viewer,"    umin=%g (minimum iterate parameter)\n",ctx->umin);
     63:      if (ctx->compute_err) {
     64:        PetscViewerASCIIPrintf(viewer,"    freq_err=%D (frequency for computing err)\n",ctx->compute_err_freq);
     65:      }
     66:   } else {
     67:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by SNES matrix free Jorge",((PetscObject)viewer)->type_name);
     68:   }
     69:   return(0);
     70: }

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

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

    101:   PetscObjectGetComm((PetscObject)mat,&comm);
    102:   MatShellGetContext(mat,(void **)&ctx);
    103:   snes = ctx->snes;
    104:   w    = ctx->w;
    105:   umin = ctx->umin;

    107:   SNESGetSolution(snes,&U);
    108:   eval_fct = SNESComputeFunction;
    109:   SNESGetFunction(snes,&F,PETSC_NULL,PETSC_NULL);

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

    114:     /* Use Jorge's method to compute h */
    115:     if (ctx->jorge) {
    116:       DiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);

    118:     /* Use the Brown/Saad method to compute h */
    119:     } else {
    120:       /* Compute error if desired */
    121:       SNESGetIterationNumber(snes,&iter);
    122:       if ((ctx->need_err) || ((ctx->compute_err_freq) && (ctx->compute_err_iter != iter) && (!((iter-1)%ctx->compute_err_freq)))) {
    123:         /* Use Jorge's method to compute noise */
    124:         DiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);
    125:         ctx->error_rel = sqrt(noise);
    126:         PetscLogInfo((snes,"SNESMatrixFreeMult2_Private: Using Jorge's noise: noise=%g, sqrt(noise)=%g, h_more=%g\n",noise,ctx->error_rel,h));
    127:         ctx->compute_err_iter = iter;
    128:         ctx->need_err = PETSC_FALSE;
    129:       }

    131:       VecDotBegin(U,a,&dot);
    132:       VecNormBegin(a,NORM_1,&sum);
    133:       VecNormBegin(a,NORM_2,&norm);
    134:       VecDotEnd(U,a,&dot);
    135:       VecNormEnd(a,NORM_1,&sum);
    136:       VecNormEnd(a,NORM_2,&norm);

    139:       /* Safeguard for step sizes too small */
    140:       if (sum == 0.0) {dot = 1.0; norm = 1.0;}
    141: #if defined(PETSC_USE_COMPLEX)
    142:       else if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum;
    143:       else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum;
    144: #else
    145:       else if (dot < umin*sum && dot >= 0.0) dot = umin*sum;
    146:       else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum;
    147: #endif
    148:       h = PetscRealPart(ctx->error_rel*dot/(norm*norm));
    149:     }
    150:   } else {
    151:     h = ctx->h;
    152:   }
    153:   if (!ctx->jorge || !ctx->need_h) {PetscLogInfo((snes,"SNESMatrixFreeMult2_Private: h = %g\n",h));}

    155:   /* Evaluate function at F(u + ha) */
    156:   hs = h;
    157:   VecWAXPY(w,hs,a,U);
    158:   eval_fct(snes,w,y);
    159:   VecAXPY(y,mone,F);
    160:   hs = 1.0/hs;
    161:   VecScale(y,hs);
    162:   if (ctx->sp) {MatNullSpaceRemove(ctx->sp,y,PETSC_NULL);}

    164:   PetscLogEventEnd(MATSNESMF_Mult,a,y,0,0);
    165:   return(0);
    166: }

    170: /*@C
    171:    SNESMatrixFreeCreate2 - Creates a matrix-free matrix
    172:    context for use with a SNES solver.  This matrix can be used as
    173:    the Jacobian argument for the routine SNESSetJacobian().

    175:    Input Parameters:
    176: .  snes - the SNES context
    177: .  x - vector where SNES solution is to be stored.

    179:    Output Parameter:
    180: .  J - the matrix-free matrix

    182:    Level: advanced

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

    189: $       J(u)*a = [J(u+h*a) - J(u)]/h,
    190: $   where by default
    191: $        h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
    192: $          = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
    193: $   where
    194: $        error_rel = square root of relative error in
    195: $                    function evaluation
    196: $        umin = minimum iterate parameter
    197: $   Alternatively, the differencing parameter, h, can be set using
    198: $   Jorge's nifty new strategy if one specifies the option
    199: $          -snes_mf_jorge

    201:    The user can set these parameters via MatSNESMFSetFunctionError().
    202:    See the nonlinear solvers chapter of the users manual for details.

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

    207:    Options Database Keys:
    208: $  -snes_mf_err <error_rel>
    209: $  -snes_mf_unim <umin>
    210: $  -snes_mf_compute_err
    211: $  -snes_mf_freq_err <freq>
    212: $  -snes_mf_jorge

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

    216: .seealso: MatDestroy(), MatSNESMFSetFunctionError()
    217: @*/
    218: PetscErrorCode PETSCSNES_DLLEXPORT SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J)
    219: {
    220:   MPI_Comm       comm;
    221:   MFCtx_Private  *mfctx;
    223:   PetscInt       n,nloc;
    224:   PetscTruth     flg;
    225:   char           p[64];

    228:   PetscNew(MFCtx_Private,&mfctx);
    229:   PetscLogObjectMemory(snes,sizeof(MFCtx_Private));
    230:   mfctx->sp   = 0;
    231:   mfctx->snes = snes;
    232:   mfctx->error_rel        = PETSC_SQRT_MACHINE_EPSILON;
    233:   mfctx->umin             = 1.e-6;
    234:   mfctx->h                = 0.0;
    235:   mfctx->need_h           = PETSC_TRUE;
    236:   mfctx->need_err         = PETSC_FALSE;
    237:   mfctx->compute_err      = PETSC_FALSE;
    238:   mfctx->compute_err_freq = 0;
    239:   mfctx->compute_err_iter = -1;
    240:   PetscOptionsGetReal(snes->prefix,"-snes_mf_err",&mfctx->error_rel,PETSC_NULL);
    241:   PetscOptionsGetReal(snes->prefix,"-snes_mf_umin",&mfctx->umin,PETSC_NULL);
    242:   PetscOptionsHasName(snes->prefix,"-snes_mf_jorge",&mfctx->jorge);
    243:   PetscOptionsHasName(snes->prefix,"-snes_mf_compute_err",&mfctx->compute_err);
    244:   PetscOptionsGetInt(snes->prefix,"-snes_mf_freq_err",&mfctx->compute_err_freq,&flg);
    245:   if (flg) {
    246:     if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0;
    247:     mfctx->compute_err = PETSC_TRUE;
    248:   }
    249:   if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE;
    250:   if (mfctx->jorge || mfctx->compute_err) {
    251:     DiffParameterCreate_More(snes,x,&mfctx->data);
    252:   } else mfctx->data = 0;

    254:   PetscOptionsHasName(PETSC_NULL,"-help",&flg);
    255:   PetscStrcpy(p,"-");
    256:   if (snes->prefix) PetscStrcat(p,snes->prefix);
    257:   if (flg) {
    258:     PetscPrintf(snes->comm," Matrix-free Options (via SNES):\n");
    259:     PetscPrintf(snes->comm,"   %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n",p,mfctx->error_rel);
    260:     PetscPrintf(snes->comm,"   %ssnes_mf_umin <umin>: see users manual (default %g)\n",p,mfctx->umin);
    261:     PetscPrintf(snes->comm,"   %ssnes_mf_jorge: use Jorge More's method\n",p);
    262:     PetscPrintf(snes->comm,"   %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);
    263:     PetscPrintf(snes->comm,"   %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);
    264:     PetscPrintf(snes->comm,"   %ssnes_mf_noise_file <file>: set file for printing noise info\n",p);
    265:   }
    266:   VecDuplicate(x,&mfctx->w);
    267:   PetscObjectGetComm((PetscObject)x,&comm);
    268:   VecGetSize(x,&n);
    269:   VecGetLocalSize(x,&nloc);
    270:   MatCreate(comm,J);
    271:   MatSetSizes(*J,nloc,n,n,n);
    272:   MatSetType(*J,MATSHELL);
    273:   MatShellSetContext(*J,mfctx);
    274:   MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))SNESMatrixFreeMult2_Private);
    275:   MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))SNESMatrixFreeDestroy2_Private);
    276:   MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))SNESMatrixFreeView2_Private);

    278:   PetscLogObjectParent(*J,mfctx->w);
    279:   PetscLogObjectParent(snes,*J);
    280:   return(0);
    281: }

    285: /*@C
    286:    SNESDefaultMatrixFreeSetParameters2 - Sets the parameters for the approximation of
    287:    matrix-vector products using finite differences.

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

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

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

    297:    Input Parameters:
    298: +  mat - the matrix
    299: .  error_rel - relative error (should be set to the square root of
    300:      the relative error in the function evaluations)
    301: .  umin - minimum allowable u-value
    302: -  h - differencing parameter

    304:    Level: advanced

    306:    Notes:
    307:    If the user sets the parameter h directly, then this value will be used
    308:    instead of the default computation indicated above.

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

    312: .seealso: MatCreateSNESMF()
    313: @*/
    314: PetscErrorCode PETSCSNES_DLLEXPORT SNESDefaultMatrixFreeSetParameters2(Mat mat,PetscReal error,PetscReal umin,PetscReal h)
    315: {
    316:   MFCtx_Private  *ctx;

    320:   MatShellGetContext(mat,(void **)&ctx);
    321:   if (ctx) {
    322:     if (error != PETSC_DEFAULT) ctx->error_rel = error;
    323:     if (umin  != PETSC_DEFAULT) ctx->umin = umin;
    324:     if (h     != PETSC_DEFAULT) {
    325:       ctx->h = h;
    326:       ctx->need_h = PETSC_FALSE;
    327:     }
    328:   }
    329:   return(0);
    330: }

    332: PetscErrorCode PETSCSNES_DLLEXPORT SNESUnSetMatrixFreeParameter(SNES snes)
    333: {
    334:   MFCtx_Private  *ctx;
    336:   Mat            mat;

    339:   SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL,PETSC_NULL);
    340:   MatShellGetContext(mat,(void **)&ctx);
    341:   if (ctx) ctx->need_h = PETSC_TRUE;
    342:   return(0);
    343: }
    344: