Actual source code: mffd.c
petsc-3.11.4 2019-09-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 MatInitializePackage().
36: Level: developer
38: .keywords: Vec, initialize, package
39: .seealso: PetscInitialize()
40: @*/
41: PetscErrorCode MatMFFDInitializePackage(void)
42: {
43: char logList[256];
44: PetscBool opt,pkg;
48: if (MatMFFDPackageInitialized) return(0);
49: MatMFFDPackageInitialized = PETSC_TRUE;
50: /* Register Classes */
51: PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);
52: /* Register Constructors */
53: MatMFFDRegisterAll();
54: /* Register Events */
55: PetscLogEventRegister("MatMult MF",MATMFFD_CLASSID,&MATMFFD_Mult);
56: /* Process info exclusions */
57: PetscOptionsGetString(NULL,NULL,"-info_exclude",logList,sizeof(logList),&opt);
58: if (opt) {
59: PetscStrInList("matmffd",logList,',',&pkg);
60: if (pkg) {PetscInfoDeactivateClass(MATMFFD_CLASSID);}
61: }
62: /* Process summary exclusions */
63: PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
64: if (opt) {
65: PetscStrInList("matmffd",logList,',',&pkg);
66: if (pkg) {PetscLogEventExcludeClass(MATMFFD_CLASSID);}
67: }
68: /* Register package finalizer */
69: PetscRegisterFinalize(MatMFFDFinalizePackage);
70: return(0);
71: }
73: /*@C
74: MatMFFDSetType - Sets the method that is used to compute the
75: differencing parameter for finite differene matrix-free formulations.
77: Input Parameters:
78: + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
79: or MatSetType(mat,MATMFFD);
80: - ftype - the type requested, either MATMFFD_WP or MATMFFD_DS
82: Level: advanced
84: Notes:
85: For example, such routines can compute h for use in
86: Jacobian-vector products of the form
88: F(x+ha) - F(x)
89: F'(u)a ~= ----------------
90: h
92: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction(), MatCreateMFFD()
93: @*/
94: PetscErrorCode MatMFFDSetType(Mat mat,MatMFFDType ftype)
95: {
96: PetscErrorCode ierr,(*r)(MatMFFD);
97: MatMFFD ctx = (MatMFFD)mat->data;
98: PetscBool match;
104: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
105: if (!match) return(0);
107: /* already set, so just return */
108: PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
109: if (match) return(0);
111: /* destroy the old one if it exists */
112: if (ctx->ops->destroy) {
113: (*ctx->ops->destroy)(ctx);
114: }
116: PetscFunctionListFind(MatMFFDList,ftype,&r);
117: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
118: (*r)(ctx);
119: PetscObjectChangeTypeName((PetscObject)ctx,ftype);
120: return(0);
121: }
123: static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);
125: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
126: static PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
127: {
128: MatMFFD ctx = (MatMFFD)mat->data;
131: ctx->funcisetbase = func;
132: /* allow users to compose their own getdiagonal and allow MatHasOperation
133: to return false if the two functions pointers are not set */
134: if (!mat->ops->getdiagonal && func) {
135: mat->ops->getdiagonal = MatGetDiagonal_MFFD;
136: }
137: return(0);
138: }
140: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
141: static PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
142: {
143: MatMFFD ctx = (MatMFFD)mat->data;
146: ctx->funci = funci;
147: /* allow users to compose their own getdiagonal and allow MatHasOperation
148: to return false if the two functions pointers are not set */
149: if (!mat->ops->getdiagonal && funci) {
150: mat->ops->getdiagonal = MatGetDiagonal_MFFD;
151: }
152: return(0);
153: }
155: static PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
156: {
157: MatMFFD ctx = (MatMFFD)J->data;
160: ctx->ncurrenth = 0;
161: return(0);
162: }
164: /*@C
165: MatMFFDRegister - Adds a method to the MatMFFD registry.
167: Not Collective
169: Input Parameters:
170: + name_solver - name of a new user-defined compute-h module
171: - routine_create - routine to create method context
173: Level: developer
175: Notes:
176: MatMFFDRegister() may be called multiple times to add several user-defined solvers.
178: Sample usage:
179: .vb
180: MatMFFDRegister("my_h",MyHCreate);
181: .ve
183: Then, your solver can be chosen with the procedural interface via
184: $ MatMFFDSetType(mfctx,"my_h")
185: or at runtime via the option
186: $ -mat_mffd_type my_h
188: .keywords: MatMFFD, register
190: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
191: @*/
192: PetscErrorCode MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
193: {
197: MatInitializePackage();
198: PetscFunctionListAdd(&MatMFFDList,sname,function);
199: return(0);
200: }
202: /* ----------------------------------------------------------------------------------------*/
203: static PetscErrorCode MatDestroy_MFFD(Mat mat)
204: {
206: MatMFFD ctx = (MatMFFD)mat->data;
209: VecDestroy(&ctx->w);
210: VecDestroy(&ctx->drscale);
211: VecDestroy(&ctx->dlscale);
212: VecDestroy(&ctx->dshift);
213: VecDestroy(&ctx->dshiftw);
214: VecDestroy(&ctx->current_u);
215: if (ctx->current_f_allocated) {
216: VecDestroy(&ctx->current_f);
217: }
218: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
219: PetscHeaderDestroy(&ctx);
220: mat->data = 0;
222: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
223: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
224: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
225: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
226: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
227: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
228: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
229: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
230: return(0);
231: }
233: /*
234: MatMFFDView_MFFD - Views matrix-free parameters.
236: */
237: static PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
238: {
240: MatMFFD ctx = (MatMFFD)J->data;
241: PetscBool iascii, viewbase, viewfunction;
242: const char *prefix;
245: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
246: if (iascii) {
247: PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
248: PetscViewerASCIIPushTab(viewer);
249: PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
250: if (!((PetscObject)ctx)->type_name) {
251: PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
252: } else {
253: PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
254: }
255: #if defined(PETSC_USE_COMPLEX)
256: if (ctx->usecomplex) {
257: PetscViewerASCIIPrintf(viewer,"Using Lyness complex number trick to compute the matrix-vector product\n");
258: }
259: #endif
260: if (ctx->ops->view) {
261: (*ctx->ops->view)(ctx,viewer);
262: }
263: PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);
265: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_base", &viewbase);
266: if (viewbase) {
267: PetscViewerASCIIPrintf(viewer, "Base:\n");
268: VecView(ctx->current_u, viewer);
269: }
270: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_function", &viewfunction);
271: if (viewfunction) {
272: PetscViewerASCIIPrintf(viewer, "Function:\n");
273: VecView(ctx->current_f, viewer);
274: }
275: PetscViewerASCIIPopTab(viewer);
276: }
277: return(0);
278: }
280: /*
281: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
282: allows the user to indicate the beginning of a new linear solve by calling
283: MatAssemblyXXX() on the matrix free matrix. This then allows the
284: MatCreateMFFD_WP() to properly compute ||U|| only the first time
285: in the linear solver rather than every time.
287: This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
288: it must be labeled as PETSC_EXTERN
289: */
290: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
291: {
293: MatMFFD j = (MatMFFD)J->data;
296: MatMFFDResetHHistory(J);
297: j->vshift = 0.0;
298: j->vscale = 1.0;
299: return(0);
300: }
302: /*
303: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
305: y ~= (F(u + ha) - F(u))/h,
306: where F = nonlinear function, as set by SNESSetFunction()
307: u = current iterate
308: h = difference interval
309: */
310: static PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
311: {
312: MatMFFD ctx = (MatMFFD)mat->data;
313: PetscScalar h;
314: Vec w,U,F;
316: PetscBool zeroa;
319: 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");
320: /* We log matrix-free matrix-vector products separately, so that we can
321: separate the performance monitoring from the cases that use conventional
322: storage. We may eventually modify event logging to associate events
323: with particular objects, hence alleviating the more general problem. */
324: PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);
326: w = ctx->w;
327: U = ctx->current_u;
328: F = ctx->current_f;
329: /*
330: Compute differencing parameter
331: */
332: if (!((PetscObject)ctx)->type_name) {
333: MatMFFDSetType(mat,MATMFFD_WP);
334: MatSetFromOptions(mat);
335: }
336: (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
337: if (zeroa) {
338: VecSet(y,0.0);
339: return(0);
340: }
342: if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
343: if (ctx->checkh) {
344: (*ctx->checkh)(ctx->checkhctx,U,a,&h);
345: }
347: /* keep a record of the current differencing parameter h */
348: ctx->currenth = h;
349: #if defined(PETSC_USE_COMPLEX)
350: PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
351: #else
352: PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
353: #endif
354: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
355: ctx->historyh[ctx->ncurrenth] = h;
356: }
357: ctx->ncurrenth++;
359: #if defined(PETSC_USE_COMPLEX)
360: if (ctx->usecomplex) h = PETSC_i*h;
361: #endif
362:
363: /* w = u + ha */
364: if (ctx->drscale) {
365: VecPointwiseMult(ctx->drscale,a,U);
366: VecAYPX(U,h,w);
367: } else {
368: VecWAXPY(w,h,a,U);
369: }
371: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
372: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
373: (*ctx->func)(ctx->funcctx,U,F);
374: }
375: (*ctx->func)(ctx->funcctx,w,y);
377: #if defined(PETSC_USE_COMPLEX)
378: if (ctx->usecomplex) {
379: VecImaginaryPart(y);
380: h = PetscImaginaryPart(h);
381: } else {
382: VecAXPY(y,-1.0,F);
383: }
384: #else
385: VecAXPY(y,-1.0,F);
386: #endif
387: VecScale(y,1.0/h);
389: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
390: VecAXPBY(y,ctx->vshift,ctx->vscale,a);
391: }
392: if (ctx->dlscale) {
393: VecPointwiseMult(y,ctx->dlscale,y);
394: }
395: if (ctx->dshift) {
396: if (!ctx->dshiftw) {
397: VecDuplicate(y,&ctx->dshiftw);
398: }
399: VecPointwiseMult(ctx->dshift,a,ctx->dshiftw);
400: VecAXPY(y,1.0,ctx->dshiftw);
401: }
403: if (mat->nullsp) {MatNullSpaceRemove(mat->nullsp,y);}
405: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
406: return(0);
407: }
409: /*
410: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
412: y ~= (F(u + ha) - F(u))/h,
413: where F = nonlinear function, as set by SNESSetFunction()
414: u = current iterate
415: h = difference interval
416: */
417: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
418: {
419: MatMFFD ctx = (MatMFFD)mat->data;
420: PetscScalar h,*aa,*ww,v;
421: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
422: Vec w,U;
424: PetscInt i,rstart,rend;
427: if (!ctx->func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunction() first");
428: if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioni() first");
429: if (!ctx->funcisetbase) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioniBase() first");
430: w = ctx->w;
431: U = ctx->current_u;
432: (*ctx->func)(ctx->funcctx,U,a);
433: (*ctx->funcisetbase)(ctx->funcctx,U);
434: VecCopy(U,w);
436: VecGetOwnershipRange(a,&rstart,&rend);
437: VecGetArray(a,&aa);
438: for (i=rstart; i<rend; i++) {
439: VecGetArray(w,&ww);
440: h = ww[i-rstart];
441: if (h == 0.0) h = 1.0;
442: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
443: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
444: h *= epsilon;
446: ww[i-rstart] += h;
447: VecRestoreArray(w,&ww);
448: (*ctx->funci)(ctx->funcctx,i,w,&v);
449: aa[i-rstart] = (v - aa[i-rstart])/h;
451: /* possibly shift and scale result */
452: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
453: aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
454: }
456: VecGetArray(w,&ww);
457: ww[i-rstart] -= h;
458: VecRestoreArray(w,&ww);
459: }
460: VecRestoreArray(a,&aa);
461: return(0);
462: }
464: static PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
465: {
466: MatMFFD aij = (MatMFFD)mat->data;
470: if (ll && !aij->dlscale) {
471: VecDuplicate(ll,&aij->dlscale);
472: }
473: if (rr && !aij->drscale) {
474: VecDuplicate(rr,&aij->drscale);
475: }
476: if (ll) {
477: VecCopy(ll,aij->dlscale);
478: }
479: if (rr) {
480: VecCopy(rr,aij->drscale);
481: }
482: return(0);
483: }
485: static PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
486: {
487: MatMFFD aij = (MatMFFD)mat->data;
491: if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
492: if (!aij->dshift) {
493: VecDuplicate(ll,&aij->dshift);
494: }
495: VecAXPY(aij->dshift,1.0,ll);
496: return(0);
497: }
499: static PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
500: {
501: MatMFFD shell = (MatMFFD)Y->data;
504: shell->vshift += a;
505: return(0);
506: }
508: static PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
509: {
510: MatMFFD shell = (MatMFFD)Y->data;
513: shell->vscale *= a;
514: return(0);
515: }
517: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
518: {
520: MatMFFD ctx = (MatMFFD)J->data;
523: MatMFFDResetHHistory(J);
524: if (!ctx->current_u) {
525: VecDuplicate(U,&ctx->current_u);
526: VecLockReadPush(ctx->current_u);
527: }
528: VecLockReadPop(ctx->current_u);
529: VecCopy(U,ctx->current_u);
530: VecLockReadPush(ctx->current_u);
531: if (F) {
532: if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
533: ctx->current_f = F;
534: ctx->current_f_allocated = PETSC_FALSE;
535: } else if (!ctx->current_f_allocated) {
536: MatCreateVecs(J,NULL,&ctx->current_f);
538: ctx->current_f_allocated = PETSC_TRUE;
539: }
540: if (!ctx->w) {
541: VecDuplicate(ctx->current_u,&ctx->w);
542: }
543: J->assembled = PETSC_TRUE;
544: return(0);
545: }
547: typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
549: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
550: {
551: MatMFFD ctx = (MatMFFD)J->data;
554: ctx->checkh = fun;
555: ctx->checkhctx = ectx;
556: return(0);
557: }
559: /*@C
560: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
561: MatMFFD options in the database.
563: Collective on Mat
565: Input Parameter:
566: + A - the Mat context
567: - prefix - the prefix to prepend to all option names
569: Notes:
570: A hyphen (-) must NOT be given at the beginning of the prefix name.
571: The first character of all runtime options is AUTOMATICALLY the hyphen.
573: Level: advanced
575: .keywords: SNES, matrix-free, parameters
577: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
578: @*/
579: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
581: {
582: MatMFFD mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;
588: PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
589: return(0);
590: }
592: static PetscErrorCode MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
593: {
594: MatMFFD mfctx = (MatMFFD)mat->data;
596: PetscBool flg;
597: char ftype[256];
602: PetscObjectOptionsBegin((PetscObject)mfctx);
603: PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
604: if (flg) {
605: MatMFFDSetType(mat,ftype);
606: }
608: PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
609: PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
611: flg = PETSC_FALSE;
612: PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
613: if (flg) {
614: MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
615: }
616: #if defined(PETSC_USE_COMPLEX)
617: PetscOptionsBool("-mat_mffd_complex","Use Lyness complex number trick to compute the matrix-vector product","None",mfctx->usecomplex,&mfctx->usecomplex,NULL);
618: #endif
619: if (mfctx->ops->setfromoptions) {
620: (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
621: }
622: PetscOptionsEnd();
623: return(0);
624: }
626: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
627: {
628: MatMFFD ctx = (MatMFFD)mat->data;
631: ctx->recomputeperiod = period;
632: return(0);
633: }
635: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
636: {
637: MatMFFD ctx = (MatMFFD)mat->data;
640: ctx->func = func;
641: ctx->funcctx = funcctx;
642: return(0);
643: }
645: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
646: {
647: MatMFFD ctx = (MatMFFD)mat->data;
650: if (error != PETSC_DEFAULT) ctx->error_rel = error;
651: return(0);
652: }
654: static PetscErrorCode MatMissingDiagonal_MFFD(Mat A,PetscBool *missing,PetscInt *d)
655: {
657: *missing = PETSC_FALSE;
658: return(0);
659: }
661: /*MC
662: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
664: Level: advanced
666: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),
667: MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
668: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
669: MatMFFDGetH(),
670: M*/
671: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
672: {
673: MatMFFD mfctx;
677: MatMFFDInitializePackage();
679: PetscHeaderCreate(mfctx,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),MatDestroy_MFFD,MatView_MFFD);
681: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
682: mfctx->recomputeperiod = 1;
683: mfctx->count = 0;
684: mfctx->currenth = 0.0;
685: mfctx->historyh = NULL;
686: mfctx->ncurrenth = 0;
687: mfctx->maxcurrenth = 0;
688: ((PetscObject)mfctx)->type_name = 0;
690: mfctx->vshift = 0.0;
691: mfctx->vscale = 1.0;
693: /*
694: Create the empty data structure to contain compute-h routines.
695: These will be filled in below from the command line options or
696: a later call with MatMFFDSetType() or if that is not called
697: then it will default in the first use of MatMult_MFFD()
698: */
699: mfctx->ops->compute = 0;
700: mfctx->ops->destroy = 0;
701: mfctx->ops->view = 0;
702: mfctx->ops->setfromoptions = 0;
703: mfctx->hctx = 0;
705: mfctx->func = 0;
706: mfctx->funcctx = 0;
707: mfctx->w = NULL;
709: A->data = mfctx;
711: A->ops->mult = MatMult_MFFD;
712: A->ops->destroy = MatDestroy_MFFD;
713: A->ops->view = MatView_MFFD;
714: A->ops->assemblyend = MatAssemblyEnd_MFFD;
715: A->ops->scale = MatScale_MFFD;
716: A->ops->shift = MatShift_MFFD;
717: A->ops->diagonalscale = MatDiagonalScale_MFFD;
718: A->ops->diagonalset = MatDiagonalSet_MFFD;
719: A->ops->setfromoptions = MatSetFromOptions_MFFD;
720: A->ops->missingdiagonal = MatMissingDiagonal_MFFD;
721: A->assembled = PETSC_TRUE;
723: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
724: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
725: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
726: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
727: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
728: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
729: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
730: PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
732: mfctx->mat = A;
734: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
735: return(0);
736: }
738: /*@
739: MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
741: Collective on Vec
743: Input Parameters:
744: + comm - MPI communicator
745: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
746: This value should be the same as the local size used in creating the
747: y vector for the matrix-vector product y = Ax.
748: . n - This value should be the same as the local size used in creating the
749: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
750: calculated if N is given) For square matrices n is almost always m.
751: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
752: - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
755: Output Parameter:
756: . J - the matrix-free matrix
758: Options Database Keys: call MatSetFromOptions() to trigger these
759: + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
760: . -mat_mffd_err - square root of estimated relative error in function evaluation
761: . -mat_mffd_period - how often h is recomputed, defaults to 1, everytime
762: . -mat_mffd_check_positivity - possibly decrease h until U + h*a has only positive values
763: - -mat_mffd_complex - use the Lyness trick with complex numbers to compute the matrix-vector product instead of differencing
764: (requires real valued functions but that PETSc be configured for complex numbers)
767: Level: advanced
769: Notes:
770: The matrix-free matrix context merely contains the function pointers
771: and work space for performing finite difference approximations of
772: Jacobian-vector products, F'(u)*a,
774: The default code uses the following approach to compute h
776: .vb
777: F'(u)*a = [F(u+h*a) - F(u)]/h where
778: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
779: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
780: where
781: error_rel = square root of relative error in function evaluation
782: umin = minimum iterate parameter
783: .ve
785: You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
786: preconditioner matrix
788: The user can set the error_rel via MatMFFDSetFunctionError() and
789: umin via MatMFFDDSSetUmin(); see Users-Manual: ch_snes for details.
791: The user should call MatDestroy() when finished with the matrix-free
792: matrix context.
794: Options Database Keys:
795: + -mat_mffd_err <error_rel> - Sets error_rel
796: . -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
797: - -mat_mffd_check_positivity
799: .keywords: default, matrix-free, create, matrix
801: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
802: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
803: MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
805: @*/
806: PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
807: {
811: MatCreate(comm,J);
812: MatSetSizes(*J,m,n,M,N);
813: MatSetType(*J,MATMFFD);
814: MatSetUp(*J);
815: return(0);
816: }
818: /*@
819: MatMFFDGetH - Gets the last value that was used as the differencing
820: parameter.
822: Not Collective
824: Input Parameters:
825: . mat - the matrix obtained with MatCreateSNESMF()
827: Output Paramter:
828: . h - the differencing step size
830: Level: advanced
832: .keywords: SNES, matrix-free, parameters
834: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
835: @*/
836: PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
837: {
838: MatMFFD ctx = (MatMFFD)mat->data;
840: PetscBool match;
845: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
846: if (!match) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
848: *h = ctx->currenth;
849: return(0);
850: }
852: /*@C
853: MatMFFDSetFunction - Sets the function used in applying the matrix free.
855: Logically Collective on Mat
857: Input Parameters:
858: + mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
859: . func - the function to use
860: - funcctx - optional function context passed to function
862: Calling Sequence of func:
863: $ func (void *funcctx, Vec x, Vec f)
865: + funcctx - user provided context
866: . x - input vector
867: - f - computed output function
869: Level: advanced
871: Notes:
872: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
873: matrix inside your compute Jacobian routine
875: If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
877: .keywords: SNES, matrix-free, function
879: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
880: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
881: @*/
882: PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
883: {
888: PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
889: return(0);
890: }
892: /*@C
893: MatMFFDSetFunctioni - Sets the function for a single component
895: Logically Collective on Mat
897: Input Parameters:
898: + mat - the matrix free matrix created via MatCreateSNESMF()
899: - funci - the function to use
901: Level: advanced
903: Notes:
904: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
905: matrix inside your compute Jacobian routine.
906: This function is necessary to compute the diagonal of the matrix.
907: funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.
909: .keywords: SNES, matrix-free, function
911: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
913: @*/
914: PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
915: {
920: PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
921: return(0);
922: }
924: /*@C
925: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
927: Logically Collective on Mat
929: Input Parameters:
930: + mat - the matrix free matrix created via MatCreateSNESMF()
931: - func - the function to use
933: Level: advanced
935: Notes:
936: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
937: matrix inside your compute Jacobian routine.
938: This function is necessary to compute the diagonal of the matrix.
941: .keywords: SNES, matrix-free, function
943: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
944: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
945: @*/
946: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
947: {
952: PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
953: return(0);
954: }
956: /*@
957: MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime
959: Logically Collective on Mat
961: Input Parameters:
962: + mat - the matrix free matrix created via MatCreateSNESMF()
963: - period - 1 for everytime, 2 for every second etc
965: Options Database Keys:
966: + -mat_mffd_period <period>
968: Level: advanced
971: .keywords: SNES, matrix-free, parameters
973: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
974: MatMFFDSetHHistory(), MatMFFDResetHHistory()
975: @*/
976: PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
977: {
983: PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
984: return(0);
985: }
987: /*@
988: MatMFFDSetFunctionError - Sets the error_rel for the approximation of
989: matrix-vector products using finite differences.
991: Logically Collective on Mat
993: Input Parameters:
994: + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
995: - error_rel - relative error (should be set to the square root of
996: the relative error in the function evaluations)
998: Options Database Keys:
999: + -mat_mffd_err <error_rel> - Sets error_rel
1001: Level: advanced
1003: Notes:
1004: The default matrix-free matrix-vector product routine computes
1005: .vb
1006: F'(u)*a = [F(u+h*a) - F(u)]/h where
1007: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
1008: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
1009: .ve
1011: .keywords: SNES, matrix-free, parameters
1013: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
1014: MatMFFDSetHHistory(), MatMFFDResetHHistory()
1015: @*/
1016: PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
1017: {
1023: PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1024: return(0);
1025: }
1027: /*@
1028: MatMFFDSetHHistory - Sets an array to collect a history of the
1029: differencing values (h) computed for the matrix-free product.
1031: Logically Collective on Mat
1033: Input Parameters:
1034: + J - the matrix-free matrix context
1035: . histroy - space to hold the history
1036: - nhistory - number of entries in history, if more entries are generated than
1037: nhistory, then the later ones are discarded
1039: Level: advanced
1041: Notes:
1042: Use MatMFFDResetHHistory() to reset the history counter and collect
1043: a new batch of differencing parameters, h.
1045: .keywords: SNES, matrix-free, h history, differencing history
1047: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1048: MatMFFDResetHHistory(), MatMFFDSetFunctionError()
1050: @*/
1051: PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1052: {
1053: MatMFFD ctx = (MatMFFD)J->data;
1055: PetscBool match;
1061: PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1062: if (!match) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1063: ctx->historyh = history;
1064: ctx->maxcurrenth = nhistory;
1065: ctx->currenth = 0.;
1066: return(0);
1067: }
1069: /*@
1070: MatMFFDResetHHistory - Resets the counter to zero to begin
1071: collecting a new set of differencing histories.
1073: Logically Collective on Mat
1075: Input Parameters:
1076: . J - the matrix-free matrix context
1078: Level: advanced
1080: Notes:
1081: Use MatMFFDSetHHistory() to create the original history counter.
1083: .keywords: SNES, matrix-free, h history, differencing history
1085: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1086: MatMFFDSetHHistory(), MatMFFDSetFunctionError()
1088: @*/
1089: PetscErrorCode MatMFFDResetHHistory(Mat J)
1090: {
1095: PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1096: return(0);
1097: }
1099: /*@
1100: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1101: Jacobian are computed
1103: Logically Collective on Mat
1105: Input Parameters:
1106: + J - the MatMFFD matrix
1107: . U - the vector
1108: - F - (optional) vector that contains F(u) if it has been already computed
1110: Notes:
1111: This is rarely used directly
1113: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
1114: point during the first MatMult() after each call to MatMFFDSetBase().
1116: Level: advanced
1118: @*/
1119: PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
1120: {
1127: PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1128: return(0);
1129: }
1131: /*@C
1132: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1133: it to satisfy some criteria
1135: Logically Collective on Mat
1137: Input Parameters:
1138: + J - the MatMFFD matrix
1139: . fun - the function that checks h
1140: - ctx - any context needed by the function
1142: Options Database Keys:
1143: . -mat_mffd_check_positivity
1145: Level: advanced
1147: Notes:
1148: For example, MatMFFDCheckPositivity() insures that all entries
1149: of U + h*a are non-negative
1151: The function you provide is called after the default h has been computed and allows you to
1152: modify it.
1154: .seealso: MatMFFDCheckPositivity()
1155: @*/
1156: PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1157: {
1162: PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1163: return(0);
1164: }
1166: /*@
1167: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1168: zero, decreases h until this is satisfied.
1170: Logically Collective on Vec
1172: Input Parameters:
1173: + U - base vector that is added to
1174: . a - vector that is added
1175: . h - scaling factor on a
1176: - dummy - context variable (unused)
1178: Options Database Keys:
1179: . -mat_mffd_check_positivity
1181: Level: advanced
1183: Notes:
1184: This is rarely used directly, rather it is passed as an argument to
1185: MatMFFDSetCheckh()
1187: .seealso: MatMFFDSetCheckh()
1188: @*/
1189: PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1190: {
1191: PetscReal val, minval;
1192: PetscScalar *u_vec, *a_vec;
1194: PetscInt i,n;
1195: MPI_Comm comm;
1201: PetscObjectGetComm((PetscObject)U,&comm);
1202: VecGetArray(U,&u_vec);
1203: VecGetArray(a,&a_vec);
1204: VecGetLocalSize(U,&n);
1205: minval = PetscAbsScalar(*h)*PetscRealConstant(1.01);
1206: for (i=0; i<n; i++) {
1207: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1208: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1209: if (val < minval) minval = val;
1210: }
1211: }
1212: VecRestoreArray(U,&u_vec);
1213: VecRestoreArray(a,&a_vec);
1214: MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1215: if (val <= PetscAbsScalar(*h)) {
1216: PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1217: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1218: else *h = -0.99*val;
1219: }
1220: return(0);
1221: }