Actual source code: mffd.c

      1: #include <../src/mat/impls/shell/shell.h>
      2: #include <../src/mat/impls/mffd/mffdimpl.h>

      4: PetscFunctionList MatMFFDList              = NULL;
      5: PetscBool         MatMFFDRegisterAllCalled = PETSC_FALSE;

      7: PetscClassId  MATMFFD_CLASSID;
      8: PetscLogEvent MATMFFD_Mult;

     10: static PetscBool MatMFFDPackageInitialized = PETSC_FALSE;

     12: /*@C
     13:   MatMFFDFinalizePackage - This function destroys everything in the MATMFFD` package. It is
     14:   called from `PetscFinalize()`.

     16:   Level: developer

     18: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `PetscFinalize()`, `MatCreateMFFD()`, `MatCreateSNESMF()`
     19: @*/
     20: PetscErrorCode MatMFFDFinalizePackage(void)
     21: {
     22:   PetscFunctionBegin;
     23:   PetscCall(PetscFunctionListDestroy(&MatMFFDList));
     24:   MatMFFDPackageInitialized = PETSC_FALSE;
     25:   MatMFFDRegisterAllCalled  = PETSC_FALSE;
     26:   PetscFunctionReturn(PETSC_SUCCESS);
     27: }

     29: /*@C
     30:   MatMFFDInitializePackage - This function initializes everything in the MATMFFD` package. It is called
     31:   from `MatInitializePackage()`.

     33:   Level: developer

     35: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `PetscInitialize()`
     36: @*/
     37: PetscErrorCode MatMFFDInitializePackage(void)
     38: {
     39:   char      logList[256];
     40:   PetscBool opt, pkg;

     42:   PetscFunctionBegin;
     43:   if (MatMFFDPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
     44:   MatMFFDPackageInitialized = PETSC_TRUE;
     45:   /* Register Classes */
     46:   PetscCall(PetscClassIdRegister("MatMFFD", &MATMFFD_CLASSID));
     47:   /* Register Constructors */
     48:   PetscCall(MatMFFDRegisterAll());
     49:   /* Register Events */
     50:   PetscCall(PetscLogEventRegister("MatMult MF", MATMFFD_CLASSID, &MATMFFD_Mult));
     51:   /* Process Info */
     52:   {
     53:     PetscClassId classids[1];

     55:     classids[0] = MATMFFD_CLASSID;
     56:     PetscCall(PetscInfoProcessClass("matmffd", 1, classids));
     57:   }
     58:   /* Process summary exclusions */
     59:   PetscCall(PetscOptionsGetString(NULL, NULL, "-log_exclude", logList, sizeof(logList), &opt));
     60:   if (opt) {
     61:     PetscCall(PetscStrInList("matmffd", logList, ',', &pkg));
     62:     if (pkg) PetscCall(PetscLogEventExcludeClass(MATMFFD_CLASSID));
     63:   }
     64:   /* Register package finalizer */
     65:   PetscCall(PetscRegisterFinalize(MatMFFDFinalizePackage));
     66:   PetscFunctionReturn(PETSC_SUCCESS);
     67: }

     69: static PetscErrorCode MatMFFDSetType_MFFD(Mat mat, MatMFFDType ftype)
     70: {
     71:   MatMFFD   ctx;
     72:   PetscBool match;
     73:   PetscErrorCode (*r)(MatMFFD);

     75:   PetscFunctionBegin;
     77:   PetscAssertPointer(ftype, 2);
     78:   PetscCall(MatShellGetContext(mat, &ctx));

     80:   /* already set, so just return */
     81:   PetscCall(PetscObjectTypeCompare((PetscObject)ctx, ftype, &match));
     82:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

     84:   /* destroy the old one if it exists */
     85:   PetscTryTypeMethod(ctx, destroy);

     87:   PetscCall(PetscFunctionListFind(MatMFFDList, ftype, &r));
     88:   PetscCheck(r, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown MatMFFD type %s given", ftype);
     89:   PetscCall((*r)(ctx));
     90:   PetscCall(PetscObjectChangeTypeName((PetscObject)ctx, ftype));
     91:   PetscFunctionReturn(PETSC_SUCCESS);
     92: }

     94: /*@
     95:   MatMFFDSetType - Sets the method that is used to compute the
     96:   differencing parameter for finite difference matrix-free formulations.

     98:   Input Parameters:
     99: + mat   - the "matrix-free" matrix created via `MatCreateSNESMF()`, or `MatCreateMFFD()`
    100:           or `MatSetType`(mat,`MATMFFD`);
    101: - ftype - the type requested, either `MATMFFD_WP` or `MATMFFD_DS`

    103:   Level: advanced

    105:   Note:
    106:   For example, such routines can compute `h` for use in
    107:   Jacobian-vector products of the form
    108: .vb

    110:                         F(x+ha) - F(x)
    111:           F'(u)a  ~=  ----------------
    112:                               h
    113: .ve

    115: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MATMFFD_WP`, `MATMFFD_DS`, `MatCreateSNESMF()`, `MatMFFDRegister()`, `MatMFFDSetFunction()`, `MatCreateMFFD()`
    116: @*/
    117: PetscErrorCode MatMFFDSetType(Mat mat, MatMFFDType ftype)
    118: {
    119:   PetscFunctionBegin;
    121:   PetscAssertPointer(ftype, 2);
    122:   PetscTryMethod(mat, "MatMFFDSetType_C", (Mat, MatMFFDType), (mat, ftype));
    123:   PetscFunctionReturn(PETSC_SUCCESS);
    124: }

    126: static PetscErrorCode MatGetDiagonal_MFFD(Mat, Vec);

    128: typedef PetscErrorCode (*FCN1)(void *, Vec); /* force argument to next function to not be extern C*/
    129: static PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat, FCN1 func)
    130: {
    131:   MatMFFD ctx;

    133:   PetscFunctionBegin;
    134:   PetscCall(MatShellGetContext(mat, &ctx));
    135:   ctx->funcisetbase = func;
    136:   PetscFunctionReturn(PETSC_SUCCESS);
    137: }

    139: typedef PetscErrorCode (*FCN2)(void *, PetscInt, Vec, PetscScalar *); /* force argument to next function to not be extern C*/
    140: static PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat, FCN2 funci)
    141: {
    142:   MatMFFD ctx;

    144:   PetscFunctionBegin;
    145:   PetscCall(MatShellGetContext(mat, &ctx));
    146:   ctx->funci = funci;
    147:   PetscCall(MatShellSetOperation(mat, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_MFFD));
    148:   PetscFunctionReturn(PETSC_SUCCESS);
    149: }

    151: static PetscErrorCode MatMFFDGetH_MFFD(Mat mat, PetscScalar *h)
    152: {
    153:   MatMFFD ctx;

    155:   PetscFunctionBegin;
    156:   PetscCall(MatShellGetContext(mat, &ctx));
    157:   *h = ctx->currenth;
    158:   PetscFunctionReturn(PETSC_SUCCESS);
    159: }

    161: static PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
    162: {
    163:   MatMFFD ctx;

    165:   PetscFunctionBegin;
    166:   PetscCall(MatShellGetContext(J, &ctx));
    167:   ctx->ncurrenth = 0;
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: /*@C
    172:   MatMFFDRegister - Adds a method to the `MATMFFD` registry.

    174:   Not Collective, No Fortran Support

    176:   Input Parameters:
    177: + sname    - name of a new user-defined compute-h module
    178: - function - routine to create method context

    180:   Level: developer

    182:   Note:
    183:   `MatMFFDRegister()` may be called multiple times to add several user-defined solvers.

    185:   Example Usage:
    186: .vb
    187:    MatMFFDRegister("my_h", MyHCreate);
    188: .ve

    190:   Then, your solver can be chosen with the procedural interface via `MatMFFDSetType`(mfctx, "my_h")` or at runtime via the option
    191:   `-mat_mffd_type my_h`

    193: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDRegisterAll()`, `MatMFFDRegisterDestroy()`
    194:  @*/
    195: PetscErrorCode MatMFFDRegister(const char sname[], PetscErrorCode (*function)(MatMFFD))
    196: {
    197:   PetscFunctionBegin;
    198:   PetscCall(MatInitializePackage());
    199:   PetscCall(PetscFunctionListAdd(&MatMFFDList, sname, function));
    200:   PetscFunctionReturn(PETSC_SUCCESS);
    201: }

    203: static PetscErrorCode MatDestroy_MFFD(Mat mat)
    204: {
    205:   MatMFFD ctx;

    207:   PetscFunctionBegin;
    208:   PetscCall(MatShellGetContext(mat, &ctx));
    209:   PetscCall(VecDestroy(&ctx->w));
    210:   PetscCall(VecDestroy(&ctx->current_u));
    211:   if (ctx->current_f_allocated) PetscCall(VecDestroy(&ctx->current_f));
    212:   PetscTryTypeMethod(ctx, destroy);
    213:   PetscCall(PetscHeaderDestroy(&ctx));

    215:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetBase_C", NULL));
    216:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunctioniBase_C", NULL));
    217:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunctioni_C", NULL));
    218:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunction_C", NULL));
    219:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunctionError_C", NULL));
    220:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetCheckh_C", NULL));
    221:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetPeriod_C", NULL));
    222:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDResetHHistory_C", NULL));
    223:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetHHistory_C", NULL));
    224:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetType_C", NULL));
    225:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMFFDGetH_C", NULL));
    226:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSNESMFSetReuseBase_C", NULL));
    227:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSNESMFGetReuseBase_C", NULL));
    228:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContext_C", NULL));
    229:   PetscFunctionReturn(PETSC_SUCCESS);
    230: }

    232: /*
    233:    MatMFFDView_MFFD - Views matrix-free parameters.

    235: */
    236: static PetscErrorCode MatView_MFFD(Mat J, PetscViewer viewer)
    237: {
    238:   MatMFFD     ctx;
    239:   PetscBool   iascii, viewbase, viewfunction;
    240:   const char *prefix;

    242:   PetscFunctionBegin;
    243:   PetscCall(MatShellGetContext(J, &ctx));
    244:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    245:   if (iascii) {
    246:     PetscCall(PetscViewerASCIIPrintf(viewer, "Matrix-free approximation:\n"));
    247:     PetscCall(PetscViewerASCIIPushTab(viewer));
    248:     PetscCall(PetscViewerASCIIPrintf(viewer, "err=%g (relative error in function evaluation)\n", (double)ctx->error_rel));
    249:     if (!((PetscObject)ctx)->type_name) {
    250:       PetscCall(PetscViewerASCIIPrintf(viewer, "The compute h routine has not yet been set\n"));
    251:     } else {
    252:       PetscCall(PetscViewerASCIIPrintf(viewer, "Using %s compute h routine\n", ((PetscObject)ctx)->type_name));
    253:     }
    254: #if defined(PETSC_USE_COMPLEX)
    255:     if (ctx->usecomplex) PetscCall(PetscViewerASCIIPrintf(viewer, "Using Lyness complex number trick to compute the matrix-vector product\n"));
    256: #endif
    257:     PetscTryTypeMethod(ctx, view, viewer);
    258:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)J, &prefix));

    260:     PetscCall(PetscOptionsHasName(((PetscObject)J)->options, prefix, "-mat_mffd_view_base", &viewbase));
    261:     if (viewbase) {
    262:       PetscCall(PetscViewerASCIIPrintf(viewer, "Base:\n"));
    263:       PetscCall(VecView(ctx->current_u, viewer));
    264:     }
    265:     PetscCall(PetscOptionsHasName(((PetscObject)J)->options, prefix, "-mat_mffd_view_function", &viewfunction));
    266:     if (viewfunction) {
    267:       PetscCall(PetscViewerASCIIPrintf(viewer, "Function:\n"));
    268:       PetscCall(VecView(ctx->current_f, viewer));
    269:     }
    270:     PetscCall(PetscViewerASCIIPopTab(viewer));
    271:   }
    272:   PetscFunctionReturn(PETSC_SUCCESS);
    273: }

    275: /*
    276:    MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
    277:    allows the user to indicate the beginning of a new linear solve by calling
    278:    MatAssemblyXXX() on the matrix-free matrix. This then allows the
    279:    MatCreateMFFD_WP() to properly compute ||U|| only the first time
    280:    in the linear solver rather than every time.

    282:    This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
    283:    it must be labeled as PETSC_EXTERN
    284: */
    285: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J, MatAssemblyType mt)
    286: {
    287:   MatMFFD j;

    289:   PetscFunctionBegin;
    290:   PetscCall(MatShellGetContext(J, &j));
    291:   PetscCall(MatMFFDResetHHistory(J));
    292:   PetscFunctionReturn(PETSC_SUCCESS);
    293: }

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

    298:         y ~= (F(u + ha) - F(u))/h,
    299:   where F = nonlinear function, as set by SNESSetFunction()
    300:         u = current iterate
    301:         h = difference interval
    302: */
    303: static PetscErrorCode MatMult_MFFD(Mat mat, Vec a, Vec y)
    304: {
    305:   MatMFFD     ctx;
    306:   PetscScalar h;
    307:   Vec         w, U, F;
    308:   PetscBool   zeroa;

    310:   PetscFunctionBegin;
    311:   PetscCall(MatShellGetContext(mat, &ctx));
    312:   PetscCheck(ctx->current_u, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "MatMFFDSetBase() has not been called, this is often caused by forgetting to call MatAssemblyBegin/End on the first Mat in the SNES compute function");
    313:   /* We log matrix-free matrix-vector products separately, so that we can
    314:      separate the performance monitoring from the cases that use conventional
    315:      storage.  We may eventually modify event logging to associate events
    316:      with particular objects, hence alleviating the more general problem. */
    317:   PetscCall(PetscLogEventBegin(MATMFFD_Mult, a, y, 0, 0));

    319:   w = ctx->w;
    320:   U = ctx->current_u;
    321:   F = ctx->current_f;
    322:   /*
    323:       Compute differencing parameter
    324:   */
    325:   if (!((PetscObject)ctx)->type_name) {
    326:     PetscCall(MatMFFDSetType(mat, MATMFFD_WP));
    327:     PetscCall(MatSetFromOptions(mat));
    328:   }
    329:   PetscUseTypeMethod(ctx, compute, U, a, &h, &zeroa);
    330:   if (zeroa) {
    331:     PetscCall(VecSet(y, 0.0));
    332:     PetscCall(PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0));
    333:     PetscFunctionReturn(PETSC_SUCCESS);
    334:   }

    336:   PetscCheck(!mat->erroriffailure || !PetscIsInfOrNanScalar(h), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Computed Nan differencing parameter h");
    337:   if (ctx->checkh) PetscCall((*ctx->checkh)(ctx->checkhctx, U, a, &h));

    339:   /* keep a record of the current differencing parameter h */
    340:   ctx->currenth = h;
    341: #if defined(PETSC_USE_COMPLEX)
    342:   PetscCall(PetscInfo(mat, "Current differencing parameter: %g + %g i\n", (double)PetscRealPart(h), (double)PetscImaginaryPart(h)));
    343: #else
    344:   PetscCall(PetscInfo(mat, "Current differencing parameter: %15.12e\n", (double)PetscRealPart(h)));
    345: #endif
    346:   if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) ctx->historyh[ctx->ncurrenth] = h;
    347:   ctx->ncurrenth++;

    349: #if defined(PETSC_USE_COMPLEX)
    350:   if (ctx->usecomplex) h = PETSC_i * h;
    351: #endif

    353:   /* w = u + ha */
    354:   PetscCall(VecWAXPY(w, h, a, U));

    356:   /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
    357:   if (ctx->ncurrenth == 1 && ctx->current_f_allocated) PetscCall((*ctx->func)(ctx->funcctx, U, F));
    358:   PetscCall((*ctx->func)(ctx->funcctx, w, y));

    360: #if defined(PETSC_USE_COMPLEX)
    361:   if (ctx->usecomplex) {
    362:     PetscCall(VecImaginaryPart(y));
    363:     h = PetscImaginaryPart(h);
    364:   } else {
    365:     PetscCall(VecAXPY(y, -1.0, F));
    366:   }
    367: #else
    368:   PetscCall(VecAXPY(y, -1.0, F));
    369: #endif
    370:   PetscCall(VecScale(y, 1.0 / h));
    371:   if (mat->nullsp) PetscCall(MatNullSpaceRemove(mat->nullsp, y));

    373:   PetscCall(PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0));
    374:   PetscFunctionReturn(PETSC_SUCCESS);
    375: }

    377: /*
    378:   MatGetDiagonal_MFFD - Gets the diagonal for a matrix-free matrix

    380:         y ~= (F(u + ha) - F(u))/h,
    381:   where F = nonlinear function, as set by SNESSetFunction()
    382:         u = current iterate
    383:         h = difference interval
    384: */
    385: static PetscErrorCode MatGetDiagonal_MFFD(Mat mat, Vec a)
    386: {
    387:   MatMFFD     ctx;
    388:   PetscScalar h, *aa, *ww, v;
    389:   PetscReal   epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
    390:   Vec         w, U;
    391:   PetscInt    i, rstart, rend;

    393:   PetscFunctionBegin;
    394:   PetscCall(MatShellGetContext(mat, &ctx));
    395:   PetscCheck(ctx->func, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Requires calling MatMFFDSetFunction() first");
    396:   PetscCheck(ctx->funci, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Requires calling MatMFFDSetFunctioni() first");
    397:   w = ctx->w;
    398:   U = ctx->current_u;
    399:   PetscCall((*ctx->func)(ctx->funcctx, U, a));
    400:   if (ctx->funcisetbase) PetscCall((*ctx->funcisetbase)(ctx->funcctx, U));
    401:   PetscCall(VecCopy(U, w));

    403:   PetscCall(VecGetOwnershipRange(a, &rstart, &rend));
    404:   PetscCall(VecGetArray(a, &aa));
    405:   for (i = rstart; i < rend; i++) {
    406:     PetscCall(VecGetArray(w, &ww));
    407:     h = ww[i - rstart];
    408:     if (h == 0.0) h = 1.0;
    409:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
    410:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
    411:     h *= epsilon;

    413:     ww[i - rstart] += h;
    414:     PetscCall(VecRestoreArray(w, &ww));
    415:     PetscCall((*ctx->funci)(ctx->funcctx, i, w, &v));
    416:     aa[i - rstart] = (v - aa[i - rstart]) / h;

    418:     PetscCall(VecGetArray(w, &ww));
    419:     ww[i - rstart] -= h;
    420:     PetscCall(VecRestoreArray(w, &ww));
    421:   }
    422:   PetscCall(VecRestoreArray(a, &aa));
    423:   PetscFunctionReturn(PETSC_SUCCESS);
    424: }

    426: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J, Vec U, Vec F)
    427: {
    428:   MatMFFD ctx;

    430:   PetscFunctionBegin;
    431:   PetscCall(MatShellGetContext(J, &ctx));
    432:   PetscCall(MatMFFDResetHHistory(J));
    433:   if (!ctx->current_u) {
    434:     PetscCall(VecDuplicate(U, &ctx->current_u));
    435:     PetscCall(VecLockReadPush(ctx->current_u));
    436:   }
    437:   PetscCall(VecLockReadPop(ctx->current_u));
    438:   PetscCall(VecCopy(U, ctx->current_u));
    439:   PetscCall(VecLockReadPush(ctx->current_u));
    440:   if (F) {
    441:     if (ctx->current_f_allocated) PetscCall(VecDestroy(&ctx->current_f));
    442:     ctx->current_f           = F;
    443:     ctx->current_f_allocated = PETSC_FALSE;
    444:   } else if (!ctx->current_f_allocated) {
    445:     PetscCall(MatCreateVecs(J, NULL, &ctx->current_f));
    446:     ctx->current_f_allocated = PETSC_TRUE;
    447:   }
    448:   if (!ctx->w) PetscCall(VecDuplicate(ctx->current_u, &ctx->w));
    449:   J->assembled = PETSC_TRUE;
    450:   PetscFunctionReturn(PETSC_SUCCESS);
    451: }

    453: typedef PetscErrorCode (*FCN3)(void *, Vec, Vec, PetscScalar *); /* force argument to next function to not be extern C*/

    455: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J, FCN3 fun, void *ectx)
    456: {
    457:   MatMFFD ctx;

    459:   PetscFunctionBegin;
    460:   PetscCall(MatShellGetContext(J, &ctx));
    461:   ctx->checkh    = fun;
    462:   ctx->checkhctx = ectx;
    463:   PetscFunctionReturn(PETSC_SUCCESS);
    464: }

    466: /*@
    467:   MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
    468:   MATMFFD` options in the database.

    470:   Collective

    472:   Input Parameters:
    473: + mat    - the `MATMFFD` context
    474: - prefix - the prefix to prepend to all option names

    476:   Note:
    477:   A hyphen (-) must NOT be given at the beginning of the prefix name.
    478:   The first character of all runtime options is AUTOMATICALLY the hyphen.

    480:   Level: advanced

    482: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatSetFromOptions()`, `MatCreateSNESMF()`, `MatCreateMFFD()`
    483: @*/
    484: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat, const char prefix[])
    485: {
    486:   MatMFFD mfctx;

    488:   PetscFunctionBegin;
    490:   PetscCall(MatShellGetContext(mat, &mfctx));
    492:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mfctx, prefix));
    493:   PetscFunctionReturn(PETSC_SUCCESS);
    494: }

    496: static PetscErrorCode MatSetFromOptions_MFFD(Mat mat, PetscOptionItems *PetscOptionsObject)
    497: {
    498:   MatMFFD   mfctx;
    499:   PetscBool flg;
    500:   char      ftype[256];

    502:   PetscFunctionBegin;
    503:   PetscCall(MatShellGetContext(mat, &mfctx));
    505:   PetscObjectOptionsBegin((PetscObject)mfctx);
    506:   PetscCall(PetscOptionsFList("-mat_mffd_type", "Matrix free type", "MatMFFDSetType", MatMFFDList, ((PetscObject)mfctx)->type_name, ftype, 256, &flg));
    507:   if (flg) PetscCall(MatMFFDSetType(mat, ftype));

    509:   PetscCall(PetscOptionsReal("-mat_mffd_err", "set sqrt relative error in function", "MatMFFDSetFunctionError", mfctx->error_rel, &mfctx->error_rel, NULL));
    510:   PetscCall(PetscOptionsInt("-mat_mffd_period", "how often h is recomputed", "MatMFFDSetPeriod", mfctx->recomputeperiod, &mfctx->recomputeperiod, NULL));

    512:   flg = PETSC_FALSE;
    513:   PetscCall(PetscOptionsBool("-mat_mffd_check_positivity", "Insure that U + h*a is nonnegative", "MatMFFDSetCheckh", flg, &flg, NULL));
    514:   if (flg) PetscCall(MatMFFDSetCheckh(mat, MatMFFDCheckPositivity, NULL));
    515: #if defined(PETSC_USE_COMPLEX)
    516:   PetscCall(PetscOptionsBool("-mat_mffd_complex", "Use Lyness complex number trick to compute the matrix-vector product", "None", mfctx->usecomplex, &mfctx->usecomplex, NULL));
    517: #endif
    518:   PetscTryTypeMethod(mfctx, setfromoptions, PetscOptionsObject);
    519:   PetscOptionsEnd();
    520:   PetscFunctionReturn(PETSC_SUCCESS);
    521: }

    523: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat, PetscInt period)
    524: {
    525:   MatMFFD ctx;

    527:   PetscFunctionBegin;
    528:   PetscCall(MatShellGetContext(mat, &ctx));
    529:   ctx->recomputeperiod = period;
    530:   PetscFunctionReturn(PETSC_SUCCESS);
    531: }

    533: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat, PetscErrorCode (*func)(void *, Vec, Vec), void *funcctx)
    534: {
    535:   MatMFFD ctx;

    537:   PetscFunctionBegin;
    538:   PetscCall(MatShellGetContext(mat, &ctx));
    539:   ctx->func    = func;
    540:   ctx->funcctx = funcctx;
    541:   PetscFunctionReturn(PETSC_SUCCESS);
    542: }

    544: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat, PetscReal error)
    545: {
    546:   PetscFunctionBegin;
    547:   if (error != (PetscReal)PETSC_DEFAULT) {
    548:     MatMFFD ctx;

    550:     PetscCall(MatShellGetContext(mat, &ctx));
    551:     ctx->error_rel = error;
    552:   }
    553:   PetscFunctionReturn(PETSC_SUCCESS);
    554: }

    556: static PetscErrorCode MatMFFDSetHHistory_MFFD(Mat J, PetscScalar history[], PetscInt nhistory)
    557: {
    558:   MatMFFD ctx;

    560:   PetscFunctionBegin;
    561:   PetscCall(MatShellGetContext(J, &ctx));
    562:   ctx->historyh    = history;
    563:   ctx->maxcurrenth = nhistory;
    564:   ctx->currenth    = 0.;
    565:   PetscFunctionReturn(PETSC_SUCCESS);
    566: }

    568: /*MC
    569:   MATMFFD - "mffd" - A matrix-free matrix type.

    571:   Level: advanced

    573:   Developer Notes:
    574:   This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code

    576:   Users should not MatShell... operations on this class, there is some error checking for that incorrect usage

    578: .seealso: [](ch_matrices), `Mat`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MatMFFDSetFunction()`, `MatMFFDSetType()`,
    579:           `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`
    580:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`,
    581:           `MatMFFDGetH()`,
    582: M*/
    583: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
    584: {
    585:   MatMFFD mfctx;

    587:   PetscFunctionBegin;
    588:   PetscCall(MatMFFDInitializePackage());

    590:   PetscCall(PetscHeaderCreate(mfctx, MATMFFD_CLASSID, "MatMFFD", "Matrix-free Finite Differencing", "Mat", PetscObjectComm((PetscObject)A), NULL, NULL));

    592:   mfctx->error_rel                = PETSC_SQRT_MACHINE_EPSILON;
    593:   mfctx->recomputeperiod          = 1;
    594:   mfctx->count                    = 0;
    595:   mfctx->currenth                 = 0.0;
    596:   mfctx->historyh                 = NULL;
    597:   mfctx->ncurrenth                = 0;
    598:   mfctx->maxcurrenth              = 0;
    599:   ((PetscObject)mfctx)->type_name = NULL;

    601:   /*
    602:      Create the empty data structure to contain compute-h routines.
    603:      These will be filled in below from the command line options or
    604:      a later call with MatMFFDSetType() or if that is not called
    605:      then it will default in the first use of MatMult_MFFD()
    606:   */
    607:   mfctx->ops->compute        = NULL;
    608:   mfctx->ops->destroy        = NULL;
    609:   mfctx->ops->view           = NULL;
    610:   mfctx->ops->setfromoptions = NULL;
    611:   mfctx->hctx                = NULL;

    613:   mfctx->func    = NULL;
    614:   mfctx->funcctx = NULL;
    615:   mfctx->w       = NULL;
    616:   mfctx->mat     = A;

    618:   PetscCall(MatSetType(A, MATSHELL));
    619:   PetscCall(MatShellSetContext(A, mfctx));
    620:   PetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))MatMult_MFFD));
    621:   PetscCall(MatShellSetOperation(A, MATOP_DESTROY, (void (*)(void))MatDestroy_MFFD));
    622:   PetscCall(MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MFFD));
    623:   PetscCall(MatShellSetOperation(A, MATOP_ASSEMBLY_END, (void (*)(void))MatAssemblyEnd_MFFD));
    624:   PetscCall(MatShellSetOperation(A, MATOP_SET_FROM_OPTIONS, (void (*)(void))MatSetFromOptions_MFFD));
    625:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Immutable));
    626:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
    627:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));

    629:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetBase_C", MatMFFDSetBase_MFFD));
    630:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctioniBase_C", MatMFFDSetFunctioniBase_MFFD));
    631:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctioni_C", MatMFFDSetFunctioni_MFFD));
    632:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunction_C", MatMFFDSetFunction_MFFD));
    633:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetCheckh_C", MatMFFDSetCheckh_MFFD));
    634:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetPeriod_C", MatMFFDSetPeriod_MFFD));
    635:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctionError_C", MatMFFDSetFunctionError_MFFD));
    636:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDResetHHistory_C", MatMFFDResetHHistory_MFFD));
    637:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetHHistory_C", MatMFFDSetHHistory_MFFD));
    638:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetType_C", MatMFFDSetType_MFFD));
    639:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDGetH_C", MatMFFDGetH_MFFD));
    640:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATMFFD));
    641:   PetscFunctionReturn(PETSC_SUCCESS);
    642: }

    644: /*@
    645:   MatCreateMFFD - Creates a matrix-free matrix of type `MATMFFD` that uses finite differences on a provided function to
    646:   approximately multiply a vector by the matrix (Jacobian) . See also `MatCreateSNESMF()`

    648:   Collective

    650:   Input Parameters:
    651: + comm - MPI communicator
    652: . m    - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
    653:          This value should be the same as the local size used in creating the
    654:          y vector for the matrix-vector product y = Ax.
    655: . n    - This value should be the same as the local size used in creating the
    656:          x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
    657:          calculated if `N` is given) For square matrices `n` is almost always `m`.
    658: . M    - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
    659: - N    - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)

    661:   Output Parameter:
    662: . J - the matrix-free matrix

    664:   Options Database Keys:
    665: + -mat_mffd_type             - wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
    666: . -mat_mffd_err              - square root of estimated relative error in function evaluation
    667: . -mat_mffd_period           - how often h is recomputed, defaults to 1, every time
    668: . -mat_mffd_check_positivity - possibly decrease `h` until U + h*a has only positive values
    669: . -mat_mffd_umin <umin>      - Sets umin (for default PETSc routine that computes h only)
    670: . -mat_mffd_complex          - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
    671:                                (requires real valued functions but that PETSc be configured for complex numbers)
    672: . -snes_mf                   - use the finite difference based matrix-free matrix with `SNESSolve()` and no preconditioner
    673: - -snes_mf_operator          - use the finite difference based matrix-free matrix with `SNESSolve()` but construct a preconditioner
    674:                                using the matrix passed as `pmat` to `SNESSetJacobian()`.

    676:   Level: advanced

    678:   Notes:
    679:   Use `MatMFFDSetFunction()` to provide the function that will be differenced to compute the matrix-vector product.

    681:   The matrix-free matrix context contains the function pointers
    682:   and work space for performing finite difference approximations of
    683:   Jacobian-vector products, F'(u)*a,

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

    687: .vb
    688:      F'(u)*a = [F(u+h*a) - F(u)]/h where
    689:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
    690:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
    691:  where
    692:      error_rel = square root of relative error in function evaluation
    693:      umin = minimum iterate parameter
    694: .ve

    696:   To have `SNES` use the matrix-free finite difference matrix-vector product and not provide a separate matrix
    697:   from which to compute the preconditioner (the `pmat` argument `SNESSetJacobian()`), then simply call `SNESSetJacobian()`
    698:   with `NULL` for the matrices and `MatMFFDComputeJacobian()`. Or use the options database option `-snes_mf`

    700:   The user can set `error_rel` via `MatMFFDSetFunctionError()` and `umin` via `MatMFFDDSSetUmin()`.

    702:   Use `MATSHELL` or `MatCreateShell()` to provide your own custom matrix-vector operation.

    704: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatDestroy()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`
    705:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`, `MatCreateShell()`, `MATSHELL`,
    706:           `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`
    707: @*/
    708: PetscErrorCode MatCreateMFFD(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, Mat *J)
    709: {
    710:   PetscFunctionBegin;
    711:   PetscCall(MatCreate(comm, J));
    712:   PetscCall(MatSetSizes(*J, m, n, M, N));
    713:   PetscCall(MatSetType(*J, MATMFFD));
    714:   PetscCall(MatSetUp(*J));
    715:   PetscFunctionReturn(PETSC_SUCCESS);
    716: }

    718: /*@
    719:   MatMFFDGetH - Gets the last value that was used as the differencing for a `MATMFFD` matrix
    720:   parameter.

    722:   Not Collective

    724:   Input Parameters:
    725: . mat - the `MATMFFD` matrix

    727:   Output Parameter:
    728: . h - the differencing step size

    730:   Level: advanced

    732: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDSetHHistory()`, `MatCreateMFFD()`, `MatMFFDResetHHistory()`
    733: @*/
    734: PetscErrorCode MatMFFDGetH(Mat mat, PetscScalar *h)
    735: {
    736:   PetscFunctionBegin;
    738:   PetscAssertPointer(h, 2);
    739:   PetscUseMethod(mat, "MatMFFDGetH_C", (Mat, PetscScalar *), (mat, h));
    740:   PetscFunctionReturn(PETSC_SUCCESS);
    741: }

    743: /*@C
    744:   MatMFFDSetFunction - Sets the function used in applying the matrix-free `MATMFFD` matrix.

    746:   Logically Collective

    748:   Input Parameters:
    749: + mat     - the matrix-free matrix `MATMFFD` created via `MatCreateSNESMF()` or `MatCreateMFFD()`
    750: . func    - the function to use
    751: - funcctx - optional function context passed to function

    753:   Calling sequence of `func`:
    754: + funcctx - user provided context
    755: . x       - input vector
    756: - f       - computed output function

    758:   Level: advanced

    760:   Notes:
    761:   If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix-free
    762:   matrix inside your compute Jacobian routine

    764:   If this is not set then it will use the function set with `SNESSetFunction()` if `MatCreateSNESMF()` was used.

    766: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`,
    767:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESSetFunction()`
    768: @*/
    769: PetscErrorCode MatMFFDSetFunction(Mat mat, PetscErrorCode (*func)(void *funcctx, Vec x, Vec f), void *funcctx)
    770: {
    771:   PetscFunctionBegin;
    773:   PetscTryMethod(mat, "MatMFFDSetFunction_C", (Mat, PetscErrorCode (*)(void *, Vec, Vec), void *), (mat, func, funcctx));
    774:   PetscFunctionReturn(PETSC_SUCCESS);
    775: }

    777: /*@C
    778:   MatMFFDSetFunctioni - Sets the function for a single component for a `MATMFFD` matrix

    780:   Logically Collective

    782:   Input Parameters:
    783: + mat   - the matrix-free matrix `MATMFFD`
    784: - funci - the function to use

    786:   Level: advanced

    788:   Notes:
    789:   If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix-free
    790:   matrix inside your compute Jacobian routine.

    792:   This function is necessary to compute the diagonal of the matrix.
    793:   funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.

    795: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESSetFunction()`, `MatGetDiagonal()`
    796: @*/
    797: PetscErrorCode MatMFFDSetFunctioni(Mat mat, PetscErrorCode (*funci)(void *, PetscInt, Vec, PetscScalar *))
    798: {
    799:   PetscFunctionBegin;
    801:   PetscTryMethod(mat, "MatMFFDSetFunctioni_C", (Mat, PetscErrorCode (*)(void *, PetscInt, Vec, PetscScalar *)), (mat, funci));
    802:   PetscFunctionReturn(PETSC_SUCCESS);
    803: }

    805: /*@C
    806:   MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation for a `MATMFFD` matrix

    808:   Logically Collective

    810:   Input Parameters:
    811: + mat  - the `MATMFFD` matrix-free matrix
    812: - func - the function to use

    814:   Level: advanced

    816:   Notes:
    817:   If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix-free
    818:   matrix inside your compute Jacobian routine.

    820:   This function is necessary to compute the diagonal of the matrix, used for example with `PCJACOBI`

    822: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`,
    823:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESSetFunction()`, `MatGetDiagonal()`
    824: @*/
    825: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat, PetscErrorCode (*func)(void *, Vec))
    826: {
    827:   PetscFunctionBegin;
    829:   PetscTryMethod(mat, "MatMFFDSetFunctioniBase_C", (Mat, PetscErrorCode (*)(void *, Vec)), (mat, func));
    830:   PetscFunctionReturn(PETSC_SUCCESS);
    831: }

    833: /*@
    834:   MatMFFDSetPeriod - Sets how often h is recomputed for a `MATMFFD` matrix, by default it is every time

    836:   Logically Collective

    838:   Input Parameters:
    839: + mat    - the `MATMFFD` matrix-free matrix
    840: - period - 1 for every time, 2 for every second etc

    842:   Options Database Key:
    843: . -mat_mffd_period <period> - Sets how often `h` is recomputed

    845:   Level: advanced

    847: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`,
    848:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
    849: @*/
    850: PetscErrorCode MatMFFDSetPeriod(Mat mat, PetscInt period)
    851: {
    852:   PetscFunctionBegin;
    855:   PetscTryMethod(mat, "MatMFFDSetPeriod_C", (Mat, PetscInt), (mat, period));
    856:   PetscFunctionReturn(PETSC_SUCCESS);
    857: }

    859: /*@
    860:   MatMFFDSetFunctionError - Sets the error_rel for the approximation of matrix-vector products using finite differences with the `MATMFFD` matrix

    862:   Logically Collective

    864:   Input Parameters:
    865: + mat   - the `MATMFFD` matrix-free matrix
    866: - error - relative error (should be set to the square root of the relative error in the function evaluations)

    868:   Options Database Key:
    869: . -mat_mffd_err <error_rel> - Sets error_rel

    871:   Level: advanced

    873:   Note:
    874:   The default matrix-free matrix-vector product routine computes
    875: .vb
    876:      F'(u)*a = [F(u+h*a) - F(u)]/h where
    877:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
    878:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
    879: .ve

    881: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`,
    882:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
    883: @*/
    884: PetscErrorCode MatMFFDSetFunctionError(Mat mat, PetscReal error)
    885: {
    886:   PetscFunctionBegin;
    889:   PetscTryMethod(mat, "MatMFFDSetFunctionError_C", (Mat, PetscReal), (mat, error));
    890:   PetscFunctionReturn(PETSC_SUCCESS);
    891: }

    893: /*@
    894:   MatMFFDSetHHistory - Sets an array to collect a history of the
    895:   differencing values (h) computed for the matrix-free product `MATMFFD` matrix

    897:   Logically Collective

    899:   Input Parameters:
    900: + J        - the `MATMFFD` matrix-free matrix
    901: . history  - space to hold the history
    902: - nhistory - number of entries in history, if more entries are generated than
    903:               nhistory, then the later ones are discarded

    905:   Level: advanced

    907:   Note:
    908:   Use `MatMFFDResetHHistory()` to reset the history counter and collect
    909:   a new batch of differencing parameters, h.

    911: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
    912:           `MatMFFDResetHHistory()`, `MatMFFDSetFunctionError()`
    913: @*/
    914: PetscErrorCode MatMFFDSetHHistory(Mat J, PetscScalar history[], PetscInt nhistory)
    915: {
    916:   PetscFunctionBegin;
    918:   if (history) PetscAssertPointer(history, 2);
    920:   PetscUseMethod(J, "MatMFFDSetHHistory_C", (Mat, PetscScalar[], PetscInt), (J, history, nhistory));
    921:   PetscFunctionReturn(PETSC_SUCCESS);
    922: }

    924: /*@
    925:   MatMFFDResetHHistory - Resets the counter to zero to begin
    926:   collecting a new set of differencing histories for the `MATMFFD` matrix

    928:   Logically Collective

    930:   Input Parameter:
    931: . J - the matrix-free matrix context

    933:   Level: advanced

    935:   Note:
    936:   Use `MatMFFDSetHHistory()` to create the original history counter.

    938: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
    939:           `MatMFFDSetHHistory()`, `MatMFFDSetFunctionError()`
    940: @*/
    941: PetscErrorCode MatMFFDResetHHistory(Mat J)
    942: {
    943:   PetscFunctionBegin;
    945:   PetscTryMethod(J, "MatMFFDResetHHistory_C", (Mat), (J));
    946:   PetscFunctionReturn(PETSC_SUCCESS);
    947: }

    949: /*@
    950:   MatMFFDSetBase - Sets the vector `U` at which matrix vector products of the
    951:   Jacobian are computed for the `MATMFFD` matrix

    953:   Logically Collective

    955:   Input Parameters:
    956: + J - the `MATMFFD` matrix
    957: . U - the vector
    958: - F - (optional) vector that contains F(u) if it has been already computed

    960:   Level: advanced

    962:   Notes:
    963:   This is rarely used directly

    965:   If `F` is provided then it is not recomputed. Otherwise the function is evaluated at the base
    966:   point during the first `MatMult()` after each call to `MatMFFDSetBase()`.

    968: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMult()`
    969: @*/
    970: PetscErrorCode MatMFFDSetBase(Mat J, Vec U, Vec F)
    971: {
    972:   PetscFunctionBegin;
    976:   PetscTryMethod(J, "MatMFFDSetBase_C", (Mat, Vec, Vec), (J, U, F));
    977:   PetscFunctionReturn(PETSC_SUCCESS);
    978: }

    980: /*@C
    981:   MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
    982:   it to satisfy some criteria for the `MATMFFD` matrix

    984:   Logically Collective

    986:   Input Parameters:
    987: + J   - the `MATMFFD` matrix
    988: . fun - the function that checks `h`
    989: - ctx - any context needed by the function

    991:   Options Database Keys:
    992: . -mat_mffd_check_positivity <bool> - Insure that U + h*a is non-negative

    994:   Level: advanced

    996:   Notes:
    997:   For example, `MatMFFDCheckPositivity()` insures that all entries of U + h*a are non-negative

    999:   The function you provide is called after the default `h` has been computed and allows you to
    1000:   modify it.

    1002: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDCheckPositivity()`
    1003: @*/
    1004: PetscErrorCode MatMFFDSetCheckh(Mat J, PetscErrorCode (*fun)(void *, Vec, Vec, PetscScalar *), void *ctx)
    1005: {
    1006:   PetscFunctionBegin;
    1008:   PetscTryMethod(J, "MatMFFDSetCheckh_C", (Mat, PetscErrorCode (*)(void *, Vec, Vec, PetscScalar *), void *), (J, fun, ctx));
    1009:   PetscFunctionReturn(PETSC_SUCCESS);
    1010: }

    1012: /*@
    1013:   MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
    1014:   zero, decreases h until this is satisfied for a `MATMFFD` matrix

    1016:   Logically Collective

    1018:   Input Parameters:
    1019: + U     - base vector that is added to
    1020: . a     - vector that is added
    1021: . h     - scaling factor on a
    1022: - dummy - context variable (unused)

    1024:   Options Database Keys:
    1025: . -mat_mffd_check_positivity <bool> - Insure that U + h*a is nonnegative

    1027:   Level: advanced

    1029:   Note:
    1030:   This is rarely used directly, rather it is passed as an argument to `MatMFFDSetCheckh()`

    1032: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDSetCheckh()`
    1033: @*/
    1034: PetscErrorCode MatMFFDCheckPositivity(void *dummy, Vec U, Vec a, PetscScalar *h)
    1035: {
    1036:   PetscReal    val, minval;
    1037:   PetscScalar *u_vec, *a_vec;
    1038:   PetscInt     i, n;
    1039:   MPI_Comm     comm;

    1041:   PetscFunctionBegin;
    1044:   PetscAssertPointer(h, 4);
    1045:   PetscCall(PetscObjectGetComm((PetscObject)U, &comm));
    1046:   PetscCall(VecGetArray(U, &u_vec));
    1047:   PetscCall(VecGetArray(a, &a_vec));
    1048:   PetscCall(VecGetLocalSize(U, &n));
    1049:   minval = PetscAbsScalar(*h) * PetscRealConstant(1.01);
    1050:   for (i = 0; i < n; i++) {
    1051:     if (PetscRealPart(u_vec[i] + *h * a_vec[i]) <= 0.0) {
    1052:       val = PetscAbsScalar(u_vec[i] / a_vec[i]);
    1053:       if (val < minval) minval = val;
    1054:     }
    1055:   }
    1056:   PetscCall(VecRestoreArray(U, &u_vec));
    1057:   PetscCall(VecRestoreArray(a, &a_vec));
    1058:   PetscCallMPI(MPIU_Allreduce(&minval, &val, 1, MPIU_REAL, MPIU_MIN, comm));
    1059:   if (val <= PetscAbsScalar(*h)) {
    1060:     PetscCall(PetscInfo(U, "Scaling back h from %g to %g\n", (double)PetscRealPart(*h), (double)(.99 * val)));
    1061:     if (PetscRealPart(*h) > 0.0) *h = 0.99 * val;
    1062:     else *h = -0.99 * val;
    1063:   }
    1064:   PetscFunctionReturn(PETSC_SUCCESS);
    1065: }