Actual source code: mffd.c
petsc-3.12.5 2020-03-29
2: #include <petsc/private/matimpl.h>
3: #include <../src/mat/impls/mffd/mffdimpl.h>
5: PetscFunctionList MatMFFDList = 0;
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 exclusions */
55: PetscOptionsGetString(NULL,NULL,"-info_exclude",logList,sizeof(logList),&opt);
56: if (opt) {
57: PetscStrInList("matmffd",logList,',',&pkg);
58: if (pkg) {PetscInfoDeactivateClass(MATMFFD_CLASSID);}
59: }
60: /* Process summary exclusions */
61: PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
62: if (opt) {
63: PetscStrInList("matmffd",logList,',',&pkg);
64: if (pkg) {PetscLogEventExcludeClass(MATMFFD_CLASSID);}
65: }
66: /* Register package finalizer */
67: PetscRegisterFinalize(MatMFFDFinalizePackage);
68: return(0);
69: }
71: static PetscErrorCode MatMFFDSetType_MFFD(Mat mat,MatMFFDType ftype)
72: {
73: PetscErrorCode ierr,(*r)(MatMFFD);
74: MatMFFD ctx;
75: PetscBool match;
80: MatShellGetContext(mat,&ctx);
82: /* already set, so just return */
83: PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
84: if (match) return(0);
86: /* destroy the old one if it exists */
87: if (ctx->ops->destroy) {
88: (*ctx->ops->destroy)(ctx);
89: }
91: PetscFunctionListFind(MatMFFDList,ftype,&r);
92: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
93: (*r)(ctx);
94: PetscObjectChangeTypeName((PetscObject)ctx,ftype);
95: return(0);
96: }
98: /*@C
99: MatMFFDSetType - Sets the method that is used to compute the
100: differencing parameter for finite differene matrix-free formulations.
102: Input Parameters:
103: + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
104: or MatSetType(mat,MATMFFD);
105: - ftype - the type requested, either MATMFFD_WP or MATMFFD_DS
107: Level: advanced
109: Notes:
110: For example, such routines can compute h for use in
111: Jacobian-vector products of the form
113: F(x+ha) - F(x)
114: F'(u)a ~= ----------------
115: h
117: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction(), MatCreateMFFD()
118: @*/
119: PetscErrorCode MatMFFDSetType(Mat mat,MatMFFDType ftype)
120: {
126: PetscTryMethod(mat,"MatMFFDSetType_C",(Mat,MatMFFDType),(mat,ftype));
127: return(0);
128: }
130: static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);
132: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
133: static PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
134: {
135: MatMFFD ctx;
139: MatShellGetContext(mat,&ctx);
140: ctx->funcisetbase = func;
141: return(0);
142: }
144: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
145: static PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
146: {
147: MatMFFD ctx;
151: MatShellGetContext(mat,&ctx);
152: ctx->funci = funci;
153: MatShellSetOperation(mat,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_MFFD);
154: return(0);
155: }
157: static PetscErrorCode MatMFFDGetH_MFFD(Mat mat,PetscScalar *h)
158: {
159: MatMFFD ctx;
163: MatShellGetContext(mat,&ctx);
164: *h = ctx->currenth;
165: return(0);
166: }
168: static PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
169: {
170: MatMFFD ctx;
174: MatShellGetContext(J,&ctx);
175: ctx->ncurrenth = 0;
176: return(0);
177: }
179: /*@C
180: MatMFFDRegister - Adds a method to the MatMFFD registry.
182: Not Collective
184: Input Parameters:
185: + name_solver - name of a new user-defined compute-h module
186: - routine_create - routine to create method context
188: Level: developer
190: Notes:
191: MatMFFDRegister() may be called multiple times to add several user-defined solvers.
193: Sample usage:
194: .vb
195: MatMFFDRegister("my_h",MyHCreate);
196: .ve
198: Then, your solver can be chosen with the procedural interface via
199: $ MatMFFDSetType(mfctx,"my_h")
200: or at runtime via the option
201: $ -mat_mffd_type my_h
203: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
204: @*/
205: PetscErrorCode MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
206: {
210: MatInitializePackage();
211: PetscFunctionListAdd(&MatMFFDList,sname,function);
212: return(0);
213: }
215: /* ----------------------------------------------------------------------------------------*/
216: static PetscErrorCode MatDestroy_MFFD(Mat mat)
217: {
219: MatMFFD ctx;
222: MatShellGetContext(mat,&ctx);
223: VecDestroy(&ctx->w);
224: VecDestroy(&ctx->current_u);
225: if (ctx->current_f_allocated) {
226: VecDestroy(&ctx->current_f);
227: }
228: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
229: PetscHeaderDestroy(&ctx);
231: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
232: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
233: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
234: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
235: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
236: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
237: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
238: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
239: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetHHistory_C",NULL);
240: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetType_C",NULL);
241: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDGetH_C",NULL);
242: return(0);
243: }
245: /*
246: MatMFFDView_MFFD - Views matrix-free parameters.
248: */
249: static PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
250: {
252: MatMFFD ctx;
253: PetscBool iascii, viewbase, viewfunction;
254: const char *prefix;
257: MatShellGetContext(J,&ctx);
258: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
259: if (iascii) {
260: PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
261: PetscViewerASCIIPushTab(viewer);
262: PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
263: if (!((PetscObject)ctx)->type_name) {
264: PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
265: } else {
266: PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
267: }
268: #if defined(PETSC_USE_COMPLEX)
269: if (ctx->usecomplex) {
270: PetscViewerASCIIPrintf(viewer,"Using Lyness complex number trick to compute the matrix-vector product\n");
271: }
272: #endif
273: if (ctx->ops->view) {
274: (*ctx->ops->view)(ctx,viewer);
275: }
276: PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);
278: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_base", &viewbase);
279: if (viewbase) {
280: PetscViewerASCIIPrintf(viewer, "Base:\n");
281: VecView(ctx->current_u, viewer);
282: }
283: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_function", &viewfunction);
284: if (viewfunction) {
285: PetscViewerASCIIPrintf(viewer, "Function:\n");
286: VecView(ctx->current_f, viewer);
287: }
288: PetscViewerASCIIPopTab(viewer);
289: }
290: return(0);
291: }
293: /*
294: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
295: allows the user to indicate the beginning of a new linear solve by calling
296: MatAssemblyXXX() on the matrix free matrix. This then allows the
297: MatCreateMFFD_WP() to properly compute ||U|| only the first time
298: in the linear solver rather than every time.
300: This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
301: it must be labeled as PETSC_EXTERN
302: */
303: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
304: {
306: MatMFFD j;
309: MatShellGetContext(J,&j);
310: MatMFFDResetHHistory(J);
311: return(0);
312: }
314: /*
315: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
317: y ~= (F(u + ha) - F(u))/h,
318: where F = nonlinear function, as set by SNESSetFunction()
319: u = current iterate
320: h = difference interval
321: */
322: static PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
323: {
324: MatMFFD ctx;
325: PetscScalar h;
326: Vec w,U,F;
328: PetscBool zeroa;
331: MatShellGetContext(mat,&ctx);
332: 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");
333: /* We log matrix-free matrix-vector products separately, so that we can
334: separate the performance monitoring from the cases that use conventional
335: storage. We may eventually modify event logging to associate events
336: with particular objects, hence alleviating the more general problem. */
337: PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);
339: w = ctx->w;
340: U = ctx->current_u;
341: F = ctx->current_f;
342: /*
343: Compute differencing parameter
344: */
345: if (!((PetscObject)ctx)->type_name) {
346: MatMFFDSetType(mat,MATMFFD_WP);
347: MatSetFromOptions(mat);
348: }
349: (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
350: if (zeroa) {
351: VecSet(y,0.0);
352: return(0);
353: }
355: if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
356: if (ctx->checkh) {
357: (*ctx->checkh)(ctx->checkhctx,U,a,&h);
358: }
360: /* keep a record of the current differencing parameter h */
361: ctx->currenth = h;
362: #if defined(PETSC_USE_COMPLEX)
363: PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
364: #else
365: PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
366: #endif
367: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
368: ctx->historyh[ctx->ncurrenth] = h;
369: }
370: ctx->ncurrenth++;
372: #if defined(PETSC_USE_COMPLEX)
373: if (ctx->usecomplex) h = PETSC_i*h;
374: #endif
376: /* w = u + ha */
377: VecWAXPY(w,h,a,U);
379: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
380: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
381: (*ctx->func)(ctx->funcctx,U,F);
382: }
383: (*ctx->func)(ctx->funcctx,w,y);
385: #if defined(PETSC_USE_COMPLEX)
386: if (ctx->usecomplex) {
387: VecImaginaryPart(y);
388: h = PetscImaginaryPart(h);
389: } else {
390: VecAXPY(y,-1.0,F);
391: }
392: #else
393: VecAXPY(y,-1.0,F);
394: #endif
395: VecScale(y,1.0/h);
396: if (mat->nullsp) {MatNullSpaceRemove(mat->nullsp,y);}
398: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
399: return(0);
400: }
402: /*
403: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
405: y ~= (F(u + ha) - F(u))/h,
406: where F = nonlinear function, as set by SNESSetFunction()
407: u = current iterate
408: h = difference interval
409: */
410: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
411: {
412: MatMFFD ctx;
413: PetscScalar h,*aa,*ww,v;
414: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
415: Vec w,U;
417: PetscInt i,rstart,rend;
420: MatShellGetContext(mat,&ctx);
421: if (!ctx->func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunction() first");
422: if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioni() first");
423: w = ctx->w;
424: U = ctx->current_u;
425: (*ctx->func)(ctx->funcctx,U,a);
426: if (ctx->funcisetbase) {
427: (*ctx->funcisetbase)(ctx->funcctx,U);
428: }
429: VecCopy(U,w);
431: VecGetOwnershipRange(a,&rstart,&rend);
432: VecGetArray(a,&aa);
433: for (i=rstart; i<rend; i++) {
434: VecGetArray(w,&ww);
435: h = ww[i-rstart];
436: if (h == 0.0) h = 1.0;
437: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
438: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
439: h *= epsilon;
441: ww[i-rstart] += h;
442: VecRestoreArray(w,&ww);
443: (*ctx->funci)(ctx->funcctx,i,w,&v);
444: aa[i-rstart] = (v - aa[i-rstart])/h;
446: VecGetArray(w,&ww);
447: ww[i-rstart] -= h;
448: VecRestoreArray(w,&ww);
449: }
450: VecRestoreArray(a,&aa);
451: return(0);
452: }
454: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
455: {
457: MatMFFD ctx;
460: MatShellGetContext(J,&ctx);
461: MatMFFDResetHHistory(J);
462: if (!ctx->current_u) {
463: VecDuplicate(U,&ctx->current_u);
464: VecLockReadPush(ctx->current_u);
465: }
466: VecLockReadPop(ctx->current_u);
467: VecCopy(U,ctx->current_u);
468: VecLockReadPush(ctx->current_u);
469: if (F) {
470: if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
471: ctx->current_f = F;
472: ctx->current_f_allocated = PETSC_FALSE;
473: } else if (!ctx->current_f_allocated) {
474: MatCreateVecs(J,NULL,&ctx->current_f);
476: ctx->current_f_allocated = PETSC_TRUE;
477: }
478: if (!ctx->w) {
479: VecDuplicate(ctx->current_u,&ctx->w);
480: }
481: J->assembled = PETSC_TRUE;
482: return(0);
483: }
485: typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
487: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
488: {
489: MatMFFD ctx;
493: MatShellGetContext(J,&ctx);
494: ctx->checkh = fun;
495: ctx->checkhctx = ectx;
496: return(0);
497: }
499: /*@C
500: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
501: MatMFFD options in the database.
503: Collective on Mat
505: Input Parameter:
506: + A - the Mat context
507: - prefix - the prefix to prepend to all option names
509: Notes:
510: A hyphen (-) must NOT be given at the beginning of the prefix name.
511: The first character of all runtime options is AUTOMATICALLY the hyphen.
513: Level: advanced
515: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
516: @*/
517: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
519: {
520: MatMFFD mfctx;
525: MatShellGetContext(mat,&mfctx);
527: PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
528: return(0);
529: }
531: static PetscErrorCode MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
532: {
533: MatMFFD mfctx;
535: PetscBool flg;
536: char ftype[256];
540: MatShellGetContext(mat,&mfctx);
542: PetscObjectOptionsBegin((PetscObject)mfctx);
543: PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
544: if (flg) {
545: MatMFFDSetType(mat,ftype);
546: }
548: PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
549: PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
551: flg = PETSC_FALSE;
552: PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
553: if (flg) {
554: MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
555: }
556: #if defined(PETSC_USE_COMPLEX)
557: PetscOptionsBool("-mat_mffd_complex","Use Lyness complex number trick to compute the matrix-vector product","None",mfctx->usecomplex,&mfctx->usecomplex,NULL);
558: #endif
559: if (mfctx->ops->setfromoptions) {
560: (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
561: }
562: PetscOptionsEnd();
563: return(0);
564: }
566: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
567: {
568: MatMFFD ctx;
572: MatShellGetContext(mat,&ctx);
573: ctx->recomputeperiod = period;
574: return(0);
575: }
577: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
578: {
579: MatMFFD ctx;
583: MatShellGetContext(mat,&ctx);
584: ctx->func = func;
585: ctx->funcctx = funcctx;
586: return(0);
587: }
589: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
590: {
591: MatMFFD ctx;
595: MatShellGetContext(mat,&ctx);
596: if (error != PETSC_DEFAULT) ctx->error_rel = error;
597: return(0);
598: }
600: PetscErrorCode MatMFFDSetHHistory_MFFD(Mat J,PetscScalar history[],PetscInt nhistory)
601: {
602: MatMFFD ctx;
606: MatShellGetContext(J,&ctx);
607: ctx->historyh = history;
608: ctx->maxcurrenth = nhistory;
609: ctx->currenth = 0.;
610: return(0);
611: }
613: /*MC
614: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
616: Level: advanced
618: Developers Note: This is implemented on top of MATSHELL to get support for scaling and shifting without requiring duplicate code
620: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),
621: MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
622: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
623: MatMFFDGetH(),
624: M*/
625: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
626: {
627: MatMFFD mfctx;
631: MatMFFDInitializePackage();
633: PetscHeaderCreate(mfctx,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),NULL,NULL);
635: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
636: mfctx->recomputeperiod = 1;
637: mfctx->count = 0;
638: mfctx->currenth = 0.0;
639: mfctx->historyh = NULL;
640: mfctx->ncurrenth = 0;
641: mfctx->maxcurrenth = 0;
642: ((PetscObject)mfctx)->type_name = 0;
644: /*
645: Create the empty data structure to contain compute-h routines.
646: These will be filled in below from the command line options or
647: a later call with MatMFFDSetType() or if that is not called
648: then it will default in the first use of MatMult_MFFD()
649: */
650: mfctx->ops->compute = 0;
651: mfctx->ops->destroy = 0;
652: mfctx->ops->view = 0;
653: mfctx->ops->setfromoptions = 0;
654: mfctx->hctx = 0;
656: mfctx->func = 0;
657: mfctx->funcctx = 0;
658: mfctx->w = NULL;
659: mfctx->mat = A;
661: MatSetType(A,MATSHELL);
662: MatShellSetContext(A,mfctx);
663: MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_MFFD);
664: MatShellSetOperation(A,MATOP_DESTROY,(void (*)(void))MatDestroy_MFFD);
665: MatShellSetOperation(A,MATOP_VIEW,(void (*)(void))MatView_MFFD);
666: MatShellSetOperation(A,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MFFD);
667: MatShellSetOperation(A,MATOP_SET_FROM_OPTIONS,(void (*)(void))MatSetFromOptions_MFFD);
669: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
670: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
671: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
672: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
673: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
674: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
675: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
676: PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
677: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetHHistory_C",MatMFFDSetHHistory_MFFD);
678: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetType_C",MatMFFDSetType_MFFD);
679: PetscObjectComposeFunction((PetscObject)A,"MatMFFDGetH_C",MatMFFDGetH_MFFD);
680: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
681: return(0);
682: }
684: /*@
685: MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
687: Collective on Vec
689: Input Parameters:
690: + comm - MPI communicator
691: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
692: This value should be the same as the local size used in creating the
693: y vector for the matrix-vector product y = Ax.
694: . n - This value should be the same as the local size used in creating the
695: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
696: calculated if N is given) For square matrices n is almost always m.
697: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
698: - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
701: Output Parameter:
702: . J - the matrix-free matrix
704: Options Database Keys: call MatSetFromOptions() to trigger these
705: + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
706: . -mat_mffd_err - square root of estimated relative error in function evaluation
707: . -mat_mffd_period - how often h is recomputed, defaults to 1, everytime
708: . -mat_mffd_check_positivity - possibly decrease h until U + h*a has only positive values
709: - -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
710: (requires real valued functions but that PETSc be configured for complex numbers)
713: Level: advanced
715: Notes:
716: The matrix-free matrix context merely contains the function pointers
717: and work space for performing finite difference approximations of
718: Jacobian-vector products, F'(u)*a,
720: The default code uses the following approach to compute h
722: .vb
723: F'(u)*a = [F(u+h*a) - F(u)]/h where
724: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
725: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
726: where
727: error_rel = square root of relative error in function evaluation
728: umin = minimum iterate parameter
729: .ve
731: You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
732: preconditioner matrix
734: The user can set the error_rel via MatMFFDSetFunctionError() and
735: umin via MatMFFDDSSetUmin(); see Users-Manual: ch_snes for details.
737: The user should call MatDestroy() when finished with the matrix-free
738: matrix context.
740: Options Database Keys:
741: + -mat_mffd_err <error_rel> - Sets error_rel
742: . -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
743: - -mat_mffd_check_positivity
745: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
746: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
747: MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
749: @*/
750: PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
751: {
755: MatCreate(comm,J);
756: MatSetSizes(*J,m,n,M,N);
757: MatSetType(*J,MATMFFD);
758: MatSetUp(*J);
759: return(0);
760: }
762: /*@
763: MatMFFDGetH - Gets the last value that was used as the differencing
764: parameter.
766: Not Collective
768: Input Parameters:
769: . mat - the matrix obtained with MatCreateSNESMF()
771: Output Paramter:
772: . h - the differencing step size
774: Level: advanced
776: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
777: @*/
778: PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
779: {
785: PetscUseMethod(mat,"MatMFFDGetH_C",(Mat,PetscScalar*),(mat,h));
786: return(0);
787: }
789: /*@C
790: MatMFFDSetFunction - Sets the function used in applying the matrix free.
792: Logically Collective on Mat
794: Input Parameters:
795: + mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
796: . func - the function to use
797: - funcctx - optional function context passed to function
799: Calling Sequence of func:
800: $ func (void *funcctx, Vec x, Vec f)
802: + funcctx - user provided context
803: . x - input vector
804: - f - computed output function
806: Level: advanced
808: Notes:
809: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
810: matrix inside your compute Jacobian routine
812: If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
814: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
815: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
816: @*/
817: PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
818: {
823: PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
824: return(0);
825: }
827: /*@C
828: MatMFFDSetFunctioni - Sets the function for a single component
830: Logically Collective on Mat
832: Input Parameters:
833: + mat - the matrix free matrix created via MatCreateSNESMF()
834: - funci - the function to use
836: Level: advanced
838: Notes:
839: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
840: matrix inside your compute Jacobian routine.
841: This function is necessary to compute the diagonal of the matrix.
842: funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.
844: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
846: @*/
847: PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
848: {
853: PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
854: return(0);
855: }
857: /*@C
858: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
860: Logically Collective on Mat
862: Input Parameters:
863: + mat - the matrix free matrix created via MatCreateSNESMF()
864: - func - the function to use
866: Level: advanced
868: Notes:
869: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
870: matrix inside your compute Jacobian routine.
871: This function is necessary to compute the diagonal of the matrix.
874: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
875: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
876: @*/
877: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
878: {
883: PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
884: return(0);
885: }
887: /*@
888: MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime
890: Logically Collective on Mat
892: Input Parameters:
893: + mat - the matrix free matrix created via MatCreateSNESMF()
894: - period - 1 for everytime, 2 for every second etc
896: Options Database Keys:
897: . -mat_mffd_period <period>
899: Level: advanced
902: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
903: MatMFFDSetHHistory(), MatMFFDResetHHistory()
904: @*/
905: PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
906: {
912: PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
913: return(0);
914: }
916: /*@
917: MatMFFDSetFunctionError - Sets the error_rel for the approximation of
918: matrix-vector products using finite differences.
920: Logically Collective on Mat
922: Input Parameters:
923: + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
924: - error_rel - relative error (should be set to the square root of
925: the relative error in the function evaluations)
927: Options Database Keys:
928: . -mat_mffd_err <error_rel> - Sets error_rel
930: Level: advanced
932: Notes:
933: The default matrix-free matrix-vector product routine computes
934: .vb
935: F'(u)*a = [F(u+h*a) - F(u)]/h where
936: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
937: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
938: .ve
940: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
941: MatMFFDSetHHistory(), MatMFFDResetHHistory()
942: @*/
943: PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
944: {
950: PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
951: return(0);
952: }
954: /*@
955: MatMFFDSetHHistory - Sets an array to collect a history of the
956: differencing values (h) computed for the matrix-free product.
958: Logically Collective on Mat
960: Input Parameters:
961: + J - the matrix-free matrix context
962: . histroy - space to hold the history
963: - nhistory - number of entries in history, if more entries are generated than
964: nhistory, then the later ones are discarded
966: Level: advanced
968: Notes:
969: Use MatMFFDResetHHistory() to reset the history counter and collect
970: a new batch of differencing parameters, h.
972: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
973: MatMFFDResetHHistory(), MatMFFDSetFunctionError()
975: @*/
976: PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
977: {
984: PetscUseMethod(J,"MatMFFDSetHHistory_C",(Mat,PetscScalar[],PetscInt),(J,history,nhistory));
985: return(0);
986: }
988: /*@
989: MatMFFDResetHHistory - Resets the counter to zero to begin
990: collecting a new set of differencing histories.
992: Logically Collective on Mat
994: Input Parameters:
995: . J - the matrix-free matrix context
997: Level: advanced
999: Notes:
1000: Use MatMFFDSetHHistory() to create the original history counter.
1002: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1003: MatMFFDSetHHistory(), MatMFFDSetFunctionError()
1005: @*/
1006: PetscErrorCode MatMFFDResetHHistory(Mat J)
1007: {
1012: PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1013: return(0);
1014: }
1016: /*@
1017: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1018: Jacobian are computed
1020: Logically Collective on Mat
1022: Input Parameters:
1023: + J - the MatMFFD matrix
1024: . U - the vector
1025: - F - (optional) vector that contains F(u) if it has been already computed
1027: Notes:
1028: This is rarely used directly
1030: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
1031: point during the first MatMult() after each call to MatMFFDSetBase().
1033: Level: advanced
1035: @*/
1036: PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
1037: {
1044: PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1045: return(0);
1046: }
1048: /*@C
1049: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1050: it to satisfy some criteria
1052: Logically Collective on Mat
1054: Input Parameters:
1055: + J - the MatMFFD matrix
1056: . fun - the function that checks h
1057: - ctx - any context needed by the function
1059: Options Database Keys:
1060: . -mat_mffd_check_positivity
1062: Level: advanced
1064: Notes:
1065: For example, MatMFFDCheckPositivity() insures that all entries
1066: of U + h*a are non-negative
1068: The function you provide is called after the default h has been computed and allows you to
1069: modify it.
1071: .seealso: MatMFFDCheckPositivity()
1072: @*/
1073: PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1074: {
1079: PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1080: return(0);
1081: }
1083: /*@
1084: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1085: zero, decreases h until this is satisfied.
1087: Logically Collective on Vec
1089: Input Parameters:
1090: + U - base vector that is added to
1091: . a - vector that is added
1092: . h - scaling factor on a
1093: - dummy - context variable (unused)
1095: Options Database Keys:
1096: . -mat_mffd_check_positivity
1098: Level: advanced
1100: Notes:
1101: This is rarely used directly, rather it is passed as an argument to
1102: MatMFFDSetCheckh()
1104: .seealso: MatMFFDSetCheckh()
1105: @*/
1106: PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1107: {
1108: PetscReal val, minval;
1109: PetscScalar *u_vec, *a_vec;
1111: PetscInt i,n;
1112: MPI_Comm comm;
1118: PetscObjectGetComm((PetscObject)U,&comm);
1119: VecGetArray(U,&u_vec);
1120: VecGetArray(a,&a_vec);
1121: VecGetLocalSize(U,&n);
1122: minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1123: for (i=0; i<n; i++) {
1124: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1125: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1126: if (val < minval) minval = val;
1127: }
1128: }
1129: VecRestoreArray(U,&u_vec);
1130: VecRestoreArray(a,&a_vec);
1131: MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1132: if (val <= PetscAbsScalar(*h)) {
1133: PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1134: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1135: else *h = -0.99*val;
1136: }
1137: return(0);
1138: }