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: 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: 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: if (ctx->ops->destroy) (*ctx->ops->destroy)(ctx);
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 differene 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: Notes:
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: 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: Notes:
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: 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) {
203: VecDestroy(&ctx->current_f);
204: }
205: if (ctx->ops->destroy) (*ctx->ops->destroy)(ctx);
206: PetscHeaderDestroy(&ctx);
208: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
209: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
210: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
211: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
212: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
213: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
214: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
215: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
216: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetHHistory_C",NULL);
217: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetType_C",NULL);
218: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDGetH_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) {
245: PetscViewerASCIIPrintf(viewer,"Using Lyness complex number trick to compute the matrix-vector product\n");
246: }
247: #endif
248: if (ctx->ops->view) {
249: (*ctx->ops->view)(ctx,viewer);
250: }
251: PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);
253: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_base", &viewbase);
254: if (viewbase) {
255: PetscViewerASCIIPrintf(viewer, "Base:\n");
256: VecView(ctx->current_u, viewer);
257: }
258: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_function", &viewfunction);
259: if (viewfunction) {
260: PetscViewerASCIIPrintf(viewer, "Function:\n");
261: VecView(ctx->current_f, viewer);
262: }
263: PetscViewerASCIIPopTab(viewer);
264: }
265: return 0;
266: }
268: /*
269: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
270: allows the user to indicate the beginning of a new linear solve by calling
271: MatAssemblyXXX() on the matrix free matrix. This then allows the
272: MatCreateMFFD_WP() to properly compute ||U|| only the first time
273: in the linear solver rather than every time.
275: This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
276: it must be labeled as PETSC_EXTERN
277: */
278: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
279: {
280: MatMFFD j;
282: MatShellGetContext(J,&j);
283: MatMFFDResetHHistory(J);
284: return 0;
285: }
287: /*
288: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
290: y ~= (F(u + ha) - F(u))/h,
291: where F = nonlinear function, as set by SNESSetFunction()
292: u = current iterate
293: h = difference interval
294: */
295: static PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
296: {
297: MatMFFD ctx;
298: PetscScalar h;
299: Vec w,U,F;
300: PetscBool zeroa;
302: MatShellGetContext(mat,&ctx);
304: /* We log matrix-free matrix-vector products separately, so that we can
305: separate the performance monitoring from the cases that use conventional
306: storage. We may eventually modify event logging to associate events
307: with particular objects, hence alleviating the more general problem. */
308: PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);
310: w = ctx->w;
311: U = ctx->current_u;
312: F = ctx->current_f;
313: /*
314: Compute differencing parameter
315: */
316: if (!((PetscObject)ctx)->type_name) {
317: MatMFFDSetType(mat,MATMFFD_WP);
318: MatSetFromOptions(mat);
319: }
320: (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
321: if (zeroa) {
322: VecSet(y,0.0);
323: return 0;
324: }
327: if (ctx->checkh) {
328: (*ctx->checkh)(ctx->checkhctx,U,a,&h);
329: }
331: /* keep a record of the current differencing parameter h */
332: ctx->currenth = h;
333: #if defined(PETSC_USE_COMPLEX)
334: PetscInfo(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
335: #else
336: PetscInfo(mat,"Current differencing parameter: %15.12e\n",(double)PetscRealPart(h));
337: #endif
338: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
339: ctx->historyh[ctx->ncurrenth] = h;
340: }
341: ctx->ncurrenth++;
343: #if defined(PETSC_USE_COMPLEX)
344: if (ctx->usecomplex) h = PETSC_i*h;
345: #endif
347: /* w = u + ha */
348: VecWAXPY(w,h,a,U);
350: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
351: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
352: (*ctx->func)(ctx->funcctx,U,F);
353: }
354: (*ctx->func)(ctx->funcctx,w,y);
356: #if defined(PETSC_USE_COMPLEX)
357: if (ctx->usecomplex) {
358: VecImaginaryPart(y);
359: h = PetscImaginaryPart(h);
360: } else {
361: VecAXPY(y,-1.0,F);
362: }
363: #else
364: VecAXPY(y,-1.0,F);
365: #endif
366: VecScale(y,1.0/h);
367: if (mat->nullsp) MatNullSpaceRemove(mat->nullsp,y);
369: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
370: return 0;
371: }
373: /*
374: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
376: y ~= (F(u + ha) - F(u))/h,
377: where F = nonlinear function, as set by SNESSetFunction()
378: u = current iterate
379: h = difference interval
380: */
381: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
382: {
383: MatMFFD ctx;
384: PetscScalar h,*aa,*ww,v;
385: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
386: Vec w,U;
387: PetscInt i,rstart,rend;
389: MatShellGetContext(mat,&ctx);
392: w = ctx->w;
393: U = ctx->current_u;
394: (*ctx->func)(ctx->funcctx,U,a);
395: if (ctx->funcisetbase) {
396: (*ctx->funcisetbase)(ctx->funcctx,U);
397: }
398: VecCopy(U,w);
400: VecGetOwnershipRange(a,&rstart,&rend);
401: VecGetArray(a,&aa);
402: for (i=rstart; i<rend; i++) {
403: 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: VecRestoreArray(w,&ww);
412: (*ctx->funci)(ctx->funcctx,i,w,&v);
413: aa[i-rstart] = (v - aa[i-rstart])/h;
415: VecGetArray(w,&ww);
416: ww[i-rstart] -= h;
417: VecRestoreArray(w,&ww);
418: }
419: VecRestoreArray(a,&aa);
420: return 0;
421: }
423: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
424: {
425: MatMFFD ctx;
427: MatShellGetContext(J,&ctx);
428: MatMFFDResetHHistory(J);
429: if (!ctx->current_u) {
430: VecDuplicate(U,&ctx->current_u);
431: VecLockReadPush(ctx->current_u);
432: }
433: VecLockReadPop(ctx->current_u);
434: VecCopy(U,ctx->current_u);
435: VecLockReadPush(ctx->current_u);
436: if (F) {
437: if (ctx->current_f_allocated) VecDestroy(&ctx->current_f);
438: ctx->current_f = F;
439: ctx->current_f_allocated = PETSC_FALSE;
440: } else if (!ctx->current_f_allocated) {
441: MatCreateVecs(J,NULL,&ctx->current_f);
442: ctx->current_f_allocated = PETSC_TRUE;
443: }
444: if (!ctx->w) {
445: VecDuplicate(ctx->current_u,&ctx->w);
446: }
447: J->assembled = PETSC_TRUE;
448: return 0;
449: }
451: typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
453: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
454: {
455: MatMFFD ctx;
457: MatShellGetContext(J,&ctx);
458: ctx->checkh = fun;
459: ctx->checkhctx = ectx;
460: return 0;
461: }
463: /*@C
464: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
465: MatMFFD options in the database.
467: Collective on Mat
469: Input Parameters:
470: + A - the Mat context
471: - prefix - the prefix to prepend to all option names
473: Notes:
474: A hyphen (-) must NOT be given at the beginning of the prefix name.
475: The first character of all runtime options is AUTOMATICALLY the hyphen.
477: Level: advanced
479: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
480: @*/
481: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
482: {
483: MatMFFD mfctx;
486: MatShellGetContext(mat,&mfctx);
488: PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
489: return 0;
490: }
492: static PetscErrorCode MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
493: {
494: MatMFFD mfctx;
496: PetscBool flg;
497: char ftype[256];
500: MatShellGetContext(mat,&mfctx);
502: PetscObjectOptionsBegin((PetscObject)mfctx);
503: PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
504: if (flg) {
505: MatMFFDSetType(mat,ftype);
506: }
508: PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,NULL);
509: PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,NULL);
511: flg = PETSC_FALSE;
512: PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
513: if (flg) {
514: MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,NULL);
515: }
516: #if defined(PETSC_USE_COMPLEX)
517: PetscOptionsBool("-mat_mffd_complex","Use Lyness complex number trick to compute the matrix-vector product","None",mfctx->usecomplex,&mfctx->usecomplex,NULL);
518: #endif
519: if (mfctx->ops->setfromoptions) {
520: (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
521: }
522: PetscOptionsEnd();
523: return 0;
524: }
526: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
527: {
528: MatMFFD ctx;
530: MatShellGetContext(mat,&ctx);
531: ctx->recomputeperiod = period;
532: return 0;
533: }
535: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
536: {
537: MatMFFD ctx;
539: MatShellGetContext(mat,&ctx);
540: ctx->func = func;
541: ctx->funcctx = funcctx;
542: return 0;
543: }
545: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
546: {
547: MatMFFD ctx;
549: MatShellGetContext(mat,&ctx);
550: if (error != PETSC_DEFAULT) ctx->error_rel = error;
551: return 0;
552: }
554: PetscErrorCode MatMFFDSetHHistory_MFFD(Mat J,PetscScalar history[],PetscInt nhistory)
555: {
556: MatMFFD ctx;
558: MatShellGetContext(J,&ctx);
559: ctx->historyh = history;
560: ctx->maxcurrenth = nhistory;
561: ctx->currenth = 0.;
562: return 0;
563: }
565: /*MC
566: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
568: Level: advanced
570: Developers Note: This is implemented on top of MATSHELL to get support for scaling and shifting without requiring duplicate code
572: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),
573: MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
574: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
575: MatMFFDGetH(),
576: M*/
577: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
578: {
579: MatMFFD mfctx;
581: MatMFFDInitializePackage();
583: PetscHeaderCreate(mfctx,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),NULL,NULL);
585: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
586: mfctx->recomputeperiod = 1;
587: mfctx->count = 0;
588: mfctx->currenth = 0.0;
589: mfctx->historyh = NULL;
590: mfctx->ncurrenth = 0;
591: mfctx->maxcurrenth = 0;
592: ((PetscObject)mfctx)->type_name = NULL;
594: /*
595: Create the empty data structure to contain compute-h routines.
596: These will be filled in below from the command line options or
597: a later call with MatMFFDSetType() or if that is not called
598: then it will default in the first use of MatMult_MFFD()
599: */
600: mfctx->ops->compute = NULL;
601: mfctx->ops->destroy = NULL;
602: mfctx->ops->view = NULL;
603: mfctx->ops->setfromoptions = NULL;
604: mfctx->hctx = NULL;
606: mfctx->func = NULL;
607: mfctx->funcctx = NULL;
608: mfctx->w = NULL;
609: mfctx->mat = A;
611: MatSetType(A,MATSHELL);
612: MatShellSetContext(A,mfctx);
613: MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_MFFD);
614: MatShellSetOperation(A,MATOP_DESTROY,(void (*)(void))MatDestroy_MFFD);
615: MatShellSetOperation(A,MATOP_VIEW,(void (*)(void))MatView_MFFD);
616: MatShellSetOperation(A,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MFFD);
617: MatShellSetOperation(A,MATOP_SET_FROM_OPTIONS,(void (*)(void))MatSetFromOptions_MFFD);
619: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
620: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
621: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
622: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
623: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
624: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
625: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
626: PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
627: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetHHistory_C",MatMFFDSetHHistory_MFFD);
628: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetType_C",MatMFFDSetType_MFFD);
629: PetscObjectComposeFunction((PetscObject)A,"MatMFFDGetH_C",MatMFFDGetH_MFFD);
630: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
631: return 0;
632: }
634: /*@
635: MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
637: Collective on Vec
639: Input Parameters:
640: + comm - MPI communicator
641: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
642: This value should be the same as the local size used in creating the
643: y vector for the matrix-vector product y = Ax.
644: . n - This value should be the same as the local size used in creating the
645: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
646: calculated if N is given) For square matrices n is almost always m.
647: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
648: - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
650: Output Parameter:
651: . J - the matrix-free matrix
653: Options Database Keys: call MatSetFromOptions() to trigger these
654: + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
655: . -mat_mffd_err - square root of estimated relative error in function evaluation
656: . -mat_mffd_period - how often h is recomputed, defaults to 1, every time
657: . -mat_mffd_check_positivity - possibly decrease h until U + h*a has only positive values
658: - -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
659: (requires real valued functions but that PETSc be configured for complex numbers)
661: Level: advanced
663: Notes:
664: The matrix-free matrix context merely contains the function pointers
665: and work space for performing finite difference approximations of
666: Jacobian-vector products, F'(u)*a,
668: The default code uses the following approach to compute h
670: .vb
671: F'(u)*a = [F(u+h*a) - F(u)]/h where
672: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
673: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
674: where
675: error_rel = square root of relative error in function evaluation
676: umin = minimum iterate parameter
677: .ve
679: You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
680: preconditioner matrix
682: The user can set the error_rel via MatMFFDSetFunctionError() and
683: umin via MatMFFDDSSetUmin(); see Users-Manual: ch_snes for details.
685: The user should call MatDestroy() when finished with the matrix-free
686: matrix context.
688: Options Database Keys:
689: + -mat_mffd_err <error_rel> - Sets error_rel
690: . -mat_mffd_umin <umin> - Sets umin (for default PETSc routine that computes h only)
691: - -mat_mffd_check_positivity - check positivity
693: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
694: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
695: MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
697: @*/
698: PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
699: {
700: MatCreate(comm,J);
701: MatSetSizes(*J,m,n,M,N);
702: MatSetType(*J,MATMFFD);
703: MatSetUp(*J);
704: return 0;
705: }
707: /*@
708: MatMFFDGetH - Gets the last value that was used as the differencing
709: parameter.
711: Not Collective
713: Input Parameters:
714: . mat - the matrix obtained with MatCreateSNESMF()
716: Output Parameter:
717: . h - the differencing step size
719: Level: advanced
721: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
722: @*/
723: PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
724: {
727: PetscUseMethod(mat,"MatMFFDGetH_C",(Mat,PetscScalar*),(mat,h));
728: return 0;
729: }
731: /*@C
732: MatMFFDSetFunction - Sets the function used in applying the matrix free.
734: Logically Collective on Mat
736: Input Parameters:
737: + mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
738: . func - the function to use
739: - funcctx - optional function context passed to function
741: Calling Sequence of func:
742: $ func (void *funcctx, Vec x, Vec f)
744: + funcctx - user provided context
745: . x - input vector
746: - f - computed output function
748: Level: advanced
750: Notes:
751: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
752: matrix inside your compute Jacobian routine
754: If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
756: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
757: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
758: @*/
759: PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
760: {
762: PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
763: return 0;
764: }
766: /*@C
767: MatMFFDSetFunctioni - Sets the function for a single component
769: Logically Collective on Mat
771: Input Parameters:
772: + mat - the matrix free matrix created via MatCreateSNESMF()
773: - funci - the function to use
775: Level: advanced
777: Notes:
778: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
779: matrix inside your compute Jacobian routine.
780: This function is necessary to compute the diagonal of the matrix.
781: funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.
783: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
785: @*/
786: PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
787: {
789: PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
790: return 0;
791: }
793: /*@C
794: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
796: Logically Collective on Mat
798: Input Parameters:
799: + mat - the matrix free matrix created via MatCreateSNESMF()
800: - func - the function to use
802: Level: advanced
804: Notes:
805: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
806: matrix inside your compute Jacobian routine.
807: This function is necessary to compute the diagonal of the matrix.
809: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
810: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
811: @*/
812: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
813: {
815: PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
816: return 0;
817: }
819: /*@
820: MatMFFDSetPeriod - Sets how often h is recomputed, by default it is every time
822: Logically Collective on Mat
824: Input Parameters:
825: + mat - the matrix free matrix created via MatCreateSNESMF()
826: - period - 1 for every time, 2 for every second etc
828: Options Database Keys:
829: . -mat_mffd_period <period> - Sets how often h is recomputed
831: Level: advanced
833: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
834: MatMFFDSetHHistory(), MatMFFDResetHHistory()
835: @*/
836: PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
837: {
840: PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
841: return 0;
842: }
844: /*@
845: MatMFFDSetFunctionError - Sets the error_rel for the approximation of
846: matrix-vector products using finite differences.
848: Logically Collective on Mat
850: Input Parameters:
851: + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
852: - error_rel - relative error (should be set to the square root of
853: the relative error in the function evaluations)
855: Options Database Keys:
856: . -mat_mffd_err <error_rel> - Sets error_rel
858: Level: advanced
860: Notes:
861: The default matrix-free matrix-vector product routine computes
862: .vb
863: F'(u)*a = [F(u+h*a) - F(u)]/h where
864: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
865: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
866: .ve
868: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
869: MatMFFDSetHHistory(), MatMFFDResetHHistory()
870: @*/
871: PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
872: {
875: PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
876: return 0;
877: }
879: /*@
880: MatMFFDSetHHistory - Sets an array to collect a history of the
881: differencing values (h) computed for the matrix-free product.
883: Logically Collective on Mat
885: Input Parameters:
886: + J - the matrix-free matrix context
887: . history - space to hold the history
888: - nhistory - number of entries in history, if more entries are generated than
889: nhistory, then the later ones are discarded
891: Level: advanced
893: Notes:
894: Use MatMFFDResetHHistory() to reset the history counter and collect
895: a new batch of differencing parameters, h.
897: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
898: MatMFFDResetHHistory(), MatMFFDSetFunctionError()
900: @*/
901: PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
902: {
906: PetscUseMethod(J,"MatMFFDSetHHistory_C",(Mat,PetscScalar[],PetscInt),(J,history,nhistory));
907: return 0;
908: }
910: /*@
911: MatMFFDResetHHistory - Resets the counter to zero to begin
912: collecting a new set of differencing histories.
914: Logically Collective on Mat
916: Input Parameters:
917: . J - the matrix-free matrix context
919: Level: advanced
921: Notes:
922: Use MatMFFDSetHHistory() to create the original history counter.
924: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
925: MatMFFDSetHHistory(), MatMFFDSetFunctionError()
927: @*/
928: PetscErrorCode MatMFFDResetHHistory(Mat J)
929: {
931: PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
932: return 0;
933: }
935: /*@
936: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
937: Jacobian are computed
939: Logically Collective on Mat
941: Input Parameters:
942: + J - the MatMFFD matrix
943: . U - the vector
944: - F - (optional) vector that contains F(u) if it has been already computed
946: Notes:
947: This is rarely used directly
949: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
950: point during the first MatMult() after each call to MatMFFDSetBase().
952: Level: advanced
954: @*/
955: PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
956: {
960: PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
961: return 0;
962: }
964: /*@C
965: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
966: it to satisfy some criteria
968: Logically Collective on Mat
970: Input Parameters:
971: + J - the MatMFFD matrix
972: . fun - the function that checks h
973: - ctx - any context needed by the function
975: Options Database Keys:
976: . -mat_mffd_check_positivity <bool> - Insure that U + h*a is non-negative
978: Level: advanced
980: Notes:
981: For example, MatMFFDCheckPositivity() insures that all entries
982: of U + h*a are non-negative
984: The function you provide is called after the default h has been computed and allows you to
985: modify it.
987: .seealso: MatMFFDCheckPositivity()
988: @*/
989: PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
990: {
992: PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
993: return 0;
994: }
996: /*@
997: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
998: zero, decreases h until this is satisfied.
1000: Logically Collective on Vec
1002: Input Parameters:
1003: + U - base vector that is added to
1004: . a - vector that is added
1005: . h - scaling factor on a
1006: - dummy - context variable (unused)
1008: Options Database Keys:
1009: . -mat_mffd_check_positivity <bool> - Insure that U + h*a is nonnegative
1011: Level: advanced
1013: Notes:
1014: This is rarely used directly, rather it is passed as an argument to
1015: MatMFFDSetCheckh()
1017: .seealso: MatMFFDSetCheckh()
1018: @*/
1019: PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1020: {
1021: PetscReal val, minval;
1022: PetscScalar *u_vec, *a_vec;
1023: PetscInt i,n;
1024: MPI_Comm comm;
1029: PetscObjectGetComm((PetscObject)U,&comm);
1030: VecGetArray(U,&u_vec);
1031: VecGetArray(a,&a_vec);
1032: VecGetLocalSize(U,&n);
1033: minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1034: for (i=0; i<n; i++) {
1035: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1036: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1037: if (val < minval) minval = val;
1038: }
1039: }
1040: VecRestoreArray(U,&u_vec);
1041: VecRestoreArray(a,&a_vec);
1042: MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1043: if (val <= PetscAbsScalar(*h)) {
1044: PetscInfo(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1045: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1046: else *h = -0.99*val;
1047: }
1048: return 0;
1049: }