Actual source code: mffd.c
petsc-3.14.6 2021-03-30
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: {
25: PetscFunctionListDestroy(&MatMFFDList);
26: MatMFFDPackageInitialized = PETSC_FALSE;
27: MatMFFDRegisterAllCalled = PETSC_FALSE;
28: return(0);
29: }
31: /*@C
32: MatMFFDInitializePackage - This function initializes everything in the MatMFFD package. It is called
33: from MatInitializePackage().
35: Level: developer
37: .seealso: PetscInitialize()
38: @*/
39: PetscErrorCode MatMFFDInitializePackage(void)
40: {
41: char logList[256];
42: PetscBool opt,pkg;
46: if (MatMFFDPackageInitialized) return(0);
47: MatMFFDPackageInitialized = PETSC_TRUE;
48: /* Register Classes */
49: PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);
50: /* Register Constructors */
51: MatMFFDRegisterAll();
52: /* Register Events */
53: PetscLogEventRegister("MatMult MF",MATMFFD_CLASSID,&MATMFFD_Mult);
54: /* Process Info */
55: {
56: PetscClassId classids[1];
58: classids[0] = MATMFFD_CLASSID;
59: PetscInfoProcessClass("matmffd", 1, classids);
60: }
61: /* Process summary exclusions */
62: PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
63: if (opt) {
64: PetscStrInList("matmffd",logList,',',&pkg);
65: if (pkg) {PetscLogEventExcludeClass(MATMFFD_CLASSID);}
66: }
67: /* Register package finalizer */
68: PetscRegisterFinalize(MatMFFDFinalizePackage);
69: return(0);
70: }
72: static PetscErrorCode MatMFFDSetType_MFFD(Mat mat,MatMFFDType ftype)
73: {
74: PetscErrorCode ierr,(*r)(MatMFFD);
75: MatMFFD ctx;
76: PetscBool match;
81: MatShellGetContext(mat,&ctx);
83: /* already set, so just return */
84: PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
85: if (match) return(0);
87: /* destroy the old one if it exists */
88: if (ctx->ops->destroy) {
89: (*ctx->ops->destroy)(ctx);
90: }
92: PetscFunctionListFind(MatMFFDList,ftype,&r);
93: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
94: (*r)(ctx);
95: PetscObjectChangeTypeName((PetscObject)ctx,ftype);
96: return(0);
97: }
99: /*@C
100: MatMFFDSetType - Sets the method that is used to compute the
101: differencing parameter for finite differene matrix-free formulations.
103: Input Parameters:
104: + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
105: or MatSetType(mat,MATMFFD);
106: - ftype - the type requested, either MATMFFD_WP or MATMFFD_DS
108: Level: advanced
110: Notes:
111: For example, such routines can compute h for use in
112: Jacobian-vector products of the form
114: F(x+ha) - F(x)
115: F'(u)a ~= ----------------
116: h
118: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction(), MatCreateMFFD()
119: @*/
120: PetscErrorCode MatMFFDSetType(Mat mat,MatMFFDType ftype)
121: {
127: PetscTryMethod(mat,"MatMFFDSetType_C",(Mat,MatMFFDType),(mat,ftype));
128: return(0);
129: }
131: static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);
133: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
134: static PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
135: {
136: MatMFFD ctx;
140: MatShellGetContext(mat,&ctx);
141: ctx->funcisetbase = func;
142: return(0);
143: }
145: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
146: static PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
147: {
148: MatMFFD ctx;
152: MatShellGetContext(mat,&ctx);
153: ctx->funci = funci;
154: MatShellSetOperation(mat,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_MFFD);
155: return(0);
156: }
158: static PetscErrorCode MatMFFDGetH_MFFD(Mat mat,PetscScalar *h)
159: {
160: MatMFFD ctx;
164: MatShellGetContext(mat,&ctx);
165: *h = ctx->currenth;
166: return(0);
167: }
169: static PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
170: {
171: MatMFFD ctx;
175: MatShellGetContext(J,&ctx);
176: ctx->ncurrenth = 0;
177: return(0);
178: }
180: /*@C
181: MatMFFDRegister - Adds a method to the MatMFFD registry.
183: Not Collective
185: Input Parameters:
186: + name_solver - name of a new user-defined compute-h module
187: - routine_create - routine to create method context
189: Level: developer
191: Notes:
192: MatMFFDRegister() may be called multiple times to add several user-defined solvers.
194: Sample usage:
195: .vb
196: MatMFFDRegister("my_h",MyHCreate);
197: .ve
199: Then, your solver can be chosen with the procedural interface via
200: $ MatMFFDSetType(mfctx,"my_h")
201: or at runtime via the option
202: $ -mat_mffd_type my_h
204: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
205: @*/
206: PetscErrorCode MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
207: {
211: MatInitializePackage();
212: PetscFunctionListAdd(&MatMFFDList,sname,function);
213: return(0);
214: }
216: /* ----------------------------------------------------------------------------------------*/
217: static PetscErrorCode MatDestroy_MFFD(Mat mat)
218: {
220: MatMFFD ctx;
223: MatShellGetContext(mat,&ctx);
224: VecDestroy(&ctx->w);
225: VecDestroy(&ctx->current_u);
226: if (ctx->current_f_allocated) {
227: VecDestroy(&ctx->current_f);
228: }
229: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
230: PetscHeaderDestroy(&ctx);
232: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
233: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
234: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
235: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
236: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
237: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
238: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
239: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
240: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetHHistory_C",NULL);
241: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetType_C",NULL);
242: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDGetH_C",NULL);
243: return(0);
244: }
246: /*
247: MatMFFDView_MFFD - Views matrix-free parameters.
249: */
250: static PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
251: {
253: MatMFFD ctx;
254: PetscBool iascii, viewbase, viewfunction;
255: const char *prefix;
258: MatShellGetContext(J,&ctx);
259: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
260: if (iascii) {
261: PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
262: PetscViewerASCIIPushTab(viewer);
263: PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
264: if (!((PetscObject)ctx)->type_name) {
265: PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
266: } else {
267: PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
268: }
269: #if defined(PETSC_USE_COMPLEX)
270: if (ctx->usecomplex) {
271: PetscViewerASCIIPrintf(viewer,"Using Lyness complex number trick to compute the matrix-vector product\n");
272: }
273: #endif
274: if (ctx->ops->view) {
275: (*ctx->ops->view)(ctx,viewer);
276: }
277: PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);
279: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_base", &viewbase);
280: if (viewbase) {
281: PetscViewerASCIIPrintf(viewer, "Base:\n");
282: VecView(ctx->current_u, viewer);
283: }
284: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_function", &viewfunction);
285: if (viewfunction) {
286: PetscViewerASCIIPrintf(viewer, "Function:\n");
287: VecView(ctx->current_f, viewer);
288: }
289: PetscViewerASCIIPopTab(viewer);
290: }
291: return(0);
292: }
294: /*
295: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
296: allows the user to indicate the beginning of a new linear solve by calling
297: MatAssemblyXXX() on the matrix free matrix. This then allows the
298: MatCreateMFFD_WP() to properly compute ||U|| only the first time
299: in the linear solver rather than every time.
301: This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
302: it must be labeled as PETSC_EXTERN
303: */
304: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
305: {
307: MatMFFD j;
310: MatShellGetContext(J,&j);
311: MatMFFDResetHHistory(J);
312: return(0);
313: }
315: /*
316: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
318: y ~= (F(u + ha) - F(u))/h,
319: where F = nonlinear function, as set by SNESSetFunction()
320: u = current iterate
321: h = difference interval
322: */
323: static PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
324: {
325: MatMFFD ctx;
326: PetscScalar h;
327: Vec w,U,F;
329: PetscBool zeroa;
332: MatShellGetContext(mat,&ctx);
333: if (!ctx->current_u) SETERRQ(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");
334: /* We log matrix-free matrix-vector products separately, so that we can
335: separate the performance monitoring from the cases that use conventional
336: storage. We may eventually modify event logging to associate events
337: with particular objects, hence alleviating the more general problem. */
338: PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);
340: w = ctx->w;
341: U = ctx->current_u;
342: F = ctx->current_f;
343: /*
344: Compute differencing parameter
345: */
346: if (!((PetscObject)ctx)->type_name) {
347: MatMFFDSetType(mat,MATMFFD_WP);
348: MatSetFromOptions(mat);
349: }
350: (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
351: if (zeroa) {
352: VecSet(y,0.0);
353: return(0);
354: }
356: if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
357: if (ctx->checkh) {
358: (*ctx->checkh)(ctx->checkhctx,U,a,&h);
359: }
361: /* keep a record of the current differencing parameter h */
362: ctx->currenth = h;
363: #if defined(PETSC_USE_COMPLEX)
364: PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
365: #else
366: PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
367: #endif
368: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
369: ctx->historyh[ctx->ncurrenth] = h;
370: }
371: ctx->ncurrenth++;
373: #if defined(PETSC_USE_COMPLEX)
374: if (ctx->usecomplex) h = PETSC_i*h;
375: #endif
377: /* w = u + ha */
378: VecWAXPY(w,h,a,U);
380: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
381: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
382: (*ctx->func)(ctx->funcctx,U,F);
383: }
384: (*ctx->func)(ctx->funcctx,w,y);
386: #if defined(PETSC_USE_COMPLEX)
387: if (ctx->usecomplex) {
388: VecImaginaryPart(y);
389: h = PetscImaginaryPart(h);
390: } else {
391: VecAXPY(y,-1.0,F);
392: }
393: #else
394: VecAXPY(y,-1.0,F);
395: #endif
396: VecScale(y,1.0/h);
397: if (mat->nullsp) {MatNullSpaceRemove(mat->nullsp,y);}
399: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
400: return(0);
401: }
403: /*
404: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
406: y ~= (F(u + ha) - F(u))/h,
407: where F = nonlinear function, as set by SNESSetFunction()
408: u = current iterate
409: h = difference interval
410: */
411: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
412: {
413: MatMFFD ctx;
414: PetscScalar h,*aa,*ww,v;
415: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
416: Vec w,U;
418: PetscInt i,rstart,rend;
421: MatShellGetContext(mat,&ctx);
422: if (!ctx->func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunction() first");
423: if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioni() first");
424: w = ctx->w;
425: U = ctx->current_u;
426: (*ctx->func)(ctx->funcctx,U,a);
427: if (ctx->funcisetbase) {
428: (*ctx->funcisetbase)(ctx->funcctx,U);
429: }
430: VecCopy(U,w);
432: VecGetOwnershipRange(a,&rstart,&rend);
433: VecGetArray(a,&aa);
434: for (i=rstart; i<rend; i++) {
435: VecGetArray(w,&ww);
436: h = ww[i-rstart];
437: if (h == 0.0) h = 1.0;
438: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
439: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
440: h *= epsilon;
442: ww[i-rstart] += h;
443: VecRestoreArray(w,&ww);
444: (*ctx->funci)(ctx->funcctx,i,w,&v);
445: aa[i-rstart] = (v - aa[i-rstart])/h;
447: VecGetArray(w,&ww);
448: ww[i-rstart] -= h;
449: VecRestoreArray(w,&ww);
450: }
451: VecRestoreArray(a,&aa);
452: return(0);
453: }
455: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
456: {
458: MatMFFD ctx;
461: MatShellGetContext(J,&ctx);
462: MatMFFDResetHHistory(J);
463: if (!ctx->current_u) {
464: VecDuplicate(U,&ctx->current_u);
465: VecLockReadPush(ctx->current_u);
466: }
467: VecLockReadPop(ctx->current_u);
468: VecCopy(U,ctx->current_u);
469: VecLockReadPush(ctx->current_u);
470: if (F) {
471: if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
472: ctx->current_f = F;
473: ctx->current_f_allocated = PETSC_FALSE;
474: } else if (!ctx->current_f_allocated) {
475: MatCreateVecs(J,NULL,&ctx->current_f);
477: ctx->current_f_allocated = PETSC_TRUE;
478: }
479: if (!ctx->w) {
480: VecDuplicate(ctx->current_u,&ctx->w);
481: }
482: J->assembled = PETSC_TRUE;
483: return(0);
484: }
486: typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
488: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
489: {
490: MatMFFD ctx;
494: MatShellGetContext(J,&ctx);
495: ctx->checkh = fun;
496: ctx->checkhctx = ectx;
497: return(0);
498: }
500: /*@C
501: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
502: MatMFFD options in the database.
504: Collective on Mat
506: Input Parameter:
507: + A - the Mat context
508: - prefix - the prefix to prepend to all option names
510: Notes:
511: A hyphen (-) must NOT be given at the beginning of the prefix name.
512: The first character of all runtime options is AUTOMATICALLY the hyphen.
514: Level: advanced
516: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
517: @*/
518: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
520: {
521: MatMFFD mfctx;
526: MatShellGetContext(mat,&mfctx);
528: PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
529: return(0);
530: }
532: static PetscErrorCode MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
533: {
534: MatMFFD mfctx;
536: PetscBool flg;
537: char ftype[256];
541: MatShellGetContext(mat,&mfctx);
543: PetscObjectOptionsBegin((PetscObject)mfctx);
544: PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
545: if (flg) {
546: MatMFFDSetType(mat,ftype);
547: }
549: PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,NULL);
550: PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,NULL);
552: flg = PETSC_FALSE;
553: PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
554: if (flg) {
555: MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,NULL);
556: }
557: #if defined(PETSC_USE_COMPLEX)
558: PetscOptionsBool("-mat_mffd_complex","Use Lyness complex number trick to compute the matrix-vector product","None",mfctx->usecomplex,&mfctx->usecomplex,NULL);
559: #endif
560: if (mfctx->ops->setfromoptions) {
561: (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
562: }
563: PetscOptionsEnd();
564: return(0);
565: }
567: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
568: {
569: MatMFFD ctx;
573: MatShellGetContext(mat,&ctx);
574: ctx->recomputeperiod = period;
575: return(0);
576: }
578: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
579: {
580: MatMFFD ctx;
584: MatShellGetContext(mat,&ctx);
585: ctx->func = func;
586: ctx->funcctx = funcctx;
587: return(0);
588: }
590: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
591: {
592: MatMFFD ctx;
596: MatShellGetContext(mat,&ctx);
597: if (error != PETSC_DEFAULT) ctx->error_rel = error;
598: return(0);
599: }
601: PetscErrorCode MatMFFDSetHHistory_MFFD(Mat J,PetscScalar history[],PetscInt nhistory)
602: {
603: MatMFFD ctx;
607: MatShellGetContext(J,&ctx);
608: ctx->historyh = history;
609: ctx->maxcurrenth = nhistory;
610: ctx->currenth = 0.;
611: return(0);
612: }
614: /*MC
615: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
617: Level: advanced
619: Developers Note: This is implemented on top of MATSHELL to get support for scaling and shifting without requiring duplicate code
621: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),
622: MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
623: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
624: MatMFFDGetH(),
625: M*/
626: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
627: {
628: MatMFFD mfctx;
632: MatMFFDInitializePackage();
634: PetscHeaderCreate(mfctx,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),NULL,NULL);
636: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
637: mfctx->recomputeperiod = 1;
638: mfctx->count = 0;
639: mfctx->currenth = 0.0;
640: mfctx->historyh = NULL;
641: mfctx->ncurrenth = 0;
642: mfctx->maxcurrenth = 0;
643: ((PetscObject)mfctx)->type_name = NULL;
645: /*
646: Create the empty data structure to contain compute-h routines.
647: These will be filled in below from the command line options or
648: a later call with MatMFFDSetType() or if that is not called
649: then it will default in the first use of MatMult_MFFD()
650: */
651: mfctx->ops->compute = NULL;
652: mfctx->ops->destroy = NULL;
653: mfctx->ops->view = NULL;
654: mfctx->ops->setfromoptions = NULL;
655: mfctx->hctx = NULL;
657: mfctx->func = NULL;
658: mfctx->funcctx = NULL;
659: mfctx->w = NULL;
660: mfctx->mat = A;
662: MatSetType(A,MATSHELL);
663: MatShellSetContext(A,mfctx);
664: MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_MFFD);
665: MatShellSetOperation(A,MATOP_DESTROY,(void (*)(void))MatDestroy_MFFD);
666: MatShellSetOperation(A,MATOP_VIEW,(void (*)(void))MatView_MFFD);
667: MatShellSetOperation(A,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MFFD);
668: MatShellSetOperation(A,MATOP_SET_FROM_OPTIONS,(void (*)(void))MatSetFromOptions_MFFD);
670: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
671: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
672: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
673: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
674: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
675: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
676: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
677: PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
678: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetHHistory_C",MatMFFDSetHHistory_MFFD);
679: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetType_C",MatMFFDSetType_MFFD);
680: PetscObjectComposeFunction((PetscObject)A,"MatMFFDGetH_C",MatMFFDGetH_MFFD);
681: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
682: return(0);
683: }
685: /*@
686: MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
688: Collective on Vec
690: Input Parameters:
691: + comm - MPI communicator
692: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
693: This value should be the same as the local size used in creating the
694: y vector for the matrix-vector product y = Ax.
695: . n - This value should be the same as the local size used in creating the
696: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
697: calculated if N is given) For square matrices n is almost always m.
698: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
699: - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
702: Output Parameter:
703: . J - the matrix-free matrix
705: Options Database Keys: call MatSetFromOptions() to trigger these
706: + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
707: . -mat_mffd_err - square root of estimated relative error in function evaluation
708: . -mat_mffd_period - how often h is recomputed, defaults to 1, everytime
709: . -mat_mffd_check_positivity - possibly decrease h until U + h*a has only positive values
710: - -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
711: (requires real valued functions but that PETSc be configured for complex numbers)
714: Level: advanced
716: Notes:
717: The matrix-free matrix context merely contains the function pointers
718: and work space for performing finite difference approximations of
719: Jacobian-vector products, F'(u)*a,
721: The default code uses the following approach to compute h
723: .vb
724: F'(u)*a = [F(u+h*a) - F(u)]/h where
725: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
726: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
727: where
728: error_rel = square root of relative error in function evaluation
729: umin = minimum iterate parameter
730: .ve
732: You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
733: preconditioner matrix
735: The user can set the error_rel via MatMFFDSetFunctionError() and
736: umin via MatMFFDDSSetUmin(); see Users-Manual: ch_snes for details.
738: The user should call MatDestroy() when finished with the matrix-free
739: matrix context.
741: Options Database Keys:
742: + -mat_mffd_err <error_rel> - Sets error_rel
743: . -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
744: - -mat_mffd_check_positivity
746: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
747: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
748: MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
750: @*/
751: PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
752: {
756: MatCreate(comm,J);
757: MatSetSizes(*J,m,n,M,N);
758: MatSetType(*J,MATMFFD);
759: MatSetUp(*J);
760: return(0);
761: }
763: /*@
764: MatMFFDGetH - Gets the last value that was used as the differencing
765: parameter.
767: Not Collective
769: Input Parameters:
770: . mat - the matrix obtained with MatCreateSNESMF()
772: Output Parameter:
773: . h - the differencing step size
775: Level: advanced
777: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
778: @*/
779: PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
780: {
786: PetscUseMethod(mat,"MatMFFDGetH_C",(Mat,PetscScalar*),(mat,h));
787: return(0);
788: }
790: /*@C
791: MatMFFDSetFunction - Sets the function used in applying the matrix free.
793: Logically Collective on Mat
795: Input Parameters:
796: + mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
797: . func - the function to use
798: - funcctx - optional function context passed to function
800: Calling Sequence of func:
801: $ func (void *funcctx, Vec x, Vec f)
803: + funcctx - user provided context
804: . x - input vector
805: - f - computed output function
807: Level: advanced
809: Notes:
810: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
811: matrix inside your compute Jacobian routine
813: If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
815: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
816: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
817: @*/
818: PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
819: {
824: PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
825: return(0);
826: }
828: /*@C
829: MatMFFDSetFunctioni - Sets the function for a single component
831: Logically Collective on Mat
833: Input Parameters:
834: + mat - the matrix free matrix created via MatCreateSNESMF()
835: - funci - the function to use
837: Level: advanced
839: Notes:
840: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
841: matrix inside your compute Jacobian routine.
842: This function is necessary to compute the diagonal of the matrix.
843: funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.
845: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
847: @*/
848: PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
849: {
854: PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
855: return(0);
856: }
858: /*@C
859: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
861: Logically Collective on Mat
863: Input Parameters:
864: + mat - the matrix free matrix created via MatCreateSNESMF()
865: - func - the function to use
867: Level: advanced
869: Notes:
870: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
871: matrix inside your compute Jacobian routine.
872: This function is necessary to compute the diagonal of the matrix.
875: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
876: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
877: @*/
878: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
879: {
884: PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
885: return(0);
886: }
888: /*@
889: MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime
891: Logically Collective on Mat
893: Input Parameters:
894: + mat - the matrix free matrix created via MatCreateSNESMF()
895: - period - 1 for everytime, 2 for every second etc
897: Options Database Keys:
898: . -mat_mffd_period <period>
900: Level: advanced
903: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
904: MatMFFDSetHHistory(), MatMFFDResetHHistory()
905: @*/
906: PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
907: {
913: PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
914: return(0);
915: }
917: /*@
918: MatMFFDSetFunctionError - Sets the error_rel for the approximation of
919: matrix-vector products using finite differences.
921: Logically Collective on Mat
923: Input Parameters:
924: + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
925: - error_rel - relative error (should be set to the square root of
926: the relative error in the function evaluations)
928: Options Database Keys:
929: . -mat_mffd_err <error_rel> - Sets error_rel
931: Level: advanced
933: Notes:
934: The default matrix-free matrix-vector product routine computes
935: .vb
936: F'(u)*a = [F(u+h*a) - F(u)]/h where
937: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
938: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
939: .ve
941: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
942: MatMFFDSetHHistory(), MatMFFDResetHHistory()
943: @*/
944: PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
945: {
951: PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
952: return(0);
953: }
955: /*@
956: MatMFFDSetHHistory - Sets an array to collect a history of the
957: differencing values (h) computed for the matrix-free product.
959: Logically Collective on Mat
961: Input Parameters:
962: + J - the matrix-free matrix context
963: . histroy - space to hold the history
964: - nhistory - number of entries in history, if more entries are generated than
965: nhistory, then the later ones are discarded
967: Level: advanced
969: Notes:
970: Use MatMFFDResetHHistory() to reset the history counter and collect
971: a new batch of differencing parameters, h.
973: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
974: MatMFFDResetHHistory(), MatMFFDSetFunctionError()
976: @*/
977: PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
978: {
985: PetscUseMethod(J,"MatMFFDSetHHistory_C",(Mat,PetscScalar[],PetscInt),(J,history,nhistory));
986: return(0);
987: }
989: /*@
990: MatMFFDResetHHistory - Resets the counter to zero to begin
991: collecting a new set of differencing histories.
993: Logically Collective on Mat
995: Input Parameters:
996: . J - the matrix-free matrix context
998: Level: advanced
1000: Notes:
1001: Use MatMFFDSetHHistory() to create the original history counter.
1003: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1004: MatMFFDSetHHistory(), MatMFFDSetFunctionError()
1006: @*/
1007: PetscErrorCode MatMFFDResetHHistory(Mat J)
1008: {
1013: PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1014: return(0);
1015: }
1017: /*@
1018: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1019: Jacobian are computed
1021: Logically Collective on Mat
1023: Input Parameters:
1024: + J - the MatMFFD matrix
1025: . U - the vector
1026: - F - (optional) vector that contains F(u) if it has been already computed
1028: Notes:
1029: This is rarely used directly
1031: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
1032: point during the first MatMult() after each call to MatMFFDSetBase().
1034: Level: advanced
1036: @*/
1037: PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
1038: {
1045: PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1046: return(0);
1047: }
1049: /*@C
1050: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1051: it to satisfy some criteria
1053: Logically Collective on Mat
1055: Input Parameters:
1056: + J - the MatMFFD matrix
1057: . fun - the function that checks h
1058: - ctx - any context needed by the function
1060: Options Database Keys:
1061: . -mat_mffd_check_positivity
1063: Level: advanced
1065: Notes:
1066: For example, MatMFFDCheckPositivity() insures that all entries
1067: of U + h*a are non-negative
1069: The function you provide is called after the default h has been computed and allows you to
1070: modify it.
1072: .seealso: MatMFFDCheckPositivity()
1073: @*/
1074: PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1075: {
1080: PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1081: return(0);
1082: }
1084: /*@
1085: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1086: zero, decreases h until this is satisfied.
1088: Logically Collective on Vec
1090: Input Parameters:
1091: + U - base vector that is added to
1092: . a - vector that is added
1093: . h - scaling factor on a
1094: - dummy - context variable (unused)
1096: Options Database Keys:
1097: . -mat_mffd_check_positivity
1099: Level: advanced
1101: Notes:
1102: This is rarely used directly, rather it is passed as an argument to
1103: MatMFFDSetCheckh()
1105: .seealso: MatMFFDSetCheckh()
1106: @*/
1107: PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1108: {
1109: PetscReal val, minval;
1110: PetscScalar *u_vec, *a_vec;
1112: PetscInt i,n;
1113: MPI_Comm comm;
1119: PetscObjectGetComm((PetscObject)U,&comm);
1120: VecGetArray(U,&u_vec);
1121: VecGetArray(a,&a_vec);
1122: VecGetLocalSize(U,&n);
1123: minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1124: for (i=0; i<n; i++) {
1125: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1126: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1127: if (val < minval) minval = val;
1128: }
1129: }
1130: VecRestoreArray(U,&u_vec);
1131: VecRestoreArray(a,&a_vec);
1132: MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1133: if (val <= PetscAbsScalar(*h)) {
1134: PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1135: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1136: else *h = -0.99*val;
1137: }
1138: return(0);
1139: }