Actual source code: mffd.c

  1: #include <petsc/private/matimpl.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: /*@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;
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

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
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: static 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: + mat    - 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: static 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()`, `MatMFFDResetHHistory()`
724: @*/
725: PetscErrorCode MatMFFDGetH(Mat mat, PetscScalar *h)
726: {
727:   PetscFunctionBegin;
729:   PetscAssertPointer(h, 2);
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: + funcctx - user provided context
746: . x       - input vector
747: - f       - computed output function

749:   Level: advanced

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

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

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

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

771:   Logically Collective

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

777:   Level: advanced

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

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

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

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

799:   Logically Collective

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

805:   Level: advanced

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

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

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

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

827:   Logically Collective

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

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

836:   Level: advanced

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

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

853:   Logically Collective

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

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

862:   Level: advanced

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

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

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

888:   Logically Collective

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

896:   Level: advanced

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

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

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

919:   Logically Collective

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

924:   Level: advanced

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

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

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

944:   Logically Collective

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

951:   Level: advanced

953:   Notes:
954:   This is rarely used directly

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

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

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

975:   Logically Collective

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

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

985:   Level: advanced

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

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

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

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

1007:   Logically Collective

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

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

1018:   Level: advanced

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

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

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