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: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContext_C", NULL));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*
235: MatMFFDView_MFFD - Views matrix-free parameters.
237: */
238: static PetscErrorCode MatView_MFFD(Mat J, PetscViewer viewer)
239: {
240: MatMFFD ctx;
241: PetscBool iascii, viewbase, viewfunction;
242: const char *prefix;
244: PetscFunctionBegin;
245: PetscCall(MatShellGetContext(J, &ctx));
246: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
247: if (iascii) {
248: PetscCall(PetscViewerASCIIPrintf(viewer, "Matrix-free approximation:\n"));
249: PetscCall(PetscViewerASCIIPushTab(viewer));
250: PetscCall(PetscViewerASCIIPrintf(viewer, "err=%g (relative error in function evaluation)\n", (double)ctx->error_rel));
251: if (!((PetscObject)ctx)->type_name) {
252: PetscCall(PetscViewerASCIIPrintf(viewer, "The compute h routine has not yet been set\n"));
253: } else {
254: PetscCall(PetscViewerASCIIPrintf(viewer, "Using %s compute h routine\n", ((PetscObject)ctx)->type_name));
255: }
256: #if defined(PETSC_USE_COMPLEX)
257: if (ctx->usecomplex) PetscCall(PetscViewerASCIIPrintf(viewer, "Using Lyness complex number trick to compute the matrix-vector product\n"));
258: #endif
259: PetscTryTypeMethod(ctx, view, viewer);
260: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)J, &prefix));
262: PetscCall(PetscOptionsHasName(((PetscObject)J)->options, prefix, "-mat_mffd_view_base", &viewbase));
263: if (viewbase) {
264: PetscCall(PetscViewerASCIIPrintf(viewer, "Base:\n"));
265: PetscCall(VecView(ctx->current_u, viewer));
266: }
267: PetscCall(PetscOptionsHasName(((PetscObject)J)->options, prefix, "-mat_mffd_view_function", &viewfunction));
268: if (viewfunction) {
269: PetscCall(PetscViewerASCIIPrintf(viewer, "Function:\n"));
270: PetscCall(VecView(ctx->current_f, viewer));
271: }
272: PetscCall(PetscViewerASCIIPopTab(viewer));
273: }
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: /*
278: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
279: allows the user to indicate the beginning of a new linear solve by calling
280: MatAssemblyXXX() on the matrix-free matrix. This then allows the
281: MatCreateMFFD_WP() to properly compute ||U|| only the first time
282: in the linear solver rather than every time.
284: This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
285: it must be labeled as PETSC_EXTERN
286: */
287: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J, MatAssemblyType mt)
288: {
289: MatMFFD j;
291: PetscFunctionBegin;
292: PetscCall(MatShellGetContext(J, &j));
293: PetscCall(MatMFFDResetHHistory(J));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /*
298: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
300: y ~= (F(u + ha) - F(u))/h,
301: where F = nonlinear function, as set by SNESSetFunction()
302: u = current iterate
303: h = difference interval
304: */
305: static PetscErrorCode MatMult_MFFD(Mat mat, Vec a, Vec y)
306: {
307: MatMFFD ctx;
308: PetscScalar h;
309: Vec w, U, F;
310: PetscBool zeroa;
312: PetscFunctionBegin;
313: PetscCall(MatShellGetContext(mat, &ctx));
314: 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");
315: /* We log matrix-free matrix-vector products separately, so that we can
316: separate the performance monitoring from the cases that use conventional
317: storage. We may eventually modify event logging to associate events
318: with particular objects, hence alleviating the more general problem. */
319: PetscCall(PetscLogEventBegin(MATMFFD_Mult, a, y, 0, 0));
321: w = ctx->w;
322: U = ctx->current_u;
323: F = ctx->current_f;
324: /*
325: Compute differencing parameter
326: */
327: if (!((PetscObject)ctx)->type_name) {
328: PetscCall(MatMFFDSetType(mat, MATMFFD_WP));
329: PetscCall(MatSetFromOptions(mat));
330: }
331: PetscUseTypeMethod(ctx, compute, U, a, &h, &zeroa);
332: if (zeroa) {
333: PetscCall(VecSet(y, 0.0));
334: PetscCall(PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: PetscCheck(!mat->erroriffailure || !PetscIsInfOrNanScalar(h), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Computed Nan differencing parameter h");
339: if (ctx->checkh) PetscCall((*ctx->checkh)(ctx->checkhctx, U, a, &h));
341: /* keep a record of the current differencing parameter h */
342: ctx->currenth = h;
343: #if defined(PETSC_USE_COMPLEX)
344: PetscCall(PetscInfo(mat, "Current differencing parameter: %g + %g i\n", (double)PetscRealPart(h), (double)PetscImaginaryPart(h)));
345: #else
346: PetscCall(PetscInfo(mat, "Current differencing parameter: %15.12e\n", (double)PetscRealPart(h)));
347: #endif
348: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) ctx->historyh[ctx->ncurrenth] = h;
349: ctx->ncurrenth++;
351: #if defined(PETSC_USE_COMPLEX)
352: if (ctx->usecomplex) h = PETSC_i * h;
353: #endif
355: /* w = u + ha */
356: PetscCall(VecWAXPY(w, h, a, U));
358: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
359: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) PetscCall((*ctx->func)(ctx->funcctx, U, F));
360: PetscCall((*ctx->func)(ctx->funcctx, w, y));
362: #if defined(PETSC_USE_COMPLEX)
363: if (ctx->usecomplex) {
364: PetscCall(VecImaginaryPart(y));
365: h = PetscImaginaryPart(h);
366: } else {
367: PetscCall(VecAXPY(y, -1.0, F));
368: }
369: #else
370: PetscCall(VecAXPY(y, -1.0, F));
371: #endif
372: PetscCall(VecScale(y, 1.0 / h));
373: if (mat->nullsp) PetscCall(MatNullSpaceRemove(mat->nullsp, y));
375: PetscCall(PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: /*
380: MatGetDiagonal_MFFD - Gets the diagonal for a matrix-free matrix
382: y ~= (F(u + ha) - F(u))/h,
383: where F = nonlinear function, as set by SNESSetFunction()
384: u = current iterate
385: h = difference interval
386: */
387: static PetscErrorCode MatGetDiagonal_MFFD(Mat mat, Vec a)
388: {
389: MatMFFD ctx;
390: PetscScalar h, *aa, *ww, v;
391: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
392: Vec w, U;
393: PetscInt i, rstart, rend;
395: PetscFunctionBegin;
396: PetscCall(MatShellGetContext(mat, &ctx));
397: PetscCheck(ctx->func, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Requires calling MatMFFDSetFunction() first");
398: PetscCheck(ctx->funci, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Requires calling MatMFFDSetFunctioni() first");
399: w = ctx->w;
400: U = ctx->current_u;
401: PetscCall((*ctx->func)(ctx->funcctx, U, a));
402: if (ctx->funcisetbase) PetscCall((*ctx->funcisetbase)(ctx->funcctx, U));
403: PetscCall(VecCopy(U, w));
405: PetscCall(VecGetOwnershipRange(a, &rstart, &rend));
406: PetscCall(VecGetArray(a, &aa));
407: for (i = rstart; i < rend; i++) {
408: PetscCall(VecGetArray(w, &ww));
409: h = ww[i - rstart];
410: if (h == 0.0) h = 1.0;
411: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
412: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
413: h *= epsilon;
415: ww[i - rstart] += h;
416: PetscCall(VecRestoreArray(w, &ww));
417: PetscCall((*ctx->funci)(ctx->funcctx, i, w, &v));
418: aa[i - rstart] = (v - aa[i - rstart]) / h;
420: PetscCall(VecGetArray(w, &ww));
421: ww[i - rstart] -= h;
422: PetscCall(VecRestoreArray(w, &ww));
423: }
424: PetscCall(VecRestoreArray(a, &aa));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J, Vec U, Vec F)
429: {
430: MatMFFD ctx;
432: PetscFunctionBegin;
433: PetscCall(MatShellGetContext(J, &ctx));
434: PetscCall(MatMFFDResetHHistory(J));
435: if (!ctx->current_u) {
436: PetscCall(VecDuplicate(U, &ctx->current_u));
437: PetscCall(VecLockReadPush(ctx->current_u));
438: }
439: PetscCall(VecLockReadPop(ctx->current_u));
440: PetscCall(VecCopy(U, ctx->current_u));
441: PetscCall(VecLockReadPush(ctx->current_u));
442: if (F) {
443: if (ctx->current_f_allocated) PetscCall(VecDestroy(&ctx->current_f));
444: ctx->current_f = F;
445: ctx->current_f_allocated = PETSC_FALSE;
446: } else if (!ctx->current_f_allocated) {
447: PetscCall(MatCreateVecs(J, NULL, &ctx->current_f));
448: ctx->current_f_allocated = PETSC_TRUE;
449: }
450: if (!ctx->w) PetscCall(VecDuplicate(ctx->current_u, &ctx->w));
451: J->assembled = PETSC_TRUE;
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: typedef PetscErrorCode (*FCN3)(void *, Vec, Vec, PetscScalar *); /* force argument to next function to not be extern C*/
457: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J, FCN3 fun, void *ectx)
458: {
459: MatMFFD ctx;
461: PetscFunctionBegin;
462: PetscCall(MatShellGetContext(J, &ctx));
463: ctx->checkh = fun;
464: ctx->checkhctx = ectx;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@C
469: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
470: MATMFFD` options in the database.
472: Collective
474: Input Parameters:
475: + mat - the `MATMFFD` context
476: - prefix - the prefix to prepend to all option names
478: Note:
479: A hyphen (-) must NOT be given at the beginning of the prefix name.
480: The first character of all runtime options is AUTOMATICALLY the hyphen.
482: Level: advanced
484: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatSetFromOptions()`, `MatCreateSNESMF()`, `MatCreateMFFD()`
485: @*/
486: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat, const char prefix[])
487: {
488: MatMFFD mfctx;
490: PetscFunctionBegin;
492: PetscCall(MatShellGetContext(mat, &mfctx));
494: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mfctx, prefix));
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: static PetscErrorCode MatSetFromOptions_MFFD(Mat mat, PetscOptionItems *PetscOptionsObject)
499: {
500: MatMFFD mfctx;
501: PetscBool flg;
502: char ftype[256];
504: PetscFunctionBegin;
505: PetscCall(MatShellGetContext(mat, &mfctx));
507: PetscObjectOptionsBegin((PetscObject)mfctx);
508: PetscCall(PetscOptionsFList("-mat_mffd_type", "Matrix free type", "MatMFFDSetType", MatMFFDList, ((PetscObject)mfctx)->type_name, ftype, 256, &flg));
509: if (flg) PetscCall(MatMFFDSetType(mat, ftype));
511: PetscCall(PetscOptionsReal("-mat_mffd_err", "set sqrt relative error in function", "MatMFFDSetFunctionError", mfctx->error_rel, &mfctx->error_rel, NULL));
512: PetscCall(PetscOptionsInt("-mat_mffd_period", "how often h is recomputed", "MatMFFDSetPeriod", mfctx->recomputeperiod, &mfctx->recomputeperiod, NULL));
514: flg = PETSC_FALSE;
515: PetscCall(PetscOptionsBool("-mat_mffd_check_positivity", "Insure that U + h*a is nonnegative", "MatMFFDSetCheckh", flg, &flg, NULL));
516: if (flg) PetscCall(MatMFFDSetCheckh(mat, MatMFFDCheckPositivity, NULL));
517: #if defined(PETSC_USE_COMPLEX)
518: PetscCall(PetscOptionsBool("-mat_mffd_complex", "Use Lyness complex number trick to compute the matrix-vector product", "None", mfctx->usecomplex, &mfctx->usecomplex, NULL));
519: #endif
520: PetscTryTypeMethod(mfctx, setfromoptions, PetscOptionsObject);
521: PetscOptionsEnd();
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat, PetscInt period)
526: {
527: MatMFFD ctx;
529: PetscFunctionBegin;
530: PetscCall(MatShellGetContext(mat, &ctx));
531: ctx->recomputeperiod = period;
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat, PetscErrorCode (*func)(void *, Vec, Vec), void *funcctx)
536: {
537: MatMFFD ctx;
539: PetscFunctionBegin;
540: PetscCall(MatShellGetContext(mat, &ctx));
541: ctx->func = func;
542: ctx->funcctx = funcctx;
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat, PetscReal error)
547: {
548: PetscFunctionBegin;
549: if (error != (PetscReal)PETSC_DEFAULT) {
550: MatMFFD ctx;
552: PetscCall(MatShellGetContext(mat, &ctx));
553: ctx->error_rel = error;
554: }
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: static PetscErrorCode MatMFFDSetHHistory_MFFD(Mat J, PetscScalar history[], PetscInt nhistory)
559: {
560: MatMFFD ctx;
562: PetscFunctionBegin;
563: PetscCall(MatShellGetContext(J, &ctx));
564: ctx->historyh = history;
565: ctx->maxcurrenth = nhistory;
566: ctx->currenth = 0.;
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: static PetscErrorCode MatShellSetContext_MFFD(Mat mat, void *ctx)
571: {
572: PetscFunctionBegin;
573: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the MATSHELL context for a MATMFFD, it is used by the MATMFFD");
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: static PetscErrorCode MatShellSetContextDestroy_MFFD(Mat mat, void (*f)(void))
578: {
579: PetscFunctionBegin;
580: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot set the MATSHELL context destroy for a MATMFFD, it is used by the MATMFFD");
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: /*MC
585: MATMFFD - "mffd" - A matrix-free matrix type.
587: Level: advanced
589: Developers Notes:
590: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
592: Users should not MatShell... operations on this class, there is some error checking for that incorrect usage
594: .seealso: [](ch_matrices), `Mat`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MatMFFDSetFunction()`, `MatMFFDSetType()`,
595: `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`
596: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`,
597: `MatMFFDGetH()`,
598: M*/
599: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
600: {
601: MatMFFD mfctx;
603: PetscFunctionBegin;
604: PetscCall(MatMFFDInitializePackage());
606: PetscCall(PetscHeaderCreate(mfctx, MATMFFD_CLASSID, "MatMFFD", "Matrix-free Finite Differencing", "Mat", PetscObjectComm((PetscObject)A), NULL, NULL));
608: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
609: mfctx->recomputeperiod = 1;
610: mfctx->count = 0;
611: mfctx->currenth = 0.0;
612: mfctx->historyh = NULL;
613: mfctx->ncurrenth = 0;
614: mfctx->maxcurrenth = 0;
615: ((PetscObject)mfctx)->type_name = NULL;
617: /*
618: Create the empty data structure to contain compute-h routines.
619: These will be filled in below from the command line options or
620: a later call with MatMFFDSetType() or if that is not called
621: then it will default in the first use of MatMult_MFFD()
622: */
623: mfctx->ops->compute = NULL;
624: mfctx->ops->destroy = NULL;
625: mfctx->ops->view = NULL;
626: mfctx->ops->setfromoptions = NULL;
627: mfctx->hctx = NULL;
629: mfctx->func = NULL;
630: mfctx->funcctx = NULL;
631: mfctx->w = NULL;
632: mfctx->mat = A;
634: PetscCall(MatSetType(A, MATSHELL));
635: PetscCall(MatShellSetContext(A, mfctx));
636: PetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))MatMult_MFFD));
637: PetscCall(MatShellSetOperation(A, MATOP_DESTROY, (void (*)(void))MatDestroy_MFFD));
638: PetscCall(MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MFFD));
639: PetscCall(MatShellSetOperation(A, MATOP_ASSEMBLY_END, (void (*)(void))MatAssemblyEnd_MFFD));
640: PetscCall(MatShellSetOperation(A, MATOP_SET_FROM_OPTIONS, (void (*)(void))MatSetFromOptions_MFFD));
641: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_MFFD));
642: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_MFFD));
644: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetBase_C", MatMFFDSetBase_MFFD));
645: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctioniBase_C", MatMFFDSetFunctioniBase_MFFD));
646: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctioni_C", MatMFFDSetFunctioni_MFFD));
647: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunction_C", MatMFFDSetFunction_MFFD));
648: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetCheckh_C", MatMFFDSetCheckh_MFFD));
649: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetPeriod_C", MatMFFDSetPeriod_MFFD));
650: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctionError_C", MatMFFDSetFunctionError_MFFD));
651: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDResetHHistory_C", MatMFFDResetHHistory_MFFD));
652: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetHHistory_C", MatMFFDSetHHistory_MFFD));
653: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetType_C", MatMFFDSetType_MFFD));
654: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMFFDGetH_C", MatMFFDGetH_MFFD));
655: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATMFFD));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: /*@
660: MatCreateMFFD - Creates a matrix-free matrix of type `MATMFFD`. See also `MatCreateSNESMF()`
662: Collective
664: Input Parameters:
665: + comm - MPI communicator
666: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
667: This value should be the same as the local size used in creating the
668: y vector for the matrix-vector product y = Ax.
669: . n - This value should be the same as the local size used in creating the
670: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
671: calculated if `N` is given) For square matrices `n` is almost always `m`.
672: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if `m` is given)
673: - N - number of global columns (or `PETSC_DETERMINE` to have calculated if `n` is given)
675: Output Parameter:
676: . J - the matrix-free matrix
678: Options Database Keys:
679: + -mat_mffd_type - wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
680: . -mat_mffd_err - square root of estimated relative error in function evaluation
681: . -mat_mffd_period - how often h is recomputed, defaults to 1, every time
682: . -mat_mffd_check_positivity - possibly decrease `h` until U + h*a has only positive values
683: . -mat_mffd_umin <umin> - Sets umin (for default PETSc routine that computes h only)
684: - -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
685: (requires real valued functions but that PETSc be configured for complex numbers)
687: Level: advanced
689: Notes:
690: The matrix-free matrix context merely contains the function pointers
691: and work space for performing finite difference approximations of
692: Jacobian-vector products, F'(u)*a,
694: The default code uses the following approach to compute h
696: .vb
697: F'(u)*a = [F(u+h*a) - F(u)]/h where
698: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
699: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
700: where
701: error_rel = square root of relative error in function evaluation
702: umin = minimum iterate parameter
703: .ve
705: You can call `SNESSetJacobian()` with `MatMFFDComputeJacobian()` if you are using matrix and not a different
706: preconditioner matrix
708: The user can set the error_rel via `MatMFFDSetFunctionError()` and
709: umin via `MatMFFDDSSetUmin()`.
711: The user should call `MatDestroy()` when finished with the matrix-free
712: matrix context.
714: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatDestroy()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`
715: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`,
716: `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`
717: @*/
718: PetscErrorCode MatCreateMFFD(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, Mat *J)
719: {
720: PetscFunctionBegin;
721: PetscCall(MatCreate(comm, J));
722: PetscCall(MatSetSizes(*J, m, n, M, N));
723: PetscCall(MatSetType(*J, MATMFFD));
724: PetscCall(MatSetUp(*J));
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: /*@
729: MatMFFDGetH - Gets the last value that was used as the differencing for a `MATMFFD` matrix
730: parameter.
732: Not Collective
734: Input Parameters:
735: . mat - the `MATMFFD` matrix
737: Output Parameter:
738: . h - the differencing step size
740: Level: advanced
742: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDSetHHistory()`, `MatCreateMFFD()`, `MatMFFDResetHHistory()`
743: @*/
744: PetscErrorCode MatMFFDGetH(Mat mat, PetscScalar *h)
745: {
746: PetscFunctionBegin;
748: PetscAssertPointer(h, 2);
749: PetscUseMethod(mat, "MatMFFDGetH_C", (Mat, PetscScalar *), (mat, h));
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: /*@C
754: MatMFFDSetFunction - Sets the function used in applying the matrix-free `MATMFFD` matrix.
756: Logically Collective
758: Input Parameters:
759: + mat - the matrix-free matrix `MATMFFD` created via `MatCreateSNESMF()` or `MatCreateMFFD()`
760: . func - the function to use
761: - funcctx - optional function context passed to function
763: Calling sequence of `func`:
764: + funcctx - user provided context
765: . x - input vector
766: - f - computed output function
768: Level: advanced
770: Notes:
771: If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix-free
772: matrix inside your compute Jacobian routine
774: If this is not set then it will use the function set with `SNESSetFunction()` if `MatCreateSNESMF()` was used.
776: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`,
777: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESSetFunction()`
778: @*/
779: PetscErrorCode MatMFFDSetFunction(Mat mat, PetscErrorCode (*func)(void *funcctx, Vec x, Vec f), void *funcctx)
780: {
781: PetscFunctionBegin;
783: PetscTryMethod(mat, "MatMFFDSetFunction_C", (Mat, PetscErrorCode(*)(void *, Vec, Vec), void *), (mat, func, funcctx));
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
787: /*@C
788: MatMFFDSetFunctioni - Sets the function for a single component for a `MATMFFD` matrix
790: Logically Collective
792: Input Parameters:
793: + mat - the matrix-free matrix `MATMFFD`
794: - funci - the function to use
796: Level: advanced
798: Notes:
799: If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix-free
800: matrix inside your compute Jacobian routine.
802: This function is necessary to compute the diagonal of the matrix.
803: funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.
805: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`, `MatGetDiagonal()`
806: @*/
807: PetscErrorCode MatMFFDSetFunctioni(Mat mat, PetscErrorCode (*funci)(void *, PetscInt, Vec, PetscScalar *))
808: {
809: PetscFunctionBegin;
811: PetscTryMethod(mat, "MatMFFDSetFunctioni_C", (Mat, PetscErrorCode(*)(void *, PetscInt, Vec, PetscScalar *)), (mat, funci));
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: /*@C
816: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation for a `MATMFFD` matrix
818: Logically Collective
820: Input Parameters:
821: + mat - the `MATMFFD` matrix-free matrix
822: - func - the function to use
824: Level: advanced
826: Notes:
827: If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix-free
828: matrix inside your compute Jacobian routine.
830: This function is necessary to compute the diagonal of the matrix, used for example with `PCJACOBI`
832: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`,
833: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`, `MatGetDiagonal()`
834: @*/
835: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat, PetscErrorCode (*func)(void *, Vec))
836: {
837: PetscFunctionBegin;
839: PetscTryMethod(mat, "MatMFFDSetFunctioniBase_C", (Mat, PetscErrorCode(*)(void *, Vec)), (mat, func));
840: PetscFunctionReturn(PETSC_SUCCESS);
841: }
843: /*@
844: MatMFFDSetPeriod - Sets how often h is recomputed for a `MATMFFD` matrix, by default it is every time
846: Logically Collective
848: Input Parameters:
849: + mat - the `MATMFFD` matrix-free matrix
850: - period - 1 for every time, 2 for every second etc
852: Options Database Key:
853: . -mat_mffd_period <period> - Sets how often `h` is recomputed
855: Level: advanced
857: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`,
858: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
859: @*/
860: PetscErrorCode MatMFFDSetPeriod(Mat mat, PetscInt period)
861: {
862: PetscFunctionBegin;
865: PetscTryMethod(mat, "MatMFFDSetPeriod_C", (Mat, PetscInt), (mat, period));
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: /*@
870: MatMFFDSetFunctionError - Sets the error_rel for the approximation of matrix-vector products using finite differences with the `MATMFFD` matrix
872: Logically Collective
874: Input Parameters:
875: + mat - the `MATMFFD` matrix-free matrix
876: - error - relative error (should be set to the square root of the relative error in the function evaluations)
878: Options Database Key:
879: . -mat_mffd_err <error_rel> - Sets error_rel
881: Level: advanced
883: Note:
884: The default matrix-free matrix-vector product routine computes
885: .vb
886: F'(u)*a = [F(u+h*a) - F(u)]/h where
887: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
888: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
889: .ve
891: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`,
892: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
893: @*/
894: PetscErrorCode MatMFFDSetFunctionError(Mat mat, PetscReal error)
895: {
896: PetscFunctionBegin;
899: PetscTryMethod(mat, "MatMFFDSetFunctionError_C", (Mat, PetscReal), (mat, error));
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: /*@
904: MatMFFDSetHHistory - Sets an array to collect a history of the
905: differencing values (h) computed for the matrix-free product `MATMFFD` matrix
907: Logically Collective
909: Input Parameters:
910: + J - the `MATMFFD` matrix-free matrix
911: . history - space to hold the history
912: - nhistory - number of entries in history, if more entries are generated than
913: nhistory, then the later ones are discarded
915: Level: advanced
917: Note:
918: Use `MatMFFDResetHHistory()` to reset the history counter and collect
919: a new batch of differencing parameters, h.
921: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
922: `MatMFFDResetHHistory()`, `MatMFFDSetFunctionError()`
923: @*/
924: PetscErrorCode MatMFFDSetHHistory(Mat J, PetscScalar history[], PetscInt nhistory)
925: {
926: PetscFunctionBegin;
928: if (history) PetscAssertPointer(history, 2);
930: PetscUseMethod(J, "MatMFFDSetHHistory_C", (Mat, PetscScalar[], PetscInt), (J, history, nhistory));
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: /*@
935: MatMFFDResetHHistory - Resets the counter to zero to begin
936: collecting a new set of differencing histories for the `MATMFFD` matrix
938: Logically Collective
940: Input Parameter:
941: . J - the matrix-free matrix context
943: Level: advanced
945: Note:
946: Use `MatMFFDSetHHistory()` to create the original history counter.
948: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
949: `MatMFFDSetHHistory()`, `MatMFFDSetFunctionError()`
950: @*/
951: PetscErrorCode MatMFFDResetHHistory(Mat J)
952: {
953: PetscFunctionBegin;
955: PetscTryMethod(J, "MatMFFDResetHHistory_C", (Mat), (J));
956: PetscFunctionReturn(PETSC_SUCCESS);
957: }
959: /*@
960: MatMFFDSetBase - Sets the vector `U` at which matrix vector products of the
961: Jacobian are computed for the `MATMFFD` matrix
963: Logically Collective
965: Input Parameters:
966: + J - the `MATMFFD` matrix
967: . U - the vector
968: - F - (optional) vector that contains F(u) if it has been already computed
970: Level: advanced
972: Notes:
973: This is rarely used directly
975: If `F` is provided then it is not recomputed. Otherwise the function is evaluated at the base
976: point during the first `MatMult()` after each call to `MatMFFDSetBase()`.
978: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMult()`
979: @*/
980: PetscErrorCode MatMFFDSetBase(Mat J, Vec U, Vec F)
981: {
982: PetscFunctionBegin;
986: PetscTryMethod(J, "MatMFFDSetBase_C", (Mat, Vec, Vec), (J, U, F));
987: PetscFunctionReturn(PETSC_SUCCESS);
988: }
990: /*@C
991: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
992: it to satisfy some criteria for the `MATMFFD` matrix
994: Logically Collective
996: Input Parameters:
997: + J - the `MATMFFD` matrix
998: . fun - the function that checks `h`
999: - ctx - any context needed by the function
1001: Options Database Keys:
1002: . -mat_mffd_check_positivity <bool> - Insure that U + h*a is non-negative
1004: Level: advanced
1006: Notes:
1007: For example, `MatMFFDCheckPositivity()` insures that all entries of U + h*a are non-negative
1009: The function you provide is called after the default `h` has been computed and allows you to
1010: modify it.
1012: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDCheckPositivity()`
1013: @*/
1014: PetscErrorCode MatMFFDSetCheckh(Mat J, PetscErrorCode (*fun)(void *, Vec, Vec, PetscScalar *), void *ctx)
1015: {
1016: PetscFunctionBegin;
1018: PetscTryMethod(J, "MatMFFDSetCheckh_C", (Mat, PetscErrorCode(*)(void *, Vec, Vec, PetscScalar *), void *), (J, fun, ctx));
1019: PetscFunctionReturn(PETSC_SUCCESS);
1020: }
1022: /*@
1023: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1024: zero, decreases h until this is satisfied for a `MATMFFD` matrix
1026: Logically Collective
1028: Input Parameters:
1029: + U - base vector that is added to
1030: . a - vector that is added
1031: . h - scaling factor on a
1032: - dummy - context variable (unused)
1034: Options Database Keys:
1035: . -mat_mffd_check_positivity <bool> - Insure that U + h*a is nonnegative
1037: Level: advanced
1039: Note:
1040: This is rarely used directly, rather it is passed as an argument to `MatMFFDSetCheckh()`
1042: .seealso: [](ch_matrices), `Mat`, `MATMFFD`, `MatMFFDSetCheckh()`
1043: @*/
1044: PetscErrorCode MatMFFDCheckPositivity(void *dummy, Vec U, Vec a, PetscScalar *h)
1045: {
1046: PetscReal val, minval;
1047: PetscScalar *u_vec, *a_vec;
1048: PetscInt i, n;
1049: MPI_Comm comm;
1051: PetscFunctionBegin;
1054: PetscAssertPointer(h, 4);
1055: PetscCall(PetscObjectGetComm((PetscObject)U, &comm));
1056: PetscCall(VecGetArray(U, &u_vec));
1057: PetscCall(VecGetArray(a, &a_vec));
1058: PetscCall(VecGetLocalSize(U, &n));
1059: minval = PetscAbsScalar(*h) * PetscRealConstant(1.01);
1060: for (i = 0; i < n; i++) {
1061: if (PetscRealPart(u_vec[i] + *h * a_vec[i]) <= 0.0) {
1062: val = PetscAbsScalar(u_vec[i] / a_vec[i]);
1063: if (val < minval) minval = val;
1064: }
1065: }
1066: PetscCall(VecRestoreArray(U, &u_vec));
1067: PetscCall(VecRestoreArray(a, &a_vec));
1068: PetscCall(MPIU_Allreduce(&minval, &val, 1, MPIU_REAL, MPIU_MIN, comm));
1069: if (val <= PetscAbsScalar(*h)) {
1070: PetscCall(PetscInfo(U, "Scaling back h from %g to %g\n", (double)PetscRealPart(*h), (double)(.99 * val)));
1071: if (PetscRealPart(*h) > 0.0) *h = 0.99 * val;
1072: else *h = -0.99 * val;
1073: }
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }