Actual source code: mffd.c
2: #include <petsc/private/matimpl.h>
3: #include <../src/mat/impls/mffd/mffdimpl.h>
5: PetscFunctionList MatMFFDList = NULL;
6: PetscBool MatMFFDRegisterAllCalled = PETSC_FALSE;
8: PetscClassId MATMFFD_CLASSID;
9: PetscLogEvent MATMFFD_Mult;
11: static PetscBool MatMFFDPackageInitialized = PETSC_FALSE;
12: /*@C
13: MatMFFDFinalizePackage - This function destroys everything in the MATMFFD` package. It is
14: called from `PetscFinalize()`.
16: Level: developer
18: .seealso: `MATMFFD`, `PetscFinalize()`, `MatCreateMFFD()`, `MatCreateSNESMF()`
19: @*/
20: PetscErrorCode MatMFFDFinalizePackage(void)
21: {
22: PetscFunctionListDestroy(&MatMFFDList);
23: MatMFFDPackageInitialized = PETSC_FALSE;
24: MatMFFDRegisterAllCalled = PETSC_FALSE;
25: return 0;
26: }
28: /*@C
29: MatMFFDInitializePackage - This function initializes everything in the MATMFFD` package. It is called
30: from `MatInitializePackage()`.
32: Level: developer
34: .seealso: `MATMFFD`, `PetscInitialize()`
35: @*/
36: PetscErrorCode MatMFFDInitializePackage(void)
37: {
38: char logList[256];
39: PetscBool opt, pkg;
41: if (MatMFFDPackageInitialized) return 0;
42: MatMFFDPackageInitialized = PETSC_TRUE;
43: /* Register Classes */
44: PetscClassIdRegister("MatMFFD", &MATMFFD_CLASSID);
45: /* Register Constructors */
46: MatMFFDRegisterAll();
47: /* Register Events */
48: PetscLogEventRegister("MatMult MF", MATMFFD_CLASSID, &MATMFFD_Mult);
49: /* Process Info */
50: {
51: PetscClassId classids[1];
53: classids[0] = MATMFFD_CLASSID;
54: PetscInfoProcessClass("matmffd", 1, classids);
55: }
56: /* Process summary exclusions */
57: PetscOptionsGetString(NULL, NULL, "-log_exclude", logList, sizeof(logList), &opt);
58: if (opt) {
59: PetscStrInList("matmffd", logList, ',', &pkg);
60: if (pkg) PetscLogEventExcludeClass(MATMFFD_CLASSID);
61: }
62: /* Register package finalizer */
63: PetscRegisterFinalize(MatMFFDFinalizePackage);
64: return 0;
65: }
67: static PetscErrorCode MatMFFDSetType_MFFD(Mat mat, MatMFFDType ftype)
68: {
69: MatMFFD ctx;
70: PetscBool match;
71: PetscErrorCode (*r)(MatMFFD);
75: MatShellGetContext(mat, &ctx);
77: /* already set, so just return */
78: PetscObjectTypeCompare((PetscObject)ctx, ftype, &match);
79: if (match) return 0;
81: /* destroy the old one if it exists */
82: PetscTryTypeMethod(ctx, destroy);
84: PetscFunctionListFind(MatMFFDList, ftype, &r);
86: (*r)(ctx);
87: PetscObjectChangeTypeName((PetscObject)ctx, ftype);
88: return 0;
89: }
91: /*@C
92: MatMFFDSetType - Sets the method that is used to compute the
93: differencing parameter for finite difference matrix-free formulations.
95: Input Parameters:
96: + mat - the "matrix-free" matrix created via `MatCreateSNESMF()`, or `MatCreateMFFD()`
97: or `MatSetType`(mat,`MATMFFD`);
98: - ftype - the type requested, either `MATMFFD_WP` or `MATMFFD_DS`
100: Level: advanced
102: Note:
103: For example, such routines can compute h for use in
104: Jacobian-vector products of the form
106: F(x+ha) - F(x)
107: F'(u)a ~= ----------------
108: h
110: .seealso: `MATMFFD`, `MATMFFD_WP`, `MATMFFD_DS`, `MatCreateSNESMF()`, `MatMFFDRegister()`, `MatMFFDSetFunction()`, `MatCreateMFFD()`
111: @*/
112: PetscErrorCode MatMFFDSetType(Mat mat, MatMFFDType ftype)
113: {
116: PetscTryMethod(mat, "MatMFFDSetType_C", (Mat, MatMFFDType), (mat, ftype));
117: return 0;
118: }
120: static PetscErrorCode MatGetDiagonal_MFFD(Mat, Vec);
122: typedef PetscErrorCode (*FCN1)(void *, Vec); /* force argument to next function to not be extern C*/
123: static PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat, FCN1 func)
124: {
125: MatMFFD ctx;
127: MatShellGetContext(mat, &ctx);
128: ctx->funcisetbase = func;
129: return 0;
130: }
132: typedef PetscErrorCode (*FCN2)(void *, PetscInt, Vec, PetscScalar *); /* force argument to next function to not be extern C*/
133: static PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat, FCN2 funci)
134: {
135: MatMFFD ctx;
137: MatShellGetContext(mat, &ctx);
138: ctx->funci = funci;
139: MatShellSetOperation(mat, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_MFFD);
140: return 0;
141: }
143: static PetscErrorCode MatMFFDGetH_MFFD(Mat mat, PetscScalar *h)
144: {
145: MatMFFD ctx;
147: MatShellGetContext(mat, &ctx);
148: *h = ctx->currenth;
149: return 0;
150: }
152: static PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
153: {
154: MatMFFD ctx;
156: MatShellGetContext(J, &ctx);
157: ctx->ncurrenth = 0;
158: return 0;
159: }
161: /*@C
162: MatMFFDRegister - Adds a method to the MATMFFD` registry.
164: Not Collective
166: Input Parameters:
167: + name_solver - name of a new user-defined compute-h module
168: - routine_create - routine to create method context
170: Level: developer
172: Note:
173: `MatMFFDRegister()` may be called multiple times to add several user-defined solvers.
175: Sample usage:
176: .vb
177: MatMFFDRegister("my_h",MyHCreate);
178: .ve
180: Then, your solver can be chosen with the procedural interface via
181: $ `MatMFFDSetType`(mfctx,"my_h")
182: or at runtime via the option
183: $ -mat_mffd_type my_h
185: .seealso: `MATMFFD`, `MatMFFDRegisterAll()`, `MatMFFDRegisterDestroy()`
186: @*/
187: PetscErrorCode MatMFFDRegister(const char sname[], PetscErrorCode (*function)(MatMFFD))
188: {
189: MatInitializePackage();
190: PetscFunctionListAdd(&MatMFFDList, sname, function);
191: return 0;
192: }
194: /* ----------------------------------------------------------------------------------------*/
195: static PetscErrorCode MatDestroy_MFFD(Mat mat)
196: {
197: MatMFFD ctx;
199: MatShellGetContext(mat, &ctx);
200: VecDestroy(&ctx->w);
201: VecDestroy(&ctx->current_u);
202: if (ctx->current_f_allocated) VecDestroy(&ctx->current_f);
203: PetscTryTypeMethod(ctx, destroy);
204: PetscHeaderDestroy(&ctx);
206: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetBase_C", NULL);
207: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunctioniBase_C", NULL);
208: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunctioni_C", NULL);
209: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunction_C", NULL);
210: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetFunctionError_C", NULL);
211: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetCheckh_C", NULL);
212: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetPeriod_C", NULL);
213: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDResetHHistory_C", NULL);
214: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetHHistory_C", NULL);
215: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDSetType_C", NULL);
216: PetscObjectComposeFunction((PetscObject)mat, "MatMFFDGetH_C", NULL);
217: PetscObjectComposeFunction((PetscObject)mat, "MatSNESMFSetReuseBase_C", NULL);
218: PetscObjectComposeFunction((PetscObject)mat, "MatSNESMFGetReuseBase_C", NULL);
219: return 0;
220: }
222: /*
223: MatMFFDView_MFFD - Views matrix-free parameters.
225: */
226: static PetscErrorCode MatView_MFFD(Mat J, PetscViewer viewer)
227: {
228: MatMFFD ctx;
229: PetscBool iascii, viewbase, viewfunction;
230: const char *prefix;
232: MatShellGetContext(J, &ctx);
233: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
234: if (iascii) {
235: PetscViewerASCIIPrintf(viewer, "Matrix-free approximation:\n");
236: PetscViewerASCIIPushTab(viewer);
237: PetscViewerASCIIPrintf(viewer, "err=%g (relative error in function evaluation)\n", (double)ctx->error_rel);
238: if (!((PetscObject)ctx)->type_name) {
239: PetscViewerASCIIPrintf(viewer, "The compute h routine has not yet been set\n");
240: } else {
241: PetscViewerASCIIPrintf(viewer, "Using %s compute h routine\n", ((PetscObject)ctx)->type_name);
242: }
243: #if defined(PETSC_USE_COMPLEX)
244: if (ctx->usecomplex) PetscViewerASCIIPrintf(viewer, "Using Lyness complex number trick to compute the matrix-vector product\n");
245: #endif
246: PetscTryTypeMethod(ctx, view, viewer);
247: PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);
249: PetscOptionsHasName(((PetscObject)J)->options, prefix, "-mat_mffd_view_base", &viewbase);
250: if (viewbase) {
251: PetscViewerASCIIPrintf(viewer, "Base:\n");
252: VecView(ctx->current_u, viewer);
253: }
254: PetscOptionsHasName(((PetscObject)J)->options, prefix, "-mat_mffd_view_function", &viewfunction);
255: if (viewfunction) {
256: PetscViewerASCIIPrintf(viewer, "Function:\n");
257: VecView(ctx->current_f, viewer);
258: }
259: PetscViewerASCIIPopTab(viewer);
260: }
261: return 0;
262: }
264: /*
265: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
266: allows the user to indicate the beginning of a new linear solve by calling
267: MatAssemblyXXX() on the matrix free matrix. This then allows the
268: MatCreateMFFD_WP() to properly compute ||U|| only the first time
269: in the linear solver rather than every time.
271: This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
272: it must be labeled as PETSC_EXTERN
273: */
274: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J, MatAssemblyType mt)
275: {
276: MatMFFD j;
278: MatShellGetContext(J, &j);
279: MatMFFDResetHHistory(J);
280: return 0;
281: }
283: /*
284: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
286: y ~= (F(u + ha) - F(u))/h,
287: where F = nonlinear function, as set by SNESSetFunction()
288: u = current iterate
289: h = difference interval
290: */
291: static PetscErrorCode MatMult_MFFD(Mat mat, Vec a, Vec y)
292: {
293: MatMFFD ctx;
294: PetscScalar h;
295: Vec w, U, F;
296: PetscBool zeroa;
298: MatShellGetContext(mat, &ctx);
300: /* We log matrix-free matrix-vector products separately, so that we can
301: separate the performance monitoring from the cases that use conventional
302: storage. We may eventually modify event logging to associate events
303: with particular objects, hence alleviating the more general problem. */
304: PetscLogEventBegin(MATMFFD_Mult, a, y, 0, 0);
306: w = ctx->w;
307: U = ctx->current_u;
308: F = ctx->current_f;
309: /*
310: Compute differencing parameter
311: */
312: if (!((PetscObject)ctx)->type_name) {
313: MatMFFDSetType(mat, MATMFFD_WP);
314: MatSetFromOptions(mat);
315: }
316: PetscUseTypeMethod(ctx, compute, U, a, &h, &zeroa);
317: if (zeroa) {
318: VecSet(y, 0.0);
319: PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0);
320: return 0;
321: }
324: if (ctx->checkh) (*ctx->checkh)(ctx->checkhctx, U, a, &h);
326: /* keep a record of the current differencing parameter h */
327: ctx->currenth = h;
328: #if defined(PETSC_USE_COMPLEX)
329: PetscInfo(mat, "Current differencing parameter: %g + %g i\n", (double)PetscRealPart(h), (double)PetscImaginaryPart(h));
330: #else
331: PetscInfo(mat, "Current differencing parameter: %15.12e\n", (double)PetscRealPart(h));
332: #endif
333: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) ctx->historyh[ctx->ncurrenth] = h;
334: ctx->ncurrenth++;
336: #if defined(PETSC_USE_COMPLEX)
337: if (ctx->usecomplex) h = PETSC_i * h;
338: #endif
340: /* w = u + ha */
341: VecWAXPY(w, h, a, U);
343: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
344: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) (*ctx->func)(ctx->funcctx, U, F);
345: (*ctx->func)(ctx->funcctx, w, y);
347: #if defined(PETSC_USE_COMPLEX)
348: if (ctx->usecomplex) {
349: VecImaginaryPart(y);
350: h = PetscImaginaryPart(h);
351: } else {
352: VecAXPY(y, -1.0, F);
353: }
354: #else
355: VecAXPY(y, -1.0, F);
356: #endif
357: VecScale(y, 1.0 / h);
358: if (mat->nullsp) MatNullSpaceRemove(mat->nullsp, y);
360: PetscLogEventEnd(MATMFFD_Mult, a, y, 0, 0);
361: return 0;
362: }
364: /*
365: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
367: y ~= (F(u + ha) - F(u))/h,
368: where F = nonlinear function, as set by SNESSetFunction()
369: u = current iterate
370: h = difference interval
371: */
372: PetscErrorCode MatGetDiagonal_MFFD(Mat mat, Vec a)
373: {
374: MatMFFD ctx;
375: PetscScalar h, *aa, *ww, v;
376: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON;
377: Vec w, U;
378: PetscInt i, rstart, rend;
380: MatShellGetContext(mat, &ctx);
383: w = ctx->w;
384: U = ctx->current_u;
385: (*ctx->func)(ctx->funcctx, U, a);
386: if (ctx->funcisetbase) (*ctx->funcisetbase)(ctx->funcctx, U);
387: VecCopy(U, w);
389: VecGetOwnershipRange(a, &rstart, &rend);
390: VecGetArray(a, &aa);
391: for (i = rstart; i < rend; i++) {
392: VecGetArray(w, &ww);
393: h = ww[i - rstart];
394: if (h == 0.0) h = 1.0;
395: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
396: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
397: h *= epsilon;
399: ww[i - rstart] += h;
400: VecRestoreArray(w, &ww);
401: (*ctx->funci)(ctx->funcctx, i, w, &v);
402: aa[i - rstart] = (v - aa[i - rstart]) / h;
404: VecGetArray(w, &ww);
405: ww[i - rstart] -= h;
406: VecRestoreArray(w, &ww);
407: }
408: VecRestoreArray(a, &aa);
409: return 0;
410: }
412: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J, Vec U, Vec F)
413: {
414: MatMFFD ctx;
416: MatShellGetContext(J, &ctx);
417: MatMFFDResetHHistory(J);
418: if (!ctx->current_u) {
419: VecDuplicate(U, &ctx->current_u);
420: VecLockReadPush(ctx->current_u);
421: }
422: VecLockReadPop(ctx->current_u);
423: VecCopy(U, ctx->current_u);
424: VecLockReadPush(ctx->current_u);
425: if (F) {
426: if (ctx->current_f_allocated) VecDestroy(&ctx->current_f);
427: ctx->current_f = F;
428: ctx->current_f_allocated = PETSC_FALSE;
429: } else if (!ctx->current_f_allocated) {
430: MatCreateVecs(J, NULL, &ctx->current_f);
431: ctx->current_f_allocated = PETSC_TRUE;
432: }
433: if (!ctx->w) VecDuplicate(ctx->current_u, &ctx->w);
434: J->assembled = PETSC_TRUE;
435: return 0;
436: }
438: typedef PetscErrorCode (*FCN3)(void *, Vec, Vec, PetscScalar *); /* force argument to next function to not be extern C*/
440: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J, FCN3 fun, void *ectx)
441: {
442: MatMFFD ctx;
444: MatShellGetContext(J, &ctx);
445: ctx->checkh = fun;
446: ctx->checkhctx = ectx;
447: return 0;
448: }
450: /*@C
451: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
452: MATMFFD` options in the database.
454: Collective
456: Input Parameters:
457: + A - the `MATMFFD` context
458: - prefix - the prefix to prepend to all option names
460: Note:
461: A hyphen (-) must NOT be given at the beginning of the prefix name.
462: The first character of all runtime options is AUTOMATICALLY the hyphen.
464: Level: advanced
466: .seealso: `MATMFFD`, `MatSetFromOptions()`, `MatCreateSNESMF()`, `MatCreateMFFD()`
467: @*/
468: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat, const char prefix[])
469: {
470: MatMFFD mfctx;
473: MatShellGetContext(mat, &mfctx);
475: PetscObjectSetOptionsPrefix((PetscObject)mfctx, prefix);
476: return 0;
477: }
479: static PetscErrorCode MatSetFromOptions_MFFD(Mat mat, PetscOptionItems *PetscOptionsObject)
480: {
481: MatMFFD mfctx;
482: PetscBool flg;
483: char ftype[256];
485: MatShellGetContext(mat, &mfctx);
487: PetscObjectOptionsBegin((PetscObject)mfctx);
488: PetscOptionsFList("-mat_mffd_type", "Matrix free type", "MatMFFDSetType", MatMFFDList, ((PetscObject)mfctx)->type_name, ftype, 256, &flg);
489: if (flg) MatMFFDSetType(mat, ftype);
491: PetscOptionsReal("-mat_mffd_err", "set sqrt relative error in function", "MatMFFDSetFunctionError", mfctx->error_rel, &mfctx->error_rel, NULL);
492: PetscOptionsInt("-mat_mffd_period", "how often h is recomputed", "MatMFFDSetPeriod", mfctx->recomputeperiod, &mfctx->recomputeperiod, NULL);
494: flg = PETSC_FALSE;
495: PetscOptionsBool("-mat_mffd_check_positivity", "Insure that U + h*a is nonnegative", "MatMFFDSetCheckh", flg, &flg, NULL);
496: if (flg) MatMFFDSetCheckh(mat, MatMFFDCheckPositivity, NULL);
497: #if defined(PETSC_USE_COMPLEX)
498: PetscOptionsBool("-mat_mffd_complex", "Use Lyness complex number trick to compute the matrix-vector product", "None", mfctx->usecomplex, &mfctx->usecomplex, NULL);
499: #endif
500: PetscTryTypeMethod(mfctx, setfromoptions, PetscOptionsObject);
501: PetscOptionsEnd();
502: return 0;
503: }
505: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat, PetscInt period)
506: {
507: MatMFFD ctx;
509: MatShellGetContext(mat, &ctx);
510: ctx->recomputeperiod = period;
511: return 0;
512: }
514: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat, PetscErrorCode (*func)(void *, Vec, Vec), void *funcctx)
515: {
516: MatMFFD ctx;
518: MatShellGetContext(mat, &ctx);
519: ctx->func = func;
520: ctx->funcctx = funcctx;
521: return 0;
522: }
524: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat, PetscReal error)
525: {
526: MatMFFD ctx;
528: MatShellGetContext(mat, &ctx);
529: if (error != PETSC_DEFAULT) ctx->error_rel = error;
530: return 0;
531: }
533: PetscErrorCode MatMFFDSetHHistory_MFFD(Mat J, PetscScalar history[], PetscInt nhistory)
534: {
535: MatMFFD ctx;
537: MatShellGetContext(J, &ctx);
538: ctx->historyh = history;
539: ctx->maxcurrenth = nhistory;
540: ctx->currenth = 0.;
541: return 0;
542: }
544: /*MC
545: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
547: Level: advanced
549: Developers Note: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
551: .seealso: `MatCreateMFFD()`, `MatCreateSNESMF()`, `MatMFFDSetFunction()`, `MatMFFDSetType()`,
552: `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`
553: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`,
554: `MatMFFDGetH()`,
555: M*/
556: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
557: {
558: MatMFFD mfctx;
560: MatMFFDInitializePackage();
562: PetscHeaderCreate(mfctx, MATMFFD_CLASSID, "MatMFFD", "Matrix-free Finite Differencing", "Mat", PetscObjectComm((PetscObject)A), NULL, NULL);
564: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
565: mfctx->recomputeperiod = 1;
566: mfctx->count = 0;
567: mfctx->currenth = 0.0;
568: mfctx->historyh = NULL;
569: mfctx->ncurrenth = 0;
570: mfctx->maxcurrenth = 0;
571: ((PetscObject)mfctx)->type_name = NULL;
573: /*
574: Create the empty data structure to contain compute-h routines.
575: These will be filled in below from the command line options or
576: a later call with MatMFFDSetType() or if that is not called
577: then it will default in the first use of MatMult_MFFD()
578: */
579: mfctx->ops->compute = NULL;
580: mfctx->ops->destroy = NULL;
581: mfctx->ops->view = NULL;
582: mfctx->ops->setfromoptions = NULL;
583: mfctx->hctx = NULL;
585: mfctx->func = NULL;
586: mfctx->funcctx = NULL;
587: mfctx->w = NULL;
588: mfctx->mat = A;
590: MatSetType(A, MATSHELL);
591: MatShellSetContext(A, mfctx);
592: MatShellSetOperation(A, MATOP_MULT, (void (*)(void))MatMult_MFFD);
593: MatShellSetOperation(A, MATOP_DESTROY, (void (*)(void))MatDestroy_MFFD);
594: MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MFFD);
595: MatShellSetOperation(A, MATOP_ASSEMBLY_END, (void (*)(void))MatAssemblyEnd_MFFD);
596: MatShellSetOperation(A, MATOP_SET_FROM_OPTIONS, (void (*)(void))MatSetFromOptions_MFFD);
598: PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetBase_C", MatMFFDSetBase_MFFD);
599: PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctioniBase_C", MatMFFDSetFunctioniBase_MFFD);
600: PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctioni_C", MatMFFDSetFunctioni_MFFD);
601: PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunction_C", MatMFFDSetFunction_MFFD);
602: PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetCheckh_C", MatMFFDSetCheckh_MFFD);
603: PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetPeriod_C", MatMFFDSetPeriod_MFFD);
604: PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetFunctionError_C", MatMFFDSetFunctionError_MFFD);
605: PetscObjectComposeFunction((PetscObject)A, "MatMFFDResetHHistory_C", MatMFFDResetHHistory_MFFD);
606: PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetHHistory_C", MatMFFDSetHHistory_MFFD);
607: PetscObjectComposeFunction((PetscObject)A, "MatMFFDSetType_C", MatMFFDSetType_MFFD);
608: PetscObjectComposeFunction((PetscObject)A, "MatMFFDGetH_C", MatMFFDGetH_MFFD);
609: PetscObjectChangeTypeName((PetscObject)A, MATMFFD);
610: return 0;
611: }
613: /*@
614: MatCreateMFFD - Creates a matrix-free matrix of type `MATMFFD`. See also `MatCreateSNESMF()`
616: Collective
618: Input Parameters:
619: + comm - MPI communicator
620: . m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given)
621: This value should be the same as the local size used in creating the
622: y vector for the matrix-vector product y = Ax.
623: . n - This value should be the same as the local size used in creating the
624: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
625: calculated if N is given) For square matrices n is almost always m.
626: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
627: - N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
629: Output Parameter:
630: . J - the matrix-free matrix
632: Options Database Keys:
633: + -mat_mffd_type - wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
634: . -mat_mffd_err - square root of estimated relative error in function evaluation
635: . -mat_mffd_period - how often h is recomputed, defaults to 1, every time
636: . -mat_mffd_check_positivity - possibly decrease h until U + h*a has only positive values
637: . -mat_mffd_umin <umin> - Sets umin (for default PETSc routine that computes h only)
638: - -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
639: (requires real valued functions but that PETSc be configured for complex numbers)
641: Level: advanced
643: Notes:
644: The matrix-free matrix context merely contains the function pointers
645: and work space for performing finite difference approximations of
646: Jacobian-vector products, F'(u)*a,
648: The default code uses the following approach to compute h
650: .vb
651: F'(u)*a = [F(u+h*a) - F(u)]/h where
652: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
653: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
654: where
655: error_rel = square root of relative error in function evaluation
656: umin = minimum iterate parameter
657: .ve
659: You can call `SNESSetJacobian()` with `MatMFFDComputeJacobian()` if you are using matrix and not a different
660: preconditioner matrix
662: The user can set the error_rel via `MatMFFDSetFunctionError()` and
663: umin via `MatMFFDDSSetUmin()`; see Users-Manual: ch_snes for details.
665: The user should call `MatDestroy()` when finished with the matrix-free
666: matrix context.
668: .seealso: `MATMFFD`, `MatDestroy()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`, `MatMFFDSetFunction()`
669: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateSNESMF()`,
670: `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`
671: @*/
672: PetscErrorCode MatCreateMFFD(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, Mat *J)
673: {
674: MatCreate(comm, J);
675: MatSetSizes(*J, m, n, M, N);
676: MatSetType(*J, MATMFFD);
677: MatSetUp(*J);
678: return 0;
679: }
681: /*@
682: MatMFFDGetH - Gets the last value that was used as the differencing for a `MATMFFD` matrix
683: parameter.
685: Not Collective
687: Input Parameters:
688: . mat - the `MATMFFD` matrix
690: Output Parameter:
691: . h - the differencing step size
693: Level: advanced
695: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDSetHHistory()`, `MatCreateMFFD()`, `MATMFFD`, `MatMFFDResetHHistory()`
696: @*/
697: PetscErrorCode MatMFFDGetH(Mat mat, PetscScalar *h)
698: {
701: PetscUseMethod(mat, "MatMFFDGetH_C", (Mat, PetscScalar *), (mat, h));
702: return 0;
703: }
705: /*@C
706: MatMFFDSetFunction - Sets the function used in applying the matrix free `MATMFFD` matrix.
708: Logically Collective
710: Input Parameters:
711: + mat - the matrix free matrix `MATMFFD` created via `MatCreateSNESMF()` or `MatCreateMFFD()`
712: . func - the function to use
713: - funcctx - optional function context passed to function
715: Calling Sequence of func:
716: $ func (void *funcctx, Vec x, Vec f)
718: + funcctx - user provided context
719: . x - input vector
720: - f - computed output function
722: Level: advanced
724: Notes:
725: If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix free
726: matrix inside your compute Jacobian routine
728: If this is not set then it will use the function set with `SNESSetFunction()` if `MatCreateSNESMF()` was used.
730: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`, `MATMFFD`,
731: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`
732: @*/
733: PetscErrorCode MatMFFDSetFunction(Mat mat, PetscErrorCode (*func)(void *, Vec, Vec), void *funcctx)
734: {
736: PetscTryMethod(mat, "MatMFFDSetFunction_C", (Mat, PetscErrorCode(*)(void *, Vec, Vec), void *), (mat, func, funcctx));
737: return 0;
738: }
740: /*@C
741: MatMFFDSetFunctioni - Sets the function for a single component for a `MATMFFD` matrix
743: Logically Collective
745: Input Parameters:
746: + mat - the matrix free matrix `MATMFFD`
747: - funci - the function to use
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: This function is necessary to compute the diagonal of the matrix.
756: funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.
758: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`, `MatGetDiagonal()`
759: @*/
760: PetscErrorCode MatMFFDSetFunctioni(Mat mat, PetscErrorCode (*funci)(void *, PetscInt, Vec, PetscScalar *))
761: {
763: PetscTryMethod(mat, "MatMFFDSetFunctioni_C", (Mat, PetscErrorCode(*)(void *, PetscInt, Vec, PetscScalar *)), (mat, funci));
764: return 0;
765: }
767: /*@C
768: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation for a `MATMFFD` matrix
770: Logically Collective
772: Input Parameters:
773: + mat - the `MATMFFD` matrix free matrix
774: - func - the function to use
776: Level: advanced
778: Notes:
779: If you use this you MUST call `MatAssemblyBegin()` and `MatAssemblyEnd()` on the matrix free
780: matrix inside your compute Jacobian routine.
782: This function is necessary to compute the diagonal of the matrix, used for example with `PCJACOBI`
784: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`, `MATMFFD`
785: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `SNESetFunction()`, `MatGetDiagonal()`
786: @*/
787: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat, PetscErrorCode (*func)(void *, Vec))
788: {
790: PetscTryMethod(mat, "MatMFFDSetFunctioniBase_C", (Mat, PetscErrorCode(*)(void *, Vec)), (mat, func));
791: return 0;
792: }
794: /*@
795: MatMFFDSetPeriod - Sets how often h is recomputed for a `MATMFFD` matrix, by default it is every time
797: Logically Collective
799: Input Parameters:
800: + mat - the `MATMFFD` matrix free matrix
801: - period - 1 for every time, 2 for every second etc
803: Options Database Keys:
804: . -mat_mffd_period <period> - Sets how often h is recomputed
806: Level: advanced
808: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`,
809: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
810: @*/
811: PetscErrorCode MatMFFDSetPeriod(Mat mat, PetscInt period)
812: {
815: PetscTryMethod(mat, "MatMFFDSetPeriod_C", (Mat, PetscInt), (mat, period));
816: return 0;
817: }
819: /*@
820: MatMFFDSetFunctionError - Sets the error_rel for the approximation of matrix-vector products using finite differences with the `MATMFFD` matrix
822: Logically Collective
824: Input Parameters:
825: + mat - the `MATMFFD` matrix free matrix
826: - error_rel - relative error (should be set to the square root of the relative error in the function evaluations)
828: Options Database Keys:
829: . -mat_mffd_err <error_rel> - Sets error_rel
831: Level: advanced
833: Note:
834: The default matrix-free matrix-vector product routine computes
835: .vb
836: F'(u)*a = [F(u+h*a) - F(u)]/h where
837: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
838: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
839: .ve
841: .seealso: `MATMFFD`, `MatCreateSNESMF()`, `MatMFFDGetH()`, `MatCreateMFFD()`, `MATMFFD`
842: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`
843: @*/
844: PetscErrorCode MatMFFDSetFunctionError(Mat mat, PetscReal error)
845: {
848: PetscTryMethod(mat, "MatMFFDSetFunctionError_C", (Mat, PetscReal), (mat, error));
849: return 0;
850: }
852: /*@
853: MatMFFDSetHHistory - Sets an array to collect a history of the
854: differencing values (h) computed for the matrix-free product `MATMFFD` matrix
856: Logically Collective
858: Input Parameters:
859: + J - the `MATMFFD` matrix-free matrix
860: . history - space to hold the history
861: - nhistory - number of entries in history, if more entries are generated than
862: nhistory, then the later ones are discarded
864: Level: advanced
866: Note:
867: Use `MatMFFDResetHHistory()` to reset the history counter and collect
868: a new batch of differencing parameters, h.
870: .seealso: `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
871: `MatMFFDResetHHistory()`, `MatMFFDSetFunctionError()`
872: @*/
873: PetscErrorCode MatMFFDSetHHistory(Mat J, PetscScalar history[], PetscInt nhistory)
874: {
878: PetscUseMethod(J, "MatMFFDSetHHistory_C", (Mat, PetscScalar[], PetscInt), (J, history, nhistory));
879: return 0;
880: }
882: /*@
883: MatMFFDResetHHistory - Resets the counter to zero to begin
884: collecting a new set of differencing histories for the `MATMFFD` matrix
886: Logically Collective
888: Input Parameters:
889: . J - the matrix-free matrix context
891: Level: advanced
893: Note:
894: Use `MatMFFDSetHHistory()` to create the original history counter.
896: .seealso: `MATMFFD`, `MatMFFDGetH()`, `MatCreateSNESMF()`,
897: `MatMFFDSetHHistory()`, `MatMFFDSetFunctionError()`
898: @*/
899: PetscErrorCode MatMFFDResetHHistory(Mat J)
900: {
902: PetscTryMethod(J, "MatMFFDResetHHistory_C", (Mat), (J));
903: return 0;
904: }
906: /*@
907: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
908: Jacobian are computed for the `MATMFFD` matrix
910: Logically Collective
912: Input Parameters:
913: + J - the `MATMFFD` matrix
914: . U - the vector
915: - F - (optional) vector that contains F(u) if it has been already computed
917: Notes:
918: This is rarely used directly
920: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
921: point during the first `MatMult()` after each call to `MatMFFDSetBase()`.
923: Level: advanced
925: .seealso: `MATMFFD`, `MatMult()`, `MatMFFDSetBase()`
926: @*/
927: PetscErrorCode MatMFFDSetBase(Mat J, Vec U, Vec F)
928: {
932: PetscTryMethod(J, "MatMFFDSetBase_C", (Mat, Vec, Vec), (J, U, F));
933: return 0;
934: }
936: /*@C
937: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
938: it to satisfy some criteria for the `MATMFFD` matrix
940: Logically Collective
942: Input Parameters:
943: + J - the `MATMFFD` matrix
944: . fun - the function that checks h
945: - ctx - any context needed by the function
947: Options Database Keys:
948: . -mat_mffd_check_positivity <bool> - Insure that U + h*a is non-negative
950: Level: advanced
952: Notes:
953: For example, `MatMFFDCheckPositivity()` insures that all entries
954: of U + h*a are non-negative
956: The function you provide is called after the default h has been computed and allows you to
957: modify it.
959: .seealso: `MATMFFD`, `MatMFFDCheckPositivity()`
960: @*/
961: PetscErrorCode MatMFFDSetCheckh(Mat J, PetscErrorCode (*fun)(void *, Vec, Vec, PetscScalar *), void *ctx)
962: {
964: PetscTryMethod(J, "MatMFFDSetCheckh_C", (Mat, PetscErrorCode(*)(void *, Vec, Vec, PetscScalar *), void *), (J, fun, ctx));
965: return 0;
966: }
968: /*@
969: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
970: zero, decreases h until this is satisfied for a `MATMFFD` matrix
972: Logically Collective
974: Input Parameters:
975: + U - base vector that is added to
976: . a - vector that is added
977: . h - scaling factor on a
978: - dummy - context variable (unused)
980: Options Database Keys:
981: . -mat_mffd_check_positivity <bool> - Insure that U + h*a is nonnegative
983: Level: advanced
985: Note:
986: This is rarely used directly, rather it is passed as an argument to
987: `MatMFFDSetCheckh()`
989: .seealso: `MATMFFD`, `MatMFFDSetCheckh()`
990: @*/
991: PetscErrorCode MatMFFDCheckPositivity(void *dummy, Vec U, Vec a, PetscScalar *h)
992: {
993: PetscReal val, minval;
994: PetscScalar *u_vec, *a_vec;
995: PetscInt i, n;
996: MPI_Comm comm;
1001: PetscObjectGetComm((PetscObject)U, &comm);
1002: VecGetArray(U, &u_vec);
1003: VecGetArray(a, &a_vec);
1004: VecGetLocalSize(U, &n);
1005: minval = PetscAbsScalar(*h) * PetscRealConstant(1.01);
1006: for (i = 0; i < n; i++) {
1007: if (PetscRealPart(u_vec[i] + *h * a_vec[i]) <= 0.0) {
1008: val = PetscAbsScalar(u_vec[i] / a_vec[i]);
1009: if (val < minval) minval = val;
1010: }
1011: }
1012: VecRestoreArray(U, &u_vec);
1013: VecRestoreArray(a, &a_vec);
1014: MPIU_Allreduce(&minval, &val, 1, MPIU_REAL, MPIU_MIN, comm);
1015: if (val <= PetscAbsScalar(*h)) {
1016: PetscInfo(U, "Scaling back h from %g to %g\n", (double)PetscRealPart(*h), (double)(.99 * val));
1017: if (PetscRealPart(*h) > 0.0) *h = 0.99 * val;
1018: else *h = -0.99 * val;
1019: }
1020: return 0;
1021: }