Actual source code: mffd.c
petsc-3.10.5 2019-03-28
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: .keywords: Petsc, destroy, package
19: .seealso: PetscFinalize(), MatCreateMFFD(), MatCreateSNESMF()
20: @*/
21: PetscErrorCode MatMFFDFinalizePackage(void)
22: {
26: PetscFunctionListDestroy(&MatMFFDList);
27: MatMFFDPackageInitialized = PETSC_FALSE;
28: MatMFFDRegisterAllCalled = PETSC_FALSE;
29: return(0);
30: }
32: /*@C
33: MatMFFDInitializePackage - This function initializes everything in the MatMFFD package. It is called
34: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to MatCreate_MFFD()
35: when using static libraries.
37: Level: developer
39: .keywords: Vec, initialize, package
40: .seealso: PetscInitialize()
41: @*/
42: PetscErrorCode MatMFFDInitializePackage(void)
43: {
44: char logList[256];
45: PetscBool opt,pkg;
49: if (MatMFFDPackageInitialized) return(0);
50: MatMFFDPackageInitialized = PETSC_TRUE;
51: /* Register Classes */
52: PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);
53: /* Register Constructors */
54: MatMFFDRegisterAll();
55: /* Register Events */
56: PetscLogEventRegister("MatMult MF",MATMFFD_CLASSID,&MATMFFD_Mult);
57: /* Process info exclusions */
58: PetscOptionsGetString(NULL,NULL,"-info_exclude",logList,sizeof(logList),&opt);
59: if (opt) {
60: PetscStrInList("matmffd",logList,',',&pkg);
61: if (pkg) {PetscInfoDeactivateClass(MATMFFD_CLASSID);}
62: }
63: /* Process summary exclusions */
64: PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
65: if (opt) {
66: PetscStrInList("matmffd",logList,',',&pkg);
67: if (pkg) {PetscLogEventExcludeClass(MATMFFD_CLASSID);}
68: }
69: /* Register package finalizer */
70: PetscRegisterFinalize(MatMFFDFinalizePackage);
71: return(0);
72: }
74: /*@C
75: MatMFFDSetType - Sets the method that is used to compute the
76: differencing parameter for finite differene matrix-free formulations.
78: Input Parameters:
79: + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
80: or MatSetType(mat,MATMFFD);
81: - ftype - the type requested, either MATMFFD_WP or MATMFFD_DS
83: Level: advanced
85: Notes:
86: For example, such routines can compute h for use in
87: Jacobian-vector products of the form
89: F(x+ha) - F(x)
90: F'(u)a ~= ----------------
91: h
93: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction(), MatCreateMFFD()
94: @*/
95: PetscErrorCode MatMFFDSetType(Mat mat,MatMFFDType ftype)
96: {
97: PetscErrorCode ierr,(*r)(MatMFFD);
98: MatMFFD ctx = (MatMFFD)mat->data;
99: PetscBool match;
105: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
106: if (!match) return(0);
108: /* already set, so just return */
109: PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
110: if (match) return(0);
112: /* destroy the old one if it exists */
113: if (ctx->ops->destroy) {
114: (*ctx->ops->destroy)(ctx);
115: }
117: PetscFunctionListFind(MatMFFDList,ftype,&r);
118: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
119: (*r)(ctx);
120: PetscObjectChangeTypeName((PetscObject)ctx,ftype);
121: return(0);
122: }
124: static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);
126: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
127: static PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
128: {
129: MatMFFD ctx = (MatMFFD)mat->data;
132: ctx->funcisetbase = func;
133: /* allow users to compose their own getdiagonal and allow MatHasOperation
134: to return false if the two functions pointers are not set */
135: if (!mat->ops->getdiagonal && func) {
136: mat->ops->getdiagonal = MatGetDiagonal_MFFD;
137: }
138: return(0);
139: }
141: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
142: static PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
143: {
144: MatMFFD ctx = (MatMFFD)mat->data;
147: ctx->funci = funci;
148: /* allow users to compose their own getdiagonal and allow MatHasOperation
149: to return false if the two functions pointers are not set */
150: if (!mat->ops->getdiagonal && funci) {
151: mat->ops->getdiagonal = MatGetDiagonal_MFFD;
152: }
153: return(0);
154: }
156: static PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
157: {
158: MatMFFD ctx = (MatMFFD)J->data;
161: ctx->ncurrenth = 0;
162: return(0);
163: }
165: /*@C
166: MatMFFDRegister - Adds a method to the MatMFFD registry.
168: Not Collective
170: Input Parameters:
171: + name_solver - name of a new user-defined compute-h module
172: - routine_create - routine to create method context
174: Level: developer
176: Notes:
177: MatMFFDRegister() may be called multiple times to add several user-defined solvers.
179: Sample usage:
180: .vb
181: MatMFFDRegister("my_h",MyHCreate);
182: .ve
184: Then, your solver can be chosen with the procedural interface via
185: $ MatMFFDSetType(mfctx,"my_h")
186: or at runtime via the option
187: $ -mat_mffd_type my_h
189: .keywords: MatMFFD, register
191: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
192: @*/
193: PetscErrorCode MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
194: {
198: MatInitializePackage();
199: PetscFunctionListAdd(&MatMFFDList,sname,function);
200: return(0);
201: }
203: /* ----------------------------------------------------------------------------------------*/
204: static PetscErrorCode MatDestroy_MFFD(Mat mat)
205: {
207: MatMFFD ctx = (MatMFFD)mat->data;
210: VecDestroy(&ctx->w);
211: VecDestroy(&ctx->drscale);
212: VecDestroy(&ctx->dlscale);
213: VecDestroy(&ctx->dshift);
214: VecDestroy(&ctx->dshiftw);
215: VecDestroy(&ctx->current_u);
216: if (ctx->current_f_allocated) {
217: VecDestroy(&ctx->current_f);
218: }
219: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
220: PetscHeaderDestroy(&ctx);
221: mat->data = 0;
223: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
224: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
225: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
226: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
227: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
228: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
229: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
230: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
231: return(0);
232: }
234: /*
235: MatMFFDView_MFFD - Views matrix-free parameters.
237: */
238: static PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
239: {
241: MatMFFD ctx = (MatMFFD)J->data;
242: PetscBool iascii, viewbase, viewfunction;
243: const char *prefix;
246: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
247: if (iascii) {
248: PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
249: PetscViewerASCIIPushTab(viewer);
250: PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
251: if (!((PetscObject)ctx)->type_name) {
252: PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
253: } else {
254: PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
255: }
256: if (ctx->ops->view) {
257: (*ctx->ops->view)(ctx,viewer);
258: }
259: PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);
261: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_base", &viewbase);
262: if (viewbase) {
263: PetscViewerASCIIPrintf(viewer, "Base:\n");
264: VecView(ctx->current_u, viewer);
265: }
266: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_function", &viewfunction);
267: if (viewfunction) {
268: PetscViewerASCIIPrintf(viewer, "Function:\n");
269: VecView(ctx->current_f, viewer);
270: }
271: PetscViewerASCIIPopTab(viewer);
272: }
273: return(0);
274: }
276: /*
277: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
278: allows the user to indicate the beginning of a new linear solve by calling
279: MatAssemblyXXX() on the matrix free matrix. This then allows the
280: MatCreateMFFD_WP() to properly compute ||U|| only the first time
281: in the linear solver rather than every time.
283: This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
284: it must be labeled as PETSC_EXTERN
285: */
286: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
287: {
289: MatMFFD j = (MatMFFD)J->data;
292: MatMFFDResetHHistory(J);
293: j->vshift = 0.0;
294: j->vscale = 1.0;
295: return(0);
296: }
298: /*
299: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
301: y ~= (F(u + ha) - F(u))/h,
302: where F = nonlinear function, as set by SNESSetFunction()
303: u = current iterate
304: h = difference interval
305: */
306: static PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
307: {
308: MatMFFD ctx = (MatMFFD)mat->data;
309: PetscScalar h;
310: Vec w,U,F;
312: PetscBool zeroa;
315: 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");
316: /* We log matrix-free matrix-vector products separately, so that we can
317: separate the performance monitoring from the cases that use conventional
318: storage. We may eventually modify event logging to associate events
319: with particular objects, hence alleviating the more general problem. */
320: PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);
322: w = ctx->w;
323: U = ctx->current_u;
324: F = ctx->current_f;
325: /*
326: Compute differencing parameter
327: */
328: if (!((PetscObject)ctx)->type_name) {
329: MatMFFDSetType(mat,MATMFFD_WP);
330: MatSetFromOptions(mat);
331: }
332: (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
333: if (zeroa) {
334: VecSet(y,0.0);
335: return(0);
336: }
338: if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
339: if (ctx->checkh) {
340: (*ctx->checkh)(ctx->checkhctx,U,a,&h);
341: }
343: /* keep a record of the current differencing parameter h */
344: ctx->currenth = h;
345: #if defined(PETSC_USE_COMPLEX)
346: PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
347: #else
348: PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
349: #endif
350: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
351: ctx->historyh[ctx->ncurrenth] = h;
352: }
353: ctx->ncurrenth++;
355: /* w = u + ha */
356: if (ctx->drscale) {
357: VecPointwiseMult(ctx->drscale,a,U);
358: VecAYPX(U,h,w);
359: } else {
360: VecWAXPY(w,h,a,U);
361: }
363: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
364: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
365: (*ctx->func)(ctx->funcctx,U,F);
366: }
367: (*ctx->func)(ctx->funcctx,w,y);
369: VecAXPY(y,-1.0,F);
370: VecScale(y,1.0/h);
372: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
373: VecAXPBY(y,ctx->vshift,ctx->vscale,a);
374: }
375: if (ctx->dlscale) {
376: VecPointwiseMult(y,ctx->dlscale,y);
377: }
378: if (ctx->dshift) {
379: if (!ctx->dshiftw) {
380: VecDuplicate(y,&ctx->dshiftw);
381: }
382: VecPointwiseMult(ctx->dshift,a,ctx->dshiftw);
383: VecAXPY(y,1.0,ctx->dshiftw);
384: }
386: if (mat->nullsp) {MatNullSpaceRemove(mat->nullsp,y);}
388: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
389: return(0);
390: }
392: /*
393: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
395: y ~= (F(u + ha) - F(u))/h,
396: where F = nonlinear function, as set by SNESSetFunction()
397: u = current iterate
398: h = difference interval
399: */
400: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
401: {
402: MatMFFD ctx = (MatMFFD)mat->data;
403: PetscScalar h,*aa,*ww,v;
404: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
405: Vec w,U;
407: PetscInt i,rstart,rend;
410: if (!ctx->func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunction() first");
411: if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioni() first");
412: if (!ctx->funcisetbase) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioniBase() first");
413: w = ctx->w;
414: U = ctx->current_u;
415: (*ctx->func)(ctx->funcctx,U,a);
416: (*ctx->funcisetbase)(ctx->funcctx,U);
417: VecCopy(U,w);
419: VecGetOwnershipRange(a,&rstart,&rend);
420: VecGetArray(a,&aa);
421: for (i=rstart; i<rend; i++) {
422: VecGetArray(w,&ww);
423: h = ww[i-rstart];
424: if (h == 0.0) h = 1.0;
425: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
426: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
427: h *= epsilon;
429: ww[i-rstart] += h;
430: VecRestoreArray(w,&ww);
431: (*ctx->funci)(ctx->funcctx,i,w,&v);
432: aa[i-rstart] = (v - aa[i-rstart])/h;
434: /* possibly shift and scale result */
435: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
436: aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
437: }
439: VecGetArray(w,&ww);
440: ww[i-rstart] -= h;
441: VecRestoreArray(w,&ww);
442: }
443: VecRestoreArray(a,&aa);
444: return(0);
445: }
447: static PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
448: {
449: MatMFFD aij = (MatMFFD)mat->data;
453: if (ll && !aij->dlscale) {
454: VecDuplicate(ll,&aij->dlscale);
455: }
456: if (rr && !aij->drscale) {
457: VecDuplicate(rr,&aij->drscale);
458: }
459: if (ll) {
460: VecCopy(ll,aij->dlscale);
461: }
462: if (rr) {
463: VecCopy(rr,aij->drscale);
464: }
465: return(0);
466: }
468: static PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
469: {
470: MatMFFD aij = (MatMFFD)mat->data;
474: if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
475: if (!aij->dshift) {
476: VecDuplicate(ll,&aij->dshift);
477: }
478: VecAXPY(aij->dshift,1.0,ll);
479: return(0);
480: }
482: static PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
483: {
484: MatMFFD shell = (MatMFFD)Y->data;
487: shell->vshift += a;
488: return(0);
489: }
491: static PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
492: {
493: MatMFFD shell = (MatMFFD)Y->data;
496: shell->vscale *= a;
497: return(0);
498: }
500: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
501: {
503: MatMFFD ctx = (MatMFFD)J->data;
506: MatMFFDResetHHistory(J);
507: if (!ctx->current_u) {
508: VecDuplicate(U,&ctx->current_u);
509: VecLockPush(ctx->current_u);
510: }
511: VecLockPop(ctx->current_u);
512: VecCopy(U,ctx->current_u);
513: VecLockPush(ctx->current_u);
514: if (F) {
515: if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
516: ctx->current_f = F;
517: ctx->current_f_allocated = PETSC_FALSE;
518: } else if (!ctx->current_f_allocated) {
519: MatCreateVecs(J,NULL,&ctx->current_f);
521: ctx->current_f_allocated = PETSC_TRUE;
522: }
523: if (!ctx->w) {
524: VecDuplicate(ctx->current_u,&ctx->w);
525: }
526: J->assembled = PETSC_TRUE;
527: return(0);
528: }
530: typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
532: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
533: {
534: MatMFFD ctx = (MatMFFD)J->data;
537: ctx->checkh = fun;
538: ctx->checkhctx = ectx;
539: return(0);
540: }
542: /*@C
543: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
544: MatMFFD options in the database.
546: Collective on Mat
548: Input Parameter:
549: + A - the Mat context
550: - prefix - the prefix to prepend to all option names
552: Notes:
553: A hyphen (-) must NOT be given at the beginning of the prefix name.
554: The first character of all runtime options is AUTOMATICALLY the hyphen.
556: Level: advanced
558: .keywords: SNES, matrix-free, parameters
560: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
561: @*/
562: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
564: {
565: MatMFFD mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;
571: PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
572: return(0);
573: }
575: static PetscErrorCode MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
576: {
577: MatMFFD mfctx = (MatMFFD)mat->data;
579: PetscBool flg;
580: char ftype[256];
585: PetscObjectOptionsBegin((PetscObject)mfctx);
586: PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
587: if (flg) {
588: MatMFFDSetType(mat,ftype);
589: }
591: PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
592: PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
594: flg = PETSC_FALSE;
595: PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
596: if (flg) {
597: MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
598: }
599: if (mfctx->ops->setfromoptions) {
600: (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
601: }
602: PetscOptionsEnd();
603: return(0);
604: }
606: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
607: {
608: MatMFFD ctx = (MatMFFD)mat->data;
611: ctx->recomputeperiod = period;
612: return(0);
613: }
615: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
616: {
617: MatMFFD ctx = (MatMFFD)mat->data;
620: ctx->func = func;
621: ctx->funcctx = funcctx;
622: return(0);
623: }
625: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
626: {
627: MatMFFD ctx = (MatMFFD)mat->data;
630: if (error != PETSC_DEFAULT) ctx->error_rel = error;
631: return(0);
632: }
634: static PetscErrorCode MatMissingDiagonal_MFFD(Mat A,PetscBool *missing,PetscInt *d)
635: {
637: *missing = PETSC_FALSE;
638: return(0);
639: }
641: /*MC
642: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
644: Level: advanced
646: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),
647: MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
648: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
649: MatMFFDGetH(),
650: M*/
651: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
652: {
653: MatMFFD mfctx;
657: MatMFFDInitializePackage();
659: PetscHeaderCreate(mfctx,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),MatDestroy_MFFD,MatView_MFFD);
661: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
662: mfctx->recomputeperiod = 1;
663: mfctx->count = 0;
664: mfctx->currenth = 0.0;
665: mfctx->historyh = NULL;
666: mfctx->ncurrenth = 0;
667: mfctx->maxcurrenth = 0;
668: ((PetscObject)mfctx)->type_name = 0;
670: mfctx->vshift = 0.0;
671: mfctx->vscale = 1.0;
673: /*
674: Create the empty data structure to contain compute-h routines.
675: These will be filled in below from the command line options or
676: a later call with MatMFFDSetType() or if that is not called
677: then it will default in the first use of MatMult_MFFD()
678: */
679: mfctx->ops->compute = 0;
680: mfctx->ops->destroy = 0;
681: mfctx->ops->view = 0;
682: mfctx->ops->setfromoptions = 0;
683: mfctx->hctx = 0;
685: mfctx->func = 0;
686: mfctx->funcctx = 0;
687: mfctx->w = NULL;
689: A->data = mfctx;
691: A->ops->mult = MatMult_MFFD;
692: A->ops->destroy = MatDestroy_MFFD;
693: A->ops->view = MatView_MFFD;
694: A->ops->assemblyend = MatAssemblyEnd_MFFD;
695: A->ops->scale = MatScale_MFFD;
696: A->ops->shift = MatShift_MFFD;
697: A->ops->diagonalscale = MatDiagonalScale_MFFD;
698: A->ops->diagonalset = MatDiagonalSet_MFFD;
699: A->ops->setfromoptions = MatSetFromOptions_MFFD;
700: A->ops->missingdiagonal = MatMissingDiagonal_MFFD;
701: A->assembled = PETSC_TRUE;
703: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
704: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
705: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
706: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
707: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
708: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
709: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
710: PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
712: mfctx->mat = A;
714: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
715: return(0);
716: }
718: /*@
719: MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
721: Collective on Vec
723: Input Parameters:
724: + comm - MPI communicator
725: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
726: This value should be the same as the local size used in creating the
727: y vector for the matrix-vector product y = Ax.
728: . n - This value should be the same as the local size used in creating the
729: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
730: calculated if N is given) For square matrices n is almost always m.
731: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
732: - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
735: Output Parameter:
736: . J - the matrix-free matrix
738: Options Database Keys: call MatSetFromOptions() to trigger these
739: + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
740: - -mat_mffd_err - square root of estimated relative error in function evaluation
741: - -mat_mffd_period - how often h is recomputed, defaults to 1, everytime
744: Level: advanced
746: Notes:
747: The matrix-free matrix context merely contains the function pointers
748: and work space for performing finite difference approximations of
749: Jacobian-vector products, F'(u)*a,
751: The default code uses the following approach to compute h
753: .vb
754: F'(u)*a = [F(u+h*a) - F(u)]/h where
755: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
756: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
757: where
758: error_rel = square root of relative error in function evaluation
759: umin = minimum iterate parameter
760: .ve
762: You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
763: preconditioner matrix
765: The user can set the error_rel via MatMFFDSetFunctionError() and
766: umin via MatMFFDDSSetUmin(); see Users-Manual: ch_snes for details.
768: The user should call MatDestroy() when finished with the matrix-free
769: matrix context.
771: Options Database Keys:
772: + -mat_mffd_err <error_rel> - Sets error_rel
773: . -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
774: - -mat_mffd_check_positivity
776: .keywords: default, matrix-free, create, matrix
778: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
779: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
780: MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
782: @*/
783: PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
784: {
788: MatCreate(comm,J);
789: MatSetSizes(*J,m,n,M,N);
790: MatSetType(*J,MATMFFD);
791: MatSetUp(*J);
792: return(0);
793: }
795: /*@
796: MatMFFDGetH - Gets the last value that was used as the differencing
797: parameter.
799: Not Collective
801: Input Parameters:
802: . mat - the matrix obtained with MatCreateSNESMF()
804: Output Paramter:
805: . h - the differencing step size
807: Level: advanced
809: .keywords: SNES, matrix-free, parameters
811: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
812: @*/
813: PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
814: {
815: MatMFFD ctx = (MatMFFD)mat->data;
817: PetscBool match;
822: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
823: if (!match) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
825: *h = ctx->currenth;
826: return(0);
827: }
829: /*@C
830: MatMFFDSetFunction - Sets the function used in applying the matrix free.
832: Logically Collective on Mat
834: Input Parameters:
835: + mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
836: . func - the function to use
837: - funcctx - optional function context passed to function
839: Calling Sequence of func:
840: $ func (void *funcctx, Vec x, Vec f)
842: + funcctx - user provided context
843: . x - input vector
844: - f - computed output function
846: Level: advanced
848: Notes:
849: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
850: matrix inside your compute Jacobian routine
852: If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
854: .keywords: SNES, matrix-free, function
856: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
857: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
858: @*/
859: PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
860: {
865: PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
866: return(0);
867: }
869: /*@C
870: MatMFFDSetFunctioni - Sets the function for a single component
872: Logically Collective on Mat
874: Input Parameters:
875: + mat - the matrix free matrix created via MatCreateSNESMF()
876: - funci - the function to use
878: Level: advanced
880: Notes:
881: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
882: matrix inside your compute Jacobian routine.
883: This function is necessary to compute the diagonal of the matrix.
884: funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.
886: .keywords: SNES, matrix-free, function
888: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
890: @*/
891: PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
892: {
897: PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
898: return(0);
899: }
901: /*@C
902: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
904: Logically Collective on Mat
906: Input Parameters:
907: + mat - the matrix free matrix created via MatCreateSNESMF()
908: - func - the function to use
910: Level: advanced
912: Notes:
913: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
914: matrix inside your compute Jacobian routine.
915: This function is necessary to compute the diagonal of the matrix.
918: .keywords: SNES, matrix-free, function
920: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
921: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
922: @*/
923: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
924: {
929: PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
930: return(0);
931: }
933: /*@
934: MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime
936: Logically Collective on Mat
938: Input Parameters:
939: + mat - the matrix free matrix created via MatCreateSNESMF()
940: - period - 1 for everytime, 2 for every second etc
942: Options Database Keys:
943: + -mat_mffd_period <period>
945: Level: advanced
948: .keywords: SNES, matrix-free, parameters
950: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
951: MatMFFDSetHHistory(), MatMFFDResetHHistory()
952: @*/
953: PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
954: {
960: PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
961: return(0);
962: }
964: /*@
965: MatMFFDSetFunctionError - Sets the error_rel for the approximation of
966: matrix-vector products using finite differences.
968: Logically Collective on Mat
970: Input Parameters:
971: + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
972: - error_rel - relative error (should be set to the square root of
973: the relative error in the function evaluations)
975: Options Database Keys:
976: + -mat_mffd_err <error_rel> - Sets error_rel
978: Level: advanced
980: Notes:
981: The default matrix-free matrix-vector product routine computes
982: .vb
983: F'(u)*a = [F(u+h*a) - F(u)]/h where
984: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
985: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
986: .ve
988: .keywords: SNES, matrix-free, parameters
990: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
991: MatMFFDSetHHistory(), MatMFFDResetHHistory()
992: @*/
993: PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
994: {
1000: PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1001: return(0);
1002: }
1004: /*@
1005: MatMFFDSetHHistory - Sets an array to collect a history of the
1006: differencing values (h) computed for the matrix-free product.
1008: Logically Collective on Mat
1010: Input Parameters:
1011: + J - the matrix-free matrix context
1012: . histroy - space to hold the history
1013: - nhistory - number of entries in history, if more entries are generated than
1014: nhistory, then the later ones are discarded
1016: Level: advanced
1018: Notes:
1019: Use MatMFFDResetHHistory() to reset the history counter and collect
1020: a new batch of differencing parameters, h.
1022: .keywords: SNES, matrix-free, h history, differencing history
1024: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1025: MatMFFDResetHHistory(), MatMFFDSetFunctionError()
1027: @*/
1028: PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1029: {
1030: MatMFFD ctx = (MatMFFD)J->data;
1032: PetscBool match;
1038: PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1039: if (!match) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1040: ctx->historyh = history;
1041: ctx->maxcurrenth = nhistory;
1042: ctx->currenth = 0.;
1043: return(0);
1044: }
1046: /*@
1047: MatMFFDResetHHistory - Resets the counter to zero to begin
1048: collecting a new set of differencing histories.
1050: Logically Collective on Mat
1052: Input Parameters:
1053: . J - the matrix-free matrix context
1055: Level: advanced
1057: Notes:
1058: Use MatMFFDSetHHistory() to create the original history counter.
1060: .keywords: SNES, matrix-free, h history, differencing history
1062: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1063: MatMFFDSetHHistory(), MatMFFDSetFunctionError()
1065: @*/
1066: PetscErrorCode MatMFFDResetHHistory(Mat J)
1067: {
1072: PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1073: return(0);
1074: }
1076: /*@
1077: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1078: Jacobian are computed
1080: Logically Collective on Mat
1082: Input Parameters:
1083: + J - the MatMFFD matrix
1084: . U - the vector
1085: - F - (optional) vector that contains F(u) if it has been already computed
1087: Notes:
1088: This is rarely used directly
1090: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
1091: point during the first MatMult() after each call to MatMFFDSetBase().
1093: Level: advanced
1095: @*/
1096: PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
1097: {
1104: PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1105: return(0);
1106: }
1108: /*@C
1109: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1110: it to satisfy some criteria
1112: Logically Collective on Mat
1114: Input Parameters:
1115: + J - the MatMFFD matrix
1116: . fun - the function that checks h
1117: - ctx - any context needed by the function
1119: Options Database Keys:
1120: . -mat_mffd_check_positivity
1122: Level: advanced
1124: Notes:
1125: For example, MatMFFDCheckPositivity() insures that all entries
1126: of U + h*a are non-negative
1128: The function you provide is called after the default h has been computed and allows you to
1129: modify it.
1131: .seealso: MatMFFDCheckPositivity()
1132: @*/
1133: PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1134: {
1139: PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1140: return(0);
1141: }
1143: /*@
1144: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1145: zero, decreases h until this is satisfied.
1147: Logically Collective on Vec
1149: Input Parameters:
1150: + U - base vector that is added to
1151: . a - vector that is added
1152: . h - scaling factor on a
1153: - dummy - context variable (unused)
1155: Options Database Keys:
1156: . -mat_mffd_check_positivity
1158: Level: advanced
1160: Notes:
1161: This is rarely used directly, rather it is passed as an argument to
1162: MatMFFDSetCheckh()
1164: .seealso: MatMFFDSetCheckh()
1165: @*/
1166: PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1167: {
1168: PetscReal val, minval;
1169: PetscScalar *u_vec, *a_vec;
1171: PetscInt i,n;
1172: MPI_Comm comm;
1178: PetscObjectGetComm((PetscObject)U,&comm);
1179: VecGetArray(U,&u_vec);
1180: VecGetArray(a,&a_vec);
1181: VecGetLocalSize(U,&n);
1182: minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1183: for (i=0; i<n; i++) {
1184: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1185: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1186: if (val < minval) minval = val;
1187: }
1188: }
1189: VecRestoreArray(U,&u_vec);
1190: VecRestoreArray(a,&a_vec);
1191: MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1192: if (val <= PetscAbsScalar(*h)) {
1193: PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1194: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1195: else *h = -0.99*val;
1196: }
1197: return(0);
1198: }