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, (PetscErrorCodeFn *)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 isascii, viewbase, viewfunction;
240: const char *prefix;
242: PetscFunctionBegin;
243: PetscCall(MatShellGetContext(J, &ctx));
244: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
245: if (isascii) {
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 PetscDefined(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_SINGLE_LIBRARY_VISIBILITY_INTERNAL 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 (PetscDefined(USE_COMPLEX)) PetscCall(PetscInfo(mat, "Current differencing parameter: %g + %g i\n", (double)PetscRealPart(h), (double)PetscImaginaryPart(h)));
342: else PetscCall(PetscInfo(mat, "Current differencing parameter: %15.12e\n", (double)PetscRealPart(h)));
343: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) ctx->historyh[ctx->ncurrenth] = h;
344: ctx->ncurrenth++;
346: #if PetscDefined(USE_COMPLEX)
347: if (ctx->usecomplex) h = PETSC_i * h;
348: #endif
350: /* w = u + ha */
351: PetscCall(VecWAXPY(w, h, a, U));
353: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
354: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) PetscCall((*ctx->func)(ctx->funcctx, U, F));
355: PetscCall((*ctx->func)(ctx->funcctx, w, y));
357: #if PetscDefined(USE_COMPLEX)
358: if (ctx->usecomplex) {
359: PetscCall(VecImaginaryPart(y));
360: h = PetscImaginaryPart(h);
361: } else {
362: PetscCall(VecAXPY(y, -1.0, F));
363: }
364: #else
365: PetscCall(VecAXPY(y, -1.0, F));
366: #endif
367: PetscCall(VecScale(y, 1.0 / h));
368: if (mat->nullsp) PetscCall(MatNullSpaceRemove(mat->nullsp, y));
370: PetscCall(PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: /*
375: MatGetDiagonal_MFFD - Gets the diagonal for a matrix-free matrix
377: y ~= (F(u + ha) - F(u))/h,
378: where F = nonlinear function, as set by SNESSetFunction()
379: u = current iterate
380: h = difference interval
381: */
382: static PetscErrorCode MatGetDiagonal_MFFD(Mat mat, Vec a)
383: {
384: MatMFFD ctx;
385: PetscScalar h, *aa, *ww, v;
386: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
387: Vec w, U;
388: PetscInt i, rstart, rend;
390: PetscFunctionBegin;
391: PetscCall(MatShellGetContext(mat, &ctx));
392: PetscCheck(ctx->func, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Requires calling MatMFFDSetFunction() first");
393: PetscCheck(ctx->funci, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Requires calling MatMFFDSetFunctioni() first");
394: w = ctx->w;
395: U = ctx->current_u;
396: PetscCall((*ctx->func)(ctx->funcctx, U, a));
397: if (ctx->funcisetbase) PetscCall((*ctx->funcisetbase)(ctx->funcctx, U));
398: PetscCall(VecCopy(U, w));
400: PetscCall(VecGetOwnershipRange(a, &rstart, &rend));
401: PetscCall(VecGetArray(a, &aa));
402: for (i = rstart; i < rend; i++) {
403: PetscCall(VecGetArray(w, &ww));
404: h = ww[i - rstart];
405: if (h == 0.0) h = 1.0;
406: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
407: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
408: h *= epsilon;
410: ww[i - rstart] += h;
411: PetscCall(VecRestoreArray(w, &ww));
412: PetscCall((*ctx->funci)(ctx->funcctx, i, w, &v));
413: aa[i - rstart] = (v - aa[i - rstart]) / h;
415: PetscCall(VecGetArray(w, &ww));
416: ww[i - rstart] -= h;
417: PetscCall(VecRestoreArray(w, &ww));
418: }
419: PetscCall(VecRestoreArray(a, &aa));
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode MatMFFDSetBase_MFFD(Mat J, Vec U, Vec F)
424: {
425: MatMFFD ctx;
427: PetscFunctionBegin;
428: PetscCall(MatShellGetContext(J, &ctx));
429: PetscCall(MatMFFDResetHHistory(J));
430: if (!ctx->current_u) {
431: PetscCall(VecDuplicate(U, &ctx->current_u));
432: PetscCall(VecLockReadPush(ctx->current_u));
433: }
434: PetscCall(VecLockReadPop(ctx->current_u));
435: PetscCall(VecCopy(U, ctx->current_u));
436: PetscCall(VecLockReadPush(ctx->current_u));
437: if (F) {
438: if (ctx->current_f_allocated) PetscCall(VecDestroy(&ctx->current_f));
439: ctx->current_f = F;
440: ctx->current_f_allocated = PETSC_FALSE;
441: } else if (!ctx->current_f_allocated) {
442: PetscCall(MatCreateVecs(J, NULL, &ctx->current_f));
443: ctx->current_f_allocated = PETSC_TRUE;
444: }
445: if (!ctx->w) PetscCall(VecDuplicate(ctx->current_u, &ctx->w));
446: J->assembled = PETSC_TRUE;
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: typedef PetscErrorCode (*FCN3)(void *, Vec, Vec, PetscScalar *); /* force argument to next function to not be extern C*/
451: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J, FCN3 fun, void *ectx)
452: {
453: MatMFFD ctx;
455: PetscFunctionBegin;
456: PetscCall(MatShellGetContext(J, &ctx));
457: ctx->checkh = fun;
458: ctx->checkhctx = ectx;
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: /*@
463: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
464: MATMFFD` options in the database.
466: Collective
468: Input Parameters:
469: + mat - the `MATMFFD` context
470: - prefix - the prefix to prepend to all option names
472: Note:
473: A hyphen (-) must NOT be given at the beginning of the prefix name.
474: The first character of all runtime options is AUTOMATICALLY the hyphen.
476: Level: advanced
478: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatSetFromOptions()`, `MatCreateSNESMF()`, `MatCreateMFFD()`
479: @*/
480: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat, const char prefix[])
481: {
482: MatMFFD mfctx;
484: PetscFunctionBegin;
486: PetscCall(MatShellGetContext(mat, &mfctx));
488: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mfctx, prefix));
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }
492: static PetscErrorCode MatSetFromOptions_MFFD(Mat mat, PetscOptionItems PetscOptionsObject)
493: {
494: MatMFFD mfctx;
495: PetscBool flg;
496: char ftype[256];
498: PetscFunctionBegin;
499: PetscCall(MatShellGetContext(mat, &mfctx));
501: PetscObjectOptionsBegin((PetscObject)mfctx);
502: PetscCall(PetscOptionsFList("-mat_mffd_type", "Matrix free type", "MatMFFDSetType", MatMFFDList, ((PetscObject)mfctx)->type_name, ftype, 256, &flg));
503: if (flg) PetscCall(MatMFFDSetType(mat, ftype));
505: PetscCall(PetscOptionsReal("-mat_mffd_err", "set sqrt relative error in function", "MatMFFDSetFunctionError", mfctx->error_rel, &mfctx->error_rel, NULL));
506: PetscCall(PetscOptionsInt("-mat_mffd_period", "how often h is recomputed", "MatMFFDSetPeriod", mfctx->recomputeperiod, &mfctx->recomputeperiod, NULL));
508: flg = PETSC_FALSE;
509: PetscCall(PetscOptionsBool("-mat_mffd_check_positivity", "Insure that U + h*a is nonnegative", "MatMFFDSetCheckh", flg, &flg, NULL));
510: if (flg) PetscCall(MatMFFDSetCheckh(mat, MatMFFDCheckPositivity, NULL));
511: #if PetscDefined(USE_COMPLEX)
512: PetscCall(PetscOptionsBool("-mat_mffd_complex", "Use Lyness complex number trick to compute the matrix-vector product", "None", mfctx->usecomplex, &mfctx->usecomplex, NULL));
513: #endif
514: PetscTryTypeMethod(mfctx, setfromoptions, PetscOptionsObject);
515: PetscOptionsEnd();
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat, PetscInt period)
520: {
521: MatMFFD ctx;
523: PetscFunctionBegin;
524: PetscCall(MatShellGetContext(mat, &ctx));
525: ctx->recomputeperiod = period;
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat, MatMFFDFn *func, void *funcctx)
530: {
531: MatMFFD ctx;
533: PetscFunctionBegin;
534: PetscCall(MatShellGetContext(mat, &ctx));
535: ctx->func = func;
536: ctx->funcctx = funcctx;
537: PetscFunctionReturn(PETSC_SUCCESS);
538: }
540: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat, PetscReal error)
541: {
542: PetscFunctionBegin;
543: if (error != (PetscReal)PETSC_DEFAULT) {
544: MatMFFD ctx;
546: PetscCall(MatShellGetContext(mat, &ctx));
547: ctx->error_rel = error;
548: }
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode MatMFFDSetHHistory_MFFD(Mat J, PetscScalar history[], PetscInt nhistory)
553: {
554: MatMFFD ctx;
556: PetscFunctionBegin;
557: PetscCall(MatShellGetContext(J, &ctx));
558: ctx->historyh = history;
559: ctx->maxcurrenth = nhistory;
560: ctx->currenth = 0.;
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: /*MC
565: MATMFFD - "mffd" - A matrix-free matrix type.
567: Level: advanced
569: Developer Notes:
570: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
572: Users should not MatShell... operations on this class, there is some error checking for that incorrect usage
574: .seealso: [](ch_matrices), `Mat`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MatMFFDSetFunction()`, `MatMFFDSetType()`,
575: `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`,
576: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`,
577: `MatMFFDGetH()`
578: M*/
579: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
580: {
581: MatMFFD mfctx;
583: PetscFunctionBegin;
584: PetscCall(MatMFFDInitializePackage());
586: PetscCall(PetscHeaderCreate(mfctx, MATMFFD_CLASSID, "MatMFFD", "Matrix-free Finite Differencing", "Mat", PetscObjectComm((PetscObject)A), NULL, NULL));
588: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
589: mfctx->recomputeperiod = 1;
590: mfctx->count = 0;
591: mfctx->currenth = 0.0;
592: mfctx->historyh = NULL;
593: mfctx->ncurrenth = 0;
594: mfctx->maxcurrenth = 0;
595: ((PetscObject)mfctx)->type_name = NULL;
597: /*
598: Create the empty data structure to contain compute-h routines.
599: These will be filled in below from the command line options or
600: a later call with MatMFFDSetType() or if that is not called
601: then it will default in the first use of MatMult_MFFD()
602: */
603: mfctx->ops->compute = NULL;
604: mfctx->ops->destroy = NULL;
605: mfctx->ops->view = NULL;
606: mfctx->ops->setfromoptions = NULL;
607: mfctx->hctx = NULL;
609: mfctx->func = NULL;
610: mfctx->funcctx = NULL;
611: mfctx->w = NULL;
612: mfctx->mat = A;
614: PetscCall(MatSetType(A, MATSHELL));
615: PetscCall(MatShellSetContext(A, mfctx));
616: PetscCall(MatShellSetOperation(A, MATOP_MULT, (PetscErrorCodeFn *)MatMult_MFFD));
617: PetscCall(MatShellSetOperation(A, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_MFFD));
618: PetscCall(MatShellSetOperation(A, MATOP_VIEW, (PetscErrorCodeFn *)MatView_MFFD));
619: PetscCall(MatShellSetOperation(A, MATOP_ASSEMBLY_END, (PetscErrorCodeFn *)MatAssemblyEnd_MFFD));
620: PetscCall(MatShellSetOperation(A, MATOP_SET_FROM_OPTIONS, (PetscErrorCodeFn *)MatSetFromOptions_MFFD));
621: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Immutable));
622: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
623: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
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` that uses finite differences on a provided function to
642: approximately multiply a vector by the matrix (Jacobian) . See also `MatCreateSNESMF()`
644: Collective
646: Input Parameters:
647: + comm - MPI communicator
648: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
649: This value should be the same as the local size used in creating the
650: y vector for the matrix-vector product y = Ax.
651: . n - This value should be the same as the local size used in creating the
652: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
653: calculated if `N` is given) For square matrices `n` is almost always `m`.
654: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
655: - N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
657: Output Parameter:
658: . J - the matrix-free matrix
660: Options Database Keys:
661: + -mat_mffd_type - wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
662: . -mat_mffd_err - square root of estimated relative error in function evaluation
663: . -mat_mffd_period - how often h is recomputed, defaults to 1, every time
664: . -mat_mffd_check_positivity - possibly decrease `h` until U + h*a has only positive values
665: . -mat_mffd_umin umin - Sets umin (for default PETSc routine that computes h only)
666: . -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
667: (requires real valued functions but that PETSc be configured for complex numbers)
668: . -snes_mf - use the finite difference based matrix-free matrix with `SNESSolve()` and no preconditioner
669: - -snes_mf_operator - use the finite difference based matrix-free matrix with `SNESSolve()` but construct a preconditioner
670: using the matrix passed as `pmat` to `SNESSetJacobian()`.
672: Level: advanced
674: Notes:
675: Use `MatMFFDSetFunction()` to provide the function that will be differenced to compute the matrix-vector product.
677: The matrix-free matrix context contains the function pointers
678: and work space for performing finite difference approximations of
679: Jacobian-vector products, F'(u)*a,
681: The default code uses the following approach to compute h
683: .vb
684: F'(u)*a = [F(u+h*a) - F(u)]/h where
685: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
686: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
687: where
688: error_rel = square root of relative error in function evaluation
689: umin = minimum iterate parameter
690: .ve
692: To have `SNES` use the matrix-free finite difference matrix-vector product and not provide a separate matrix
693: from which to compute the preconditioner (the `pmat` argument `SNESSetJacobian()`), then simply call `SNESSetJacobian()`
694: with `NULL` for the matrices and `MatMFFDComputeJacobian()`. Or use the options database option `-snes_mf`
696: The user can set `error_rel` via `MatMFFDSetFunctionError()` and `umin` via `MatMFFDDSSetUmin()`.
698: Use `MATSHELL` or `MatCreateShell()` to provide your own custom matrix-vector operation.
700: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatDestroy()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`,
701: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`, `MatCreateShell()`, `MATSHELL`,
702: `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`
703: @*/
704: PetscErrorCode MatCreateMFFD(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, Mat *J)
705: {
706: PetscFunctionBegin;
707: PetscCall(MatCreate(comm, J));
708: PetscCall(MatSetSizes(*J, m, n, M, N));
709: PetscCall(MatSetType(*J, MATMFFD));
710: PetscCall(MatSetUp(*J));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: /*@
715: MatMFFDGetH - Gets the last value that was used as the differencing for a `MATMFFD` matrix
716: parameter.
718: Not Collective
720: Input Parameters:
721: . mat - the `MATMFFD` matrix
723: Output Parameter:
724: . h - the differencing step size
726: Level: advanced
728: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDSetHHistory()`, `MatCreateMFFD()`, `MatMFFDResetHHistory()`
729: @*/
730: PetscErrorCode MatMFFDGetH(Mat mat, PetscScalar *h)
731: {
732: PetscFunctionBegin;
734: PetscAssertPointer(h, 2);
735: PetscUseMethod(mat, "MatMFFDGetH_C", (Mat, PetscScalar *), (mat, h));
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: /*@C
740: MatMFFDSetFunction - Sets the function used in applying the matrix-free `MATMFFD` matrix.
742: Logically Collective
744: Input Parameters:
745: + mat - the matrix-free matrix `MATMFFD` created via `MatCreateSNESMF()` or `MatCreateMFFD()`
746: . func - the function to use
747: - funcctx - optional function context passed to 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`, `MatMFFDFn`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`,
758: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESSetFunction()`
759: @*/
760: PetscErrorCode MatMFFDSetFunction(Mat mat, MatMFFDFn *func, void *funcctx)
761: {
762: PetscFunctionBegin;
764: PetscTryMethod(mat, "MatMFFDSetFunction_C", (Mat, MatMFFDFn *, void *), (mat, func, funcctx));
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: /*@C
769: MatMFFDSetFunctioni - Sets the function for computing 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`, `MatMFFDiFn`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`,
787: `SNESSetFunction()`, `MatGetDiagonal()`
788: @*/
789: PetscErrorCode MatMFFDSetFunctioni(Mat mat, MatMFFDiFn *funci)
790: {
791: PetscFunctionBegin;
793: PetscTryMethod(mat, "MatMFFDSetFunctioni_C", (Mat, MatMFFDiFn *), (mat, funci));
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: /*@C
798: MatMFFDSetFunctioniBase - Sets the function to compute 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()`,
815: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESSetFunction()`, `MatGetDiagonal()`
816: @*/
817: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat, MatMFFDiBaseFn *func)
818: {
819: PetscFunctionBegin;
821: PetscTryMethod(mat, "MatMFFDSetFunctioniBase_C", (Mat, MatMFFDiBaseFn *), (mat, func));
822: PetscFunctionReturn(PETSC_SUCCESS);
823: }
825: /*@
826: MatMFFDSetPeriod - Sets how often the step-size `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 - 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()`,
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;
910: if (history) PetscAssertPointer(history, 2);
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()`
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`, see `MatMFFDCheckhFn`
981: - ctx - any context needed by the function
983: Options Database Keys:
984: . -mat_mffd_check_positivity (true|false) - Ensure 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`, `MatMFFDCheckhFn`, `MatMFFDCheckPositivity()`
995: @*/
996: PetscErrorCode MatMFFDSetCheckh(Mat J, MatMFFDCheckhFn *fun, PetscCtx ctx)
997: {
998: PetscFunctionBegin;
1000: PetscTryMethod(J, "MatMFFDSetCheckh_C", (Mat, MatMFFDCheckhFn *, 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: + dummy - context variable (unused)
1012: . U - base vector that is added to
1013: . a - vector that is added
1014: - h - scaling factor on `a`, may be changed on output
1016: Options Database Keys:
1017: . -mat_mffd_check_positivity (true|false) - Ensure 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 n;
1031: MPI_Comm comm;
1033: PetscFunctionBegin;
1036: PetscAssertPointer(h, 4);
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 (PetscInt 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: PetscCallMPI(MPIU_Allreduce(&minval, &val, 1, MPIU_REAL, MPIU_MIN, comm));
1051: if (val <= PetscAbsScalar(*h)) {
1052: val = 0.99 * val;
1053: PetscCall(PetscInfo(U, "Scaling back h from %g to %g\n", (double)PetscRealPart(*h), (double)val));
1054: if (PetscRealPart(*h) > 0.0) *h = val;
1055: else *h = -val;
1056: }
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }