Actual source code: mffd.c


  2: #include <petsc/private/matimpl.h>
  3: #include <../src/mat/impls/mffd/mffdimpl.h>

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

  8: PetscClassId  MATMFFD_CLASSID;
  9: PetscLogEvent MATMFFD_Mult;

 11: 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;
 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, PETSC_COMM_SELF, 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: /*@C
 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;
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

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:    Sample usage:
186: .vb
187:    MatMFFDRegister("my_h", MyHCreate);
188: .ve

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

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

205: static PetscErrorCode MatDestroy_MFFD(Mat mat)
206: {
207:   MatMFFD ctx;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

471:    Collective

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

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

481:    Level: advanced

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

489:   PetscFunctionBegin;
491:   PetscCall(MatShellGetContext(mat, &mfctx));
493:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mfctx, prefix));
494:   PetscFunctionReturn(PETSC_SUCCESS);
495: }

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

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

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

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

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

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

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

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

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

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

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

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

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

572:   Level: advanced

574:   Developers Note:
575:   This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code

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

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

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

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

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

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

617:   PetscCall(MatSetType(A, MATSHELL));
618:   PetscCall(MatShellSetContext(A, mfctx));
619:   PetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))MatMult_MFFD));
620:   PetscCall(MatShellSetOperation(A, MATOP_DESTROY, (void (*)(void))MatDestroy_MFFD));
621:   PetscCall(MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MFFD));
622:   PetscCall(MatShellSetOperation(A, MATOP_ASSEMBLY_END, (void (*)(void))MatAssemblyEnd_MFFD));
623:   PetscCall(MatShellSetOperation(A, MATOP_SET_FROM_OPTIONS, (void (*)(void))MatSetFromOptions_MFFD));

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

640: /*@
641:    MatCreateMFFD - Creates a matrix-free matrix of type `MATMFFD`. See also `MatCreateSNESMF()`

643:    Collective

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

656:    Output Parameter:
657: .  J - the matrix-free matrix

659:    Options Database Keys:
660: +  -mat_mffd_type - wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
661: .  -mat_mffd_err - square root of estimated relative error in function evaluation
662: .  -mat_mffd_period - how often h is recomputed, defaults to 1, every time
663: .  -mat_mffd_check_positivity - possibly decrease `h` until U + h*a has only positive values
664: .  -mat_mffd_umin <umin> - Sets umin (for default PETSc routine that computes h only)
665: -  -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
666:                        (requires real valued functions but that PETSc be configured for complex numbers)

668:    Level: advanced

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

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

677: .vb
678:      F'(u)*a = [F(u+h*a) - F(u)]/h where
679:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
680:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
681:  where
682:      error_rel = square root of relative error in function evaluation
683:      umin = minimum iterate parameter
684: .ve

686:    You can call `SNESSetJacobian()` with `MatMFFDComputeJacobian()` if you are using matrix and not a different
687:    preconditioner matrix

689:    The user can set the error_rel via `MatMFFDSetFunctionError()` and
690:    umin via `MatMFFDDSSetUmin()`.

692:    The user should call `MatDestroy()` when finished with the matrix-free
693:    matrix context.

695: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatDestroy()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`
696:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`,
697:           `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`
698: @*/
699: PetscErrorCode MatCreateMFFD(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, Mat *J)
700: {
701:   PetscFunctionBegin;
702:   PetscCall(MatCreate(comm, J));
703:   PetscCall(MatSetSizes(*J, m, n, M, N));
704:   PetscCall(MatSetType(*J, MATMFFD));
705:   PetscCall(MatSetUp(*J));
706:   PetscFunctionReturn(PETSC_SUCCESS);
707: }

709: /*@
710:    MatMFFDGetH - Gets the last value that was used as the differencing for a `MATMFFD` matrix
711:    parameter.

713:    Not Collective

715:    Input Parameters:
716: .  mat - the `MATMFFD` matrix

718:    Output Parameter:
719: .  h - the differencing step size

721:    Level: advanced

723: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDSetHHistory()`, `MatCreateMFFD()`, `MATMFFD`, `MatMFFDResetHHistory()`
724: @*/
725: PetscErrorCode MatMFFDGetH(Mat mat, PetscScalar *h)
726: {
727:   PetscFunctionBegin;
730:   PetscUseMethod(mat, "MatMFFDGetH_C", (Mat, PetscScalar *), (mat, h));
731:   PetscFunctionReturn(PETSC_SUCCESS);
732: }

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

737:    Logically Collective

739:    Input Parameters:
740: +  mat - the matrix free matrix `MATMFFD` created via `MatCreateSNESMF()` or `MatCreateMFFD()`
741: .  func - the function to use
742: -  funcctx - optional function context passed to function

744:    Calling Sequence of `func`:
745: $  PetscErrorCode func(void *funcctx, Vec x, Vec f)
746: +  funcctx - user provided context
747: .  x - input vector
748: -  f - computed output function

750:    Level: advanced

752:    Notes:
753:     If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix free
754:     matrix inside your compute Jacobian routine

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

758: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`, `MATMFFD`,
759:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`
760: @*/
761: PetscErrorCode MatMFFDSetFunction(Mat mat, PetscErrorCode (*func)(void *, Vec, Vec), void *funcctx)
762: {
763:   PetscFunctionBegin;
765:   PetscTryMethod(mat, "MatMFFDSetFunction_C", (Mat, PetscErrorCode(*)(void *, Vec, Vec), void *), (mat, func, funcctx));
766:   PetscFunctionReturn(PETSC_SUCCESS);
767: }

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

772:    Logically Collective

774:    Input Parameters:
775: +  mat - the matrix free matrix `MATMFFD`
776: -  funci - the function to use

778:    Level: advanced

780:    Notes:
781:     If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix free
782:     matrix inside your compute Jacobian routine.

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

787: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`, `MatGetDiagonal()`
788: @*/
789: PetscErrorCode MatMFFDSetFunctioni(Mat mat, PetscErrorCode (*funci)(void *, PetscInt, Vec, PetscScalar *))
790: {
791:   PetscFunctionBegin;
793:   PetscTryMethod(mat, "MatMFFDSetFunctioni_C", (Mat, PetscErrorCode(*)(void *, PetscInt, Vec, PetscScalar *)), (mat, funci));
794:   PetscFunctionReturn(PETSC_SUCCESS);
795: }

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

800:    Logically Collective

802:    Input Parameters:
803: +  mat - the `MATMFFD` matrix free matrix
804: -  func - the function to use

806:    Level: advanced

808:    Notes:
809:     If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix free
810:     matrix inside your compute Jacobian routine.

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

814: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`, `MATMFFD`
815:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`, `MatGetDiagonal()`
816: @*/
817: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat, PetscErrorCode (*func)(void *, Vec))
818: {
819:   PetscFunctionBegin;
821:   PetscTryMethod(mat, "MatMFFDSetFunctioniBase_C", (Mat, PetscErrorCode(*)(void *, Vec)), (mat, func));
822:   PetscFunctionReturn(PETSC_SUCCESS);
823: }

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

828:    Logically Collective

830:    Input Parameters:
831: +  mat - the `MATMFFD` matrix free matrix
832: -  period - 1 for every time, 2 for every second etc

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

837:    Level: advanced

839: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`,
840:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
841: @*/
842: PetscErrorCode MatMFFDSetPeriod(Mat mat, PetscInt period)
843: {
844:   PetscFunctionBegin;
847:   PetscTryMethod(mat, "MatMFFDSetPeriod_C", (Mat, PetscInt), (mat, period));
848:   PetscFunctionReturn(PETSC_SUCCESS);
849: }

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

854:    Logically Collective

856:    Input Parameters:
857: +  mat - the `MATMFFD` matrix free matrix
858: -  error_rel - relative error (should be set to the square root of the relative error in the function evaluations)

860:    Options Database Key:
861: .  -mat_mffd_err <error_rel> - Sets error_rel

863:    Level: advanced

865:    Note:
866:    The default matrix-free matrix-vector product routine computes
867: .vb
868:      F'(u)*a = [F(u+h*a) - F(u)]/h where
869:      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
870:        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
871: .ve

873: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`, `MATMFFD`
874:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
875: @*/
876: PetscErrorCode MatMFFDSetFunctionError(Mat mat, PetscReal error)
877: {
878:   PetscFunctionBegin;
881:   PetscTryMethod(mat, "MatMFFDSetFunctionError_C", (Mat, PetscReal), (mat, error));
882:   PetscFunctionReturn(PETSC_SUCCESS);
883: }

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

889:    Logically Collective

891:    Input Parameters:
892: +  J - the `MATMFFD` matrix-free matrix
893: .  history - space to hold the history
894: -  nhistory - number of entries in history, if more entries are generated than
895:               nhistory, then the later ones are discarded

897:    Level: advanced

899:    Note:
900:    Use `MatMFFDResetHHistory()` to reset the history counter and collect
901:    a new batch of differencing parameters, h.

903: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
904:           `MatMFFDResetHHistory()`, `MatMFFDSetFunctionError()`
905: @*/
906: PetscErrorCode MatMFFDSetHHistory(Mat J, PetscScalar history[], PetscInt nhistory)
907: {
908:   PetscFunctionBegin;
912:   PetscUseMethod(J, "MatMFFDSetHHistory_C", (Mat, PetscScalar[], PetscInt), (J, history, nhistory));
913:   PetscFunctionReturn(PETSC_SUCCESS);
914: }

916: /*@
917:    MatMFFDResetHHistory - Resets the counter to zero to begin
918:    collecting a new set of differencing histories for the `MATMFFD` matrix

920:    Logically Collective

922:    Input Parameter:
923: .  J - the matrix-free matrix context

925:    Level: advanced

927:    Note:
928:    Use `MatMFFDSetHHistory()` to create the original history counter.

930: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
931:           `MatMFFDSetHHistory()`, `MatMFFDSetFunctionError()`
932: @*/
933: PetscErrorCode MatMFFDResetHHistory(Mat J)
934: {
935:   PetscFunctionBegin;
937:   PetscTryMethod(J, "MatMFFDResetHHistory_C", (Mat), (J));
938:   PetscFunctionReturn(PETSC_SUCCESS);
939: }

941: /*@
942:     MatMFFDSetBase - Sets the vector `U` at which matrix vector products of the
943:         Jacobian are computed for the `MATMFFD` matrix

945:     Logically Collective

947:     Input Parameters:
948: +   J - the `MATMFFD` matrix
949: .   U - the vector
950: -   F - (optional) vector that contains F(u) if it has been already computed

952:     Level: advanced

954:     Notes:
955:     This is rarely used directly

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

960: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMult()`, `MatMFFDSetBase()`
961: @*/
962: PetscErrorCode MatMFFDSetBase(Mat J, Vec U, Vec F)
963: {
964:   PetscFunctionBegin;
968:   PetscTryMethod(J, "MatMFFDSetBase_C", (Mat, Vec, Vec), (J, U, F));
969:   PetscFunctionReturn(PETSC_SUCCESS);
970: }

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

976:     Logically Collective

978:     Input Parameters:
979: +   J - the `MATMFFD` matrix
980: .   fun - the function that checks `h`
981: -   ctx - any context needed by the function

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

986:     Level: advanced

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

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

994: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDCheckPositivity()`
995: @*/
996: PetscErrorCode MatMFFDSetCheckh(Mat J, PetscErrorCode (*fun)(void *, Vec, Vec, PetscScalar *), void *ctx)
997: {
998:   PetscFunctionBegin;
1000:   PetscTryMethod(J, "MatMFFDSetCheckh_C", (Mat, PetscErrorCode(*)(void *, Vec, Vec, PetscScalar *), void *), (J, fun, ctx));
1001:   PetscFunctionReturn(PETSC_SUCCESS);
1002: }

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

1008:     Logically Collective

1010:     Input Parameters:
1011: +   U - base vector that is added to
1012: .   a - vector that is added
1013: .   h - scaling factor on a
1014: -   dummy - context variable (unused)

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

1019:     Level: advanced

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

1024: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDSetCheckh()`
1025: @*/
1026: PetscErrorCode MatMFFDCheckPositivity(void *dummy, Vec U, Vec a, PetscScalar *h)
1027: {
1028:   PetscReal    val, minval;
1029:   PetscScalar *u_vec, *a_vec;
1030:   PetscInt     i, n;
1031:   MPI_Comm     comm;

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