Actual source code: snesmfj.c

      1: #define PETSCSNES_DLL

     3:  #include src/mat/matimpl.h
     4:  #include src/snes/mf/snesmfj.h

      6: PetscFList MatSNESMPetscFList         = 0;
      7: PetscTruth MatSNESMFRegisterAllCalled = PETSC_FALSE;

      9: PetscCookie PETSCSNES_DLLEXPORT MATSNESMFCTX_COOKIE = 0;
     10: PetscEvent  MATSNESMF_Mult = 0;

     14: /*@C
     15:     MatSNESMFSetType - Sets the method that is used to compute the
     16:     differencing parameter for finite differene matrix-free formulations.

     18:     Input Parameters:
     19: +   mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMF()
     20:           or MatSetType(mat,MATMFFD);
     21: -   ftype - the type requested

     23:     Level: advanced

     25:     Notes:
     26:     For example, such routines can compute h for use in
     27:     Jacobian-vector products of the form

     29:                         F(x+ha) - F(x)
     30:           F'(u)a  ~=  ----------------
     31:                               h

     33: .seealso: MatCreateSNESMF(), MatSNESMFRegisterDynamic)
     34: @*/
     35: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetType(Mat mat,const MatSNESMFType ftype)
     36: {
     37:   PetscErrorCode ierr,(*r)(MatSNESMFCtx);
     38:   MatSNESMFCtx   ctx = (MatSNESMFCtx)mat->data;
     39:   PetscTruth     match;
     40:

     45:   /* already set, so just return */
     46:   PetscTypeCompare((PetscObject)ctx,ftype,&match);
     47:   if (match) return(0);

     49:   /* destroy the old one if it exists */
     50:   if (ctx->ops->destroy) {
     51:     (*ctx->ops->destroy)(ctx);
     52:   }

     54:   /* Get the function pointers for the requrested method */
     55:   if (!MatSNESMFRegisterAllCalled) {MatSNESMFRegisterAll(PETSC_NULL);}
     56:    PetscFListFind(ctx->comm,MatSNESMPetscFList,ftype,(void (**)(void)) &r);
     57:   if (!r) SETERRQ1(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatSNESMF type %s given",ftype);
     58:   (*r)(ctx);
     59:   PetscObjectChangeTypeName((PetscObject)ctx,ftype);
     60:   return(0);
     61: }

     67: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetFunctioniBase_FD(Mat mat,FCN1 func)
     68: {
     69:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

     72:   ctx->funcisetbase = func;
     73:   return(0);
     74: }

     81: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetFunctioni_FD(Mat mat,FCN2 funci)
     82: {
     83:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

     86:   ctx->funci = funci;
     87:   return(0);
     88: }

     94: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(MatSNESMFCtx))
     95: {
     97:   char           fullname[PETSC_MAX_PATH_LEN];

    100:   PetscFListConcat(path,name,fullname);
    101:   PetscFListAdd(&MatSNESMPetscFList,sname,fullname,(void (*)(void))function);
    102:   return(0);
    103: }

    108: /*@C
    109:    MatSNESMFRegisterDestroy - Frees the list of MatSNESMF methods that were
    110:    registered by MatSNESMFRegisterDynamic).

    112:    Not Collective

    114:    Level: developer

    116: .keywords: MatSNESMF, register, destroy

    118: .seealso: MatSNESMFRegisterDynamic), MatSNESMFRegisterAll()
    119: @*/
    120: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFRegisterDestroy(void)
    121: {

    125:   if (MatSNESMPetscFList) {
    126:     PetscFListDestroy(&MatSNESMPetscFList);
    127:     MatSNESMPetscFList = 0;
    128:   }
    129:   MatSNESMFRegisterAllCalled = PETSC_FALSE;
    130:   return(0);
    131: }

    133: /* ----------------------------------------------------------------------------------------*/
    136: PetscErrorCode MatDestroy_MFFD(Mat mat)
    137: {
    139:   MatSNESMFCtx   ctx = (MatSNESMFCtx)mat->data;

    142:   if (ctx->w) {
    143:     VecDestroy(ctx->w);
    144:   }
    145:   if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
    146:   if (ctx->sp) {MatNullSpaceDestroy(ctx->sp);}
    147:   PetscHeaderDestroy(ctx);

    149:   PetscObjectComposeFunction((PetscObject)mat,"MatSNESMFSetBase_C","",PETSC_NULL);
    150:   PetscObjectComposeFunction((PetscObject)mat,"MatSNESMFSetFunctioniBase_C","",PETSC_NULL);
    151:   PetscObjectComposeFunction((PetscObject)mat,"MatSNESMFSetFunctioni_C","",PETSC_NULL);
    152:   PetscObjectComposeFunction((PetscObject)mat,"MatSNESMFSetCheckh_C","",PETSC_NULL);

    154:   return(0);
    155: }

    159: /*
    160:    MatSNESMFView_MFFD - Views matrix-free parameters.

    162: */
    163: PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
    164: {
    166:   MatSNESMFCtx   ctx = (MatSNESMFCtx)J->data;
    167:   PetscTruth     iascii;

    170:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    171:   if (iascii) {
    172:      PetscViewerASCIIPrintf(viewer,"  SNES matrix-free approximation:\n");
    173:      PetscViewerASCIIPrintf(viewer,"    err=%g (relative error in function evaluation)\n",ctx->error_rel);
    174:      if (!ctx->type_name) {
    175:        PetscViewerASCIIPrintf(viewer,"    The compute h routine has not yet been set\n");
    176:      } else {
    177:        PetscViewerASCIIPrintf(viewer,"    Using %s compute h routine\n",ctx->type_name);
    178:      }
    179:      if (ctx->ops->view) {
    180:        (*ctx->ops->view)(ctx,viewer);
    181:      }
    182:   } else {
    183:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES matrix free matrix",((PetscObject)viewer)->type_name);
    184:   }
    185:   return(0);
    186: }

    190: /*
    191:    MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
    192:    allows the user to indicate the beginning of a new linear solve by calling
    193:    MatAssemblyXXX() on the matrix free matrix. This then allows the
    194:    MatSNESMFCreate_WP() to properly compute ||U|| only the first time
    195:    in the linear solver rather than every time.
    196: */
    197: PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
    198: {
    200:   MatSNESMFCtx   j = (MatSNESMFCtx)J->data;

    203:   MatSNESMFResetHHistory(J);
    204:   if (j->usesnes) {
    205:     SNESGetSolution(j->snes,&j->current_u);
    206:     SNESGetFunction(j->snes,&j->current_f,PETSC_NULL,PETSC_NULL);
    207:     if (!j->w) {
    208:       VecDuplicate(j->current_u, &j->w);
    209:     }
    210:   }
    211:   j->vshift = 0.0;
    212:   j->vscale = 1.0;
    213:   return(0);
    214: }

    218: /*
    219:   MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:

    221:         y ~= (F(u + ha) - F(u))/h,
    222:   where F = nonlinear function, as set by SNESSetFunction()
    223:         u = current iterate
    224:         h = difference interval
    225: */
    226: PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
    227: {
    228:   MatSNESMFCtx   ctx = (MatSNESMFCtx)mat->data;
    229:   SNES           snes;
    230:   PetscScalar    h,mone = -1.0;
    231:   Vec            w,U,F;
    232:   PetscErrorCode ierr,(*eval_fct)(SNES,Vec,Vec)=0;
    233:   PetscTruth     zeroa;

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

    242:   snes = ctx->snes;
    243:   w    = ctx->w;
    244:   U    = ctx->current_u;

    246:   /*
    247:       Compute differencing parameter
    248:   */
    249:   if (!ctx->ops->compute) {
    250:     MatSNESMFSetType(mat,MATSNESMF_WP);
    251:     MatSNESMFSetFromOptions(mat);
    252:   }
    253:   (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
    254:   if (zeroa) {
    255:     PetscScalar zero = 0.0;
    256:     VecSet(y,zero);
    257:     return(0);
    258:   }

    260:   if (ctx->checkh) {
    261:     (*ctx->checkh)(U,a,&h,ctx->checkhctx);
    262:   }

    264:   /* keep a record of the current differencing parameter h */
    265:   ctx->currenth = h;
    266: #if defined(PETSC_USE_COMPLEX)
    267:   PetscLogInfo((mat,"MatMult_MFFD:Current differencing parameter: %g + %g i\n",PetscRealPart(h),PetscImaginaryPart(h)));
    268: #else
    269:   PetscLogInfo((mat,"MatMult_MFFD:Current differencing parameter: %15.12e\n",h));
    270: #endif
    271:   if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
    272:     ctx->historyh[ctx->ncurrenth] = h;
    273:   }
    274:   ctx->ncurrenth++;

    276:   /* w = u + ha */
    277:   VecWAXPY(w,h,a,U);

    279:   if (ctx->usesnes) {
    280:     eval_fct = SNESComputeFunction;
    281:     F    = ctx->current_f;
    282:     if (!F) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"You must call MatAssembly() even on matrix-free matrices");
    283:     (*eval_fct)(snes,w,y);
    284:   } else {
    285:     F = ctx->funcvec;
    286:     /* compute func(U) as base for differencing */
    287:     if (ctx->ncurrenth == 1) {
    288:       (*ctx->func)(snes,U,F,ctx->funcctx);
    289:     }
    290:     (*ctx->func)(snes,w,y,ctx->funcctx);
    291:   }

    293:   VecAXPY(y,mone,F);
    294:   h    = 1.0/h;
    295:   VecScale(y,h);

    297:   VecAXPBY(y,ctx->vshift,ctx->vscale,a);

    299:   if (ctx->sp) {MatNullSpaceRemove(ctx->sp,y,PETSC_NULL);}

    301:   PetscLogEventEnd(MATSNESMF_Mult,a,y,0,0);
    302:   return(0);
    303: }

    307: /*
    308:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix

    310:         y ~= (F(u + ha) - F(u))/h,
    311:   where F = nonlinear function, as set by SNESSetFunction()
    312:         u = current iterate
    313:         h = difference interval
    314: */
    315: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
    316: {
    317:   MatSNESMFCtx   ctx = (MatSNESMFCtx)mat->data;
    318:   PetscScalar    h,*aa,*ww,v;
    319:   PetscReal      epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
    320:   Vec            w,U;
    322:   PetscInt       i,rstart,rend;

    325:   if (!ctx->funci) {
    326:     SETERRQ(PETSC_ERR_ORDER,"Requires calling MatSNESMFSetFunctioni() first");
    327:   }

    329:   w    = ctx->w;
    330:   U    = ctx->current_u;
    331:   (*ctx->func)(0,U,a,ctx->funcctx);
    332:   (*ctx->funcisetbase)(U,ctx->funcctx);
    333:   VecCopy(U,w);

    335:   VecGetOwnershipRange(a,&rstart,&rend);
    336:   VecGetArray(a,&aa);
    337:   for (i=rstart; i<rend; i++) {
    338:     VecGetArray(w,&ww);
    339:     h  = ww[i-rstart];
    340:     if (h == 0.0) h = 1.0;
    341: #if !defined(PETSC_USE_COMPLEX)
    342:     if (h < umin && h >= 0.0)      h = umin;
    343:     else if (h < 0.0 && h > -umin) h = -umin;
    344: #else
    345:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
    346:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
    347: #endif
    348:     h     *= epsilon;
    349:
    350:     ww[i-rstart] += h;
    351:     VecRestoreArray(w,&ww);
    352:     (*ctx->funci)(i,w,&v,ctx->funcctx);
    353:     aa[i-rstart]  = (v - aa[i-rstart])/h;

    355:     /* possibly shift and scale result */
    356:     aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];

    358:     VecGetArray(w,&ww);
    359:     ww[i-rstart] -= h;
    360:     VecRestoreArray(w,&ww);
    361:   }
    362:   VecRestoreArray(a,&aa);
    363:   return(0);
    364: }

    368: PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
    369: {
    370:   MatSNESMFCtx shell = (MatSNESMFCtx)Y->data;
    372:   shell->vshift += a;
    373:   return(0);
    374: }

    378: PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
    379: {
    380:   MatSNESMFCtx shell = (MatSNESMFCtx)Y->data;
    382:   shell->vscale *= a;
    383:   return(0);
    384: }

    389: /*@C
    390:    MatCreateSNESMF - Creates a matrix-free matrix context for use with
    391:    a SNES solver.  This matrix can be used as the Jacobian argument for
    392:    the routine SNESSetJacobian().

    394:    Collective on SNES and Vec

    396:    Input Parameters:
    397: +  snes - the SNES context
    398: -  x - vector where SNES solution is to be stored.

    400:    Output Parameter:
    401: .  J - the matrix-free matrix

    403:    Level: advanced

    405:    Notes:
    406:    The matrix-free matrix context merely contains the function pointers
    407:    and work space for performing finite difference approximations of
    408:    Jacobian-vector products, F'(u)*a,

    410:    The default code uses the following approach to compute h

    412: .vb
    413:      F'(u)*a = [F(u+h*a) - F(u)]/h where
    414:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
    415:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
    416:  where
    417:      error_rel = square root of relative error in function evaluation
    418:      umin = minimum iterate parameter
    419: .ve

    421:    The user can set the error_rel via MatSNESMFSetFunctionError() and
    422:    umin via MatSNESMFDefaultSetUmin(); see the nonlinear solvers chapter
    423:    of the users manual for details.

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

    428:    Options Database Keys:
    429: +  -snes_mf_err <error_rel> - Sets error_rel
    430: .  -snes_mf_unim <umin> - Sets umin (for default PETSc routine that computes h only)
    431: -  -snes_mf_ksp_monitor - KSP monitor routine that prints differencing h

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

    435: .seealso: MatDestroy(), MatSNESMFSetFunctionError(), MatSNESMFDefaultSetUmin()
    436:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(), MatCreateMF(),
    437:           MatSNESMFGetH(),MatSNESMFKSPMonitor(), MatSNESMFRegisterDynamic), MatSNESMFComputeJacobian()
    438:
    439: @*/
    440: PetscErrorCode PETSCSNES_DLLEXPORT MatCreateSNESMF(SNES snes,Vec x,Mat *J)
    441: {
    442:   MatSNESMFCtx   mfctx;

    446:   MatCreateMF(x,J);

    448:   mfctx          = (MatSNESMFCtx)(*J)->data;
    449:   mfctx->snes    = snes;
    450:   mfctx->usesnes = PETSC_TRUE;
    451:   PetscLogObjectParent(snes,*J);
    452:   return(0);
    453: }

    458: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetBase_FD(Mat J,Vec U)
    459: {
    461:   MatSNESMFCtx   ctx = (MatSNESMFCtx)J->data;

    464:   MatSNESMFResetHHistory(J);
    465:   ctx->current_u = U;
    466:   ctx->usesnes   = PETSC_FALSE;
    467:   if (!ctx->w) {
    468:     VecDuplicate(ctx->current_u, &ctx->w);
    469:   }
    470:   J->assembled = PETSC_TRUE;
    471:   return(0);
    472: }

    479: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetCheckh_FD(Mat J,FCN3 fun,void*ectx)
    480: {
    481:   MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;

    484:   ctx->checkh    = fun;
    485:   ctx->checkhctx = ectx;
    486:   return(0);
    487: }

    492: /*@
    493:    MatSNESMFSetFromOptions - Sets the MatSNESMF options from the command line
    494:    parameter.

    496:    Collective on Mat

    498:    Input Parameters:
    499: .  mat - the matrix obtained with MatCreateSNESMF()

    501:    Options Database Keys:
    502: +  -snes_mf_type - <default,wp>
    503: -  -snes_mf_err - square root of estimated relative error in function evaluation
    504: -  -snes_mf_period - how often h is recomputed, defaults to 1, everytime

    506:    Level: advanced

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

    510: .seealso: MatCreateSNESMF(),MatSNESMFSetHHistory(),
    511:           MatSNESMFResetHHistory(), MatSNESMFKSPMonitor()
    512: @*/
    513: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetFromOptions(Mat mat)
    514: {
    515:   MatSNESMFCtx   mfctx = (MatSNESMFCtx)mat->data;
    517:   PetscTruth     flg;
    518:   char           ftype[256];

    521:   if (!MatSNESMFRegisterAllCalled) {MatSNESMFRegisterAll(PETSC_NULL);}
    522:
    523:   PetscOptionsBegin(mfctx->comm,mfctx->prefix,"Set matrix free computation parameters","MatSNESMF");
    524:   PetscOptionsList("-snes_mf_type","Matrix free type","MatSNESMFSetType",MatSNESMPetscFList,mfctx->type_name,ftype,256,&flg);
    525:   if (flg) {
    526:     MatSNESMFSetType(mat,ftype);
    527:   }

    529:   PetscOptionsReal("-snes_mf_err","set sqrt relative error in function","MatSNESMFSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
    530:   PetscOptionsInt("-snes_mf_period","how often h is recomputed","MatSNESMFSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
    531:   if (mfctx->snes) {
    532:     PetscOptionsName("-snes_mf_ksp_monitor","Monitor matrix-free parameters","MatSNESMFKSPMonitor",&flg);
    533:     if (flg) {
    534:       KSP ksp;
    535:       SNESGetKSP(mfctx->snes,&ksp);
    536:       KSPSetMonitor(ksp,MatSNESMFKSPMonitor,PETSC_NULL,0);
    537:     }
    538:   }
    539:   PetscOptionsName("-snes_mf_check_positivity","Insure that U + h*a is nonnegative","MatSNESMFSetCheckh",&flg);
    540:   if (flg) {
    541:     MatSNESMFSetCheckh(mat,MatSNESMFCheckPositivity,0);
    542:   }
    543:   if (mfctx->ops->setfromoptions) {
    544:     (*mfctx->ops->setfromoptions)(mfctx);
    545:   }
    546:   PetscOptionsEnd();
    547:   return(0);
    548: }

    550: /*MC
    551:   MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.

    553:   Level: advanced

    555: .seealso: MatCreateMF(), MatCreateSNESMF()
    556: M*/
    560: PetscErrorCode PETSCSNES_DLLEXPORT MatCreate_MFFD(Mat A)
    561: {
    562:   MatSNESMFCtx mfctx;

    566: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    567:   SNESInitializePackage(PETSC_NULL);
    568: #endif

    570:   PetscHeaderCreate(mfctx,_p_MatSNESMFCtx,struct _MFOps,MATSNESMFCTX_COOKIE,0,"SNESMF",A->comm,MatDestroy_MFFD,MatView_MFFD);
    571:   mfctx->sp              = 0;
    572:   mfctx->snes            = 0;
    573:   mfctx->error_rel       = PETSC_SQRT_MACHINE_EPSILON;
    574:   mfctx->recomputeperiod = 1;
    575:   mfctx->count           = 0;
    576:   mfctx->currenth        = 0.0;
    577:   mfctx->historyh        = PETSC_NULL;
    578:   mfctx->ncurrenth       = 0;
    579:   mfctx->maxcurrenth     = 0;
    580:   mfctx->type_name       = 0;
    581:   mfctx->usesnes         = PETSC_FALSE;

    583:   mfctx->vshift          = 0.0;
    584:   mfctx->vscale          = 1.0;

    586:   /*
    587:      Create the empty data structure to contain compute-h routines.
    588:      These will be filled in below from the command line options or
    589:      a later call with MatSNESMFSetType() or if that is not called
    590:      then it will default in the first use of MatMult_MFFD()
    591:   */
    592:   mfctx->ops->compute        = 0;
    593:   mfctx->ops->destroy        = 0;
    594:   mfctx->ops->view           = 0;
    595:   mfctx->ops->setfromoptions = 0;
    596:   mfctx->hctx                = 0;

    598:   mfctx->func                = 0;
    599:   mfctx->funcctx             = 0;
    600:   mfctx->funcvec             = 0;
    601:   mfctx->w                   = PETSC_NULL;

    603:   A->data                = mfctx;

    605:   A->ops->mult           = MatMult_MFFD;
    606:   A->ops->destroy        = MatDestroy_MFFD;
    607:   A->ops->view           = MatView_MFFD;
    608:   A->ops->assemblyend    = MatAssemblyEnd_MFFD;
    609:   A->ops->getdiagonal    = MatGetDiagonal_MFFD;
    610:   A->ops->scale          = MatScale_MFFD;
    611:   A->ops->shift          = MatShift_MFFD;
    612:   A->ops->setfromoptions = MatSNESMFSetFromOptions;
    613:   A->assembled = PETSC_TRUE;

    615:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetBase_C","MatSNESMFSetBase_FD",MatSNESMFSetBase_FD);
    616:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetFunctioniBase_C","MatSNESMFSetFunctioniBase_FD",MatSNESMFSetFunctioniBase_FD);
    617:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetFunctioni_C","MatSNESMFSetFunctioni_FD",MatSNESMFSetFunctioni_FD);
    618:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetCheckh_C","MatSNESMFSetCheckh_FD",MatSNESMFSetCheckh_FD);
    619:   mfctx->mat = A;

    621:   return(0);
    622: }

    627: /*@C
    628:    MatCreateMF - Creates a matrix-free matrix. See also MatCreateSNESMF()

    630:    Collective on Vec

    632:    Input Parameters:
    633: .  x - vector that defines layout of the vectors and matrices

    635:    Output Parameter:
    636: .  J - the matrix-free matrix

    638:    Level: advanced

    640:    Notes:
    641:    The matrix-free matrix context merely contains the function pointers
    642:    and work space for performing finite difference approximations of
    643:    Jacobian-vector products, F'(u)*a,

    645:    The default code uses the following approach to compute h

    647: .vb
    648:      F'(u)*a = [F(u+h*a) - F(u)]/h where
    649:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
    650:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
    651:  where
    652:      error_rel = square root of relative error in function evaluation
    653:      umin = minimum iterate parameter
    654: .ve

    656:    The user can set the error_rel via MatSNESMFSetFunctionError() and
    657:    umin via MatSNESMFDefaultSetUmin(); see the nonlinear solvers chapter
    658:    of the users manual for details.

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

    663:    Options Database Keys:
    664: +  -snes_mf_err <error_rel> - Sets error_rel
    665: .  -snes_mf_unim <umin> - Sets umin (for default PETSc routine that computes h only)
    666: .  -snes_mf_ksp_monitor - KSP monitor routine that prints differencing h
    667: -  -snes_mf_check_positivity

    669: .keywords: default, matrix-free, create, matrix

    671: .seealso: MatDestroy(), MatSNESMFSetFunctionError(), MatSNESMFDefaultSetUmin()
    672:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(), MatCreateSNESMF(),
    673:           MatSNESMFGetH(),MatSNESMFKSPMonitor(), MatSNESMFRegisterDynamic),, MatSNESMFComputeJacobian()
    674:
    675: @*/
    676: PetscErrorCode PETSCSNES_DLLEXPORT MatCreateMF(Vec x,Mat *J)
    677: {
    678:   MPI_Comm       comm;
    680:   PetscInt       n,nloc;

    683:   PetscObjectGetComm((PetscObject)x,&comm);
    684:   VecGetSize(x,&n);
    685:   VecGetLocalSize(x,&nloc);
    686:   MatCreate(comm,J);
    687:   MatSetSizes(*J,nloc,nloc,n,n);
    688:   MatRegisterDynamic(MATMFFD,0,"MatCreate_MFFD",MatCreate_MFFD);
    689:   MatSetType(*J,MATMFFD);
    690:   return(0);
    691: }

    696: /*@
    697:    MatSNESMFGetH - Gets the last value that was used as the differencing
    698:    parameter.

    700:    Not Collective

    702:    Input Parameters:
    703: .  mat - the matrix obtained with MatCreateSNESMF()

    705:    Output Paramter:
    706: .  h - the differencing step size

    708:    Level: advanced

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

    712: .seealso: MatCreateSNESMF(),MatSNESMFSetHHistory(),
    713:           MatSNESMFResetHHistory(),MatSNESMFKSPMonitor()
    714: @*/
    715: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFGetH(Mat mat,PetscScalar *h)
    716: {
    717:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

    720:   *h = ctx->currenth;
    721:   return(0);
    722: }

    726: /*
    727:    MatSNESMFKSPMonitor - A KSP monitor for use with the default PETSc
    728:    SNES matrix free routines. Prints the differencing parameter used at
    729:    each step.
    730: */
    731: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFKSPMonitor(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
    732: {
    733:   PC             pc;
    734:   MatSNESMFCtx   ctx;
    736:   Mat            mat;
    737:   MPI_Comm       comm;
    738:   PetscTruth     nonzeroinitialguess;

    741:   PetscObjectGetComm((PetscObject)ksp,&comm);
    742:   KSPGetPC(ksp,&pc);
    743:   KSPGetInitialGuessNonzero(ksp,&nonzeroinitialguess);
    744:   PCGetOperators(pc,&mat,PETSC_NULL,PETSC_NULL);
    745:   ctx  = (MatSNESMFCtx)mat->data;

    747:   if (n > 0 || nonzeroinitialguess) {
    748: #if defined(PETSC_USE_COMPLEX)
    749:     PetscPrintf(comm,"%D KSP Residual norm %14.12e h %g + %g i\n",n,rnorm,
    750:                 PetscRealPart(ctx->currenth),PetscImaginaryPart(ctx->currenth));
    751: #else
    752:     PetscPrintf(comm,"%D KSP Residual norm %14.12e h %g \n",n,rnorm,ctx->currenth);
    753: #endif
    754:   } else {
    755:     PetscPrintf(comm,"%D KSP Residual norm %14.12e\n",n,rnorm);
    756:   }
    757:   return(0);
    758: }

    762: /*@C
    763:    MatSNESMFSetFunction - Sets the function used in applying the matrix free.

    765:    Collective on Mat

    767:    Input Parameters:
    768: +  mat - the matrix free matrix created via MatCreateSNESMF()
    769: .  v   - workspace vector
    770: .  func - the function to use
    771: -  funcctx - optional function context passed to function

    773:    Level: advanced

    775:    Notes:
    776:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
    777:     matrix inside your compute Jacobian routine

    779:     If this is not set then it will use the function set with SNESSetFunction()

    781: .keywords: SNES, matrix-free, function

    783: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
    784:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
    785:           MatSNESMFKSPMonitor(), SNESetFunction()
    786: @*/
    787: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetFunction(Mat mat,Vec v,PetscErrorCode (*func)(SNES,Vec,Vec,void *),void *funcctx)
    788: {
    789:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

    792:   ctx->func    = func;
    793:   ctx->funcctx = funcctx;
    794:   ctx->funcvec = v;
    795:   return(0);
    796: }

    800: /*@C
    801:    MatSNESMFSetFunctioni - Sets the function for a single component

    803:    Collective on Mat

    805:    Input Parameters:
    806: +  mat - the matrix free matrix created via MatCreateSNESMF()
    807: -  funci - the function to use

    809:    Level: advanced

    811:    Notes:
    812:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
    813:     matrix inside your compute Jacobian routine

    816: .keywords: SNES, matrix-free, function

    818: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
    819:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
    820:           MatSNESMFKSPMonitor(), SNESetFunction()
    821: @*/
    822: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetFunctioni(Mat mat,PetscErrorCode (*funci)(PetscInt,Vec,PetscScalar*,void *))
    823: {
    824:   PetscErrorCode ierr,(*f)(Mat,PetscErrorCode (*)(PetscInt,Vec,PetscScalar*,void *));

    828:   PetscObjectQueryFunction((PetscObject)mat,"MatSNESMFSetFunctioni_C",(void (**)(void))&f);
    829:   if (f) {
    830:     (*f)(mat,funci);
    831:   }
    832:   return(0);
    833: }

    838: /*@C
    839:    MatSNESMFSetFunctioniBase - Sets the base vector for a single component function evaluation

    841:    Collective on Mat

    843:    Input Parameters:
    844: +  mat - the matrix free matrix created via MatCreateSNESMF()
    845: -  func - the function to use

    847:    Level: advanced

    849:    Notes:
    850:     If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
    851:     matrix inside your compute Jacobian routine

    854: .keywords: SNES, matrix-free, function

    856: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
    857:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
    858:           MatSNESMFKSPMonitor(), SNESetFunction()
    859: @*/
    860: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetFunctioniBase(Mat mat,PetscErrorCode (*func)(Vec,void *))
    861: {
    862:   PetscErrorCode ierr,(*f)(Mat,PetscErrorCode (*)(Vec,void *));

    866:   PetscObjectQueryFunction((PetscObject)mat,"MatSNESMFSetFunctioniBase_C",(void (**)(void))&f);
    867:   if (f) {
    868:     (*f)(mat,func);
    869:   }
    870:   return(0);
    871: }

    876: /*@
    877:    MatSNESMFSetPeriod - Sets how often h is recomputed, by default it is everytime

    879:    Collective on Mat

    881:    Input Parameters:
    882: +  mat - the matrix free matrix created via MatCreateSNESMF()
    883: -  period - 1 for everytime, 2 for every second etc

    885:    Options Database Keys:
    886: +  -snes_mf_period <period>

    888:    Level: advanced

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

    893: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
    894:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
    895:           MatSNESMFKSPMonitor()
    896: @*/
    897: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetPeriod(Mat mat,PetscInt period)
    898: {
    899:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

    902:   ctx->recomputeperiod = period;
    903:   return(0);
    904: }

    908: /*@
    909:    MatSNESMFSetFunctionError - Sets the error_rel for the approximation of
    910:    matrix-vector products using finite differences.

    912:    Collective on Mat

    914:    Input Parameters:
    915: +  mat - the matrix free matrix created via MatCreateSNESMF()
    916: -  error_rel - relative error (should be set to the square root of
    917:                the relative error in the function evaluations)

    919:    Options Database Keys:
    920: +  -snes_mf_err <error_rel> - Sets error_rel

    922:    Level: advanced

    924:    Notes:
    925:    The default matrix-free matrix-vector product routine computes
    926: .vb
    927:      F'(u)*a = [F(u+h*a) - F(u)]/h where
    928:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
    929:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
    930: .ve

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

    934: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
    935:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
    936:           MatSNESMFKSPMonitor()
    937: @*/
    938: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetFunctionError(Mat mat,PetscReal error)
    939: {
    940:   MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;

    943:   if (error != PETSC_DEFAULT) ctx->error_rel = error;
    944:   return(0);
    945: }

    949: /*@
    950:    MatSNESMFAddNullSpace - Provides a null space that an operator is
    951:    supposed to have.  Since roundoff will create a small component in
    952:    the null space, if you know the null space you may have it
    953:    automatically removed.

    955:    Collective on Mat

    957:    Input Parameters:
    958: +  J - the matrix-free matrix context
    959: -  nullsp - object created with MatNullSpaceCreate()

    961:    Level: advanced

    963: .keywords: SNES, matrix-free, null space

    965: .seealso: MatNullSpaceCreate(), MatSNESMFGetH(), MatCreateSNESMF(),
    966:           MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
    967:           MatSNESMFKSPMonitor(), MatSNESMFErrorRel()
    968: @*/
    969: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFAddNullSpace(Mat J,MatNullSpace nullsp)
    970: {
    972:   MatSNESMFCtx   ctx = (MatSNESMFCtx)J->data;
    973:   MPI_Comm       comm;

    976:   PetscObjectGetComm((PetscObject)J,&comm);

    978:   ctx->sp = nullsp;
    979:   PetscObjectReference((PetscObject)nullsp);
    980:   return(0);
    981: }

    985: /*@
    986:    MatSNESMFSetHHistory - Sets an array to collect a history of the
    987:    differencing values (h) computed for the matrix-free product.

    989:    Collective on Mat

    991:    Input Parameters:
    992: +  J - the matrix-free matrix context
    993: .  histroy - space to hold the history
    994: -  nhistory - number of entries in history, if more entries are generated than
    995:               nhistory, then the later ones are discarded

    997:    Level: advanced

    999:    Notes:
    1000:    Use MatSNESMFResetHHistory() to reset the history counter and collect
    1001:    a new batch of differencing parameters, h.

    1003: .keywords: SNES, matrix-free, h history, differencing history

    1005: .seealso: MatSNESMFGetH(), MatCreateSNESMF(),
    1006:           MatSNESMFResetHHistory(),
    1007:           MatSNESMFKSPMonitor(), MatSNESMFSetFunctionError()

    1009: @*/
    1010: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
    1011: {
    1012:   MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;

    1015:   ctx->historyh    = history;
    1016:   ctx->maxcurrenth = nhistory;
    1017:   ctx->currenth    = 0;
    1018:   return(0);
    1019: }

    1023: /*@
    1024:    MatSNESMFResetHHistory - Resets the counter to zero to begin
    1025:    collecting a new set of differencing histories.

    1027:    Collective on Mat

    1029:    Input Parameters:
    1030: .  J - the matrix-free matrix context

    1032:    Level: advanced

    1034:    Notes:
    1035:    Use MatSNESMFSetHHistory() to create the original history counter.

    1037: .keywords: SNES, matrix-free, h history, differencing history

    1039: .seealso: MatSNESMFGetH(), MatCreateSNESMF(),
    1040:           MatSNESMFSetHHistory(),
    1041:           MatSNESMFKSPMonitor(), MatSNESMFSetFunctionError()

    1043: @*/
    1044: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFResetHHistory(Mat J)
    1045: {
    1046:   MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;

    1049:   ctx->ncurrenth    = 0;
    1050:   return(0);
    1051: }

    1055: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFComputeJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    1056: {
    1059:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    1060:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
    1061:   return(0);
    1062: }

    1066: /*@
    1067:     MatSNESMFSetBase - Sets the vector U at which matrix vector products of the
    1068:         Jacobian are computed

    1070:     Collective on Mat

    1072:     Input Parameters:
    1073: +   J - the MatSNESMF matrix
    1074: -   U - the vector

    1076:     Notes: This is rarely used directly

    1078:     Level: advanced

    1080: @*/
    1081: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetBase(Mat J,Vec U)
    1082: {
    1083:   PetscErrorCode ierr,(*f)(Mat,Vec);

    1088:   PetscObjectQueryFunction((PetscObject)J,"MatSNESMFSetBase_C",(void (**)(void))&f);
    1089:   if (f) {
    1090:     (*f)(J,U);
    1091:   }
    1092:   return(0);
    1093: }

    1097: /*@C
    1098:     MatSNESMFSetCheckh - Sets a function that checks the computed h and adjusts
    1099:         it to satisfy some criteria

    1101:     Collective on Mat

    1103:     Input Parameters:
    1104: +   J - the MatSNESMF matrix
    1105: .   fun - the function that checks h
    1106: -   ctx - any context needed by the function

    1108:     Options Database Keys:
    1109: .   -snes_mf_check_positivity

    1111:     Level: advanced

    1113:     Notes: For example, MatSNESMFSetCheckPositivity() insures that all entries
    1114:        of U + h*a are non-negative

    1116: .seealso:  MatSNESMFSetCheckPositivity()
    1117: @*/
    1118: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFSetCheckh(Mat J,PetscErrorCode (*fun)(Vec,Vec,PetscScalar*,void*),void* ctx)
    1119: {
    1120:   PetscErrorCode ierr,(*f)(Mat,PetscErrorCode (*)(Vec,Vec,PetscScalar*,void*),void*);

    1124:   PetscObjectQueryFunction((PetscObject)J,"MatSNESMFSetCheckh_C",(void (**)(void))&f);
    1125:   if (f) {
    1126:     (*f)(J,fun,ctx);
    1127:   }
    1128:   return(0);
    1129: }

    1133: /*@
    1134:     MatSNESMFCheckPositivity - Checks that all entries in U + h*a are positive or
    1135:         zero, decreases h until this is satisfied.

    1137:     Collective on Vec

    1139:     Input Parameters:
    1140: +   U - base vector that is added to
    1141: .   a - vector that is added
    1142: .   h - scaling factor on a
    1143: -   dummy - context variable (unused)

    1145:     Options Database Keys:
    1146: .   -snes_mf_check_positivity

    1148:     Level: advanced

    1150:     Notes: This is rarely used directly, rather it is passed as an argument to
    1151:            MatSNESMFSetCheckh()

    1153: .seealso:  MatSNESMFSetCheckh()
    1154: @*/
    1155: PetscErrorCode PETSCSNES_DLLEXPORT MatSNESMFCheckPositivity(Vec U,Vec a,PetscScalar *h,void *dummy)
    1156: {
    1157:   PetscReal      val, minval;
    1158:   PetscScalar    *u_vec, *a_vec;
    1160:   PetscInt       i,n;
    1161:   MPI_Comm       comm;

    1164:   PetscObjectGetComm((PetscObject)U,&comm);
    1165:   VecGetArray(U,&u_vec);
    1166:   VecGetArray(a,&a_vec);
    1167:   VecGetLocalSize(U,&n);
    1168:   minval = PetscAbsScalar(*h*1.01);
    1169:   for(i=0;i<n;i++) {
    1170:     if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
    1171:       val = PetscAbsScalar(u_vec[i]/a_vec[i]);
    1172:       if (val < minval) minval = val;
    1173:     }
    1174:   }
    1175:   VecRestoreArray(U,&u_vec);
    1176:   VecRestoreArray(a,&a_vec);
    1177:   PetscGlobalMin(&minval,&val,comm);
    1178:   if (val <= PetscAbsScalar(*h)) {
    1179:     PetscLogInfo((U,"MatSNESMFCheckPositivity: Scaling back h from %g to %g\n",PetscRealPart(*h),.99*val));
    1180:     if (PetscRealPart(*h) > 0.0) *h =  0.99*val;
    1181:     else                         *h = -0.99*val;
    1182:   }
    1183:   return(0);
    1184: }