Actual source code: mffd.c
petsc-3.8.4 2018-03-24
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: char *className;
46: PetscBool opt;
50: if (MatMFFDPackageInitialized) return(0);
51: MatMFFDPackageInitialized = PETSC_TRUE;
52: /* Register Classes */
53: PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);
54: /* Register Constructors */
55: MatMFFDRegisterAll();
56: /* Register Events */
57: PetscLogEventRegister("MatMult MF", MATMFFD_CLASSID,&MATMFFD_Mult);
59: /* Process info exclusions */
60: PetscOptionsGetString(NULL,NULL, "-info_exclude", logList, 256, &opt);
61: if (opt) {
62: PetscStrstr(logList, "matmffd", &className);
63: if (className) {
64: PetscInfoDeactivateClass(MATMFFD_CLASSID);
65: }
66: }
67: /* Process summary exclusions */
68: PetscOptionsGetString(NULL,NULL, "-log_exclude", logList, 256, &opt);
69: if (opt) {
70: PetscStrstr(logList, "matmffd", &className);
71: if (className) {
72: PetscLogEventDeactivateClass(MATMFFD_CLASSID);
73: }
74: }
75: PetscRegisterFinalize(MatMFFDFinalizePackage);
76: return(0);
77: }
79: /*@C
80: MatMFFDSetType - Sets the method that is used to compute the
81: differencing parameter for finite differene matrix-free formulations.
83: Input Parameters:
84: + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
85: or MatSetType(mat,MATMFFD);
86: - ftype - the type requested, either MATMFFD_WP or MATMFFD_DS
88: Level: advanced
90: Notes:
91: For example, such routines can compute h for use in
92: Jacobian-vector products of the form
94: F(x+ha) - F(x)
95: F'(u)a ~= ----------------
96: h
98: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction(), MatCreateMFFD()
99: @*/
100: PetscErrorCode MatMFFDSetType(Mat mat,MatMFFDType ftype)
101: {
102: PetscErrorCode ierr,(*r)(MatMFFD);
103: MatMFFD ctx = (MatMFFD)mat->data;
104: PetscBool match;
110: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
111: if (!match) return(0);
113: /* already set, so just return */
114: PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
115: if (match) return(0);
117: /* destroy the old one if it exists */
118: if (ctx->ops->destroy) {
119: (*ctx->ops->destroy)(ctx);
120: }
122: PetscFunctionListFind(MatMFFDList,ftype,&r);
123: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
124: (*r)(ctx);
125: PetscObjectChangeTypeName((PetscObject)ctx,ftype);
126: return(0);
127: }
129: static PetscErrorCode MatGetDiagonal_MFFD(Mat,Vec);
131: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
132: static PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
133: {
134: MatMFFD ctx = (MatMFFD)mat->data;
137: ctx->funcisetbase = func;
138: /* allow users to compose their own getdiagonal and allow MatHasOperation
139: to return false if the two functions pointers are not set */
140: if (!mat->ops->getdiagonal && func) {
141: mat->ops->getdiagonal = MatGetDiagonal_MFFD;
142: }
143: return(0);
144: }
146: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
147: static PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
148: {
149: MatMFFD ctx = (MatMFFD)mat->data;
152: ctx->funci = funci;
153: /* allow users to compose their own getdiagonal and allow MatHasOperation
154: to return false if the two functions pointers are not set */
155: if (!mat->ops->getdiagonal && funci) {
156: mat->ops->getdiagonal = MatGetDiagonal_MFFD;
157: }
158: return(0);
159: }
161: static PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
162: {
163: MatMFFD ctx = (MatMFFD)J->data;
166: ctx->ncurrenth = 0;
167: return(0);
168: }
170: /*@C
171: MatMFFDRegister - Adds a method to the MatMFFD registry.
173: Not Collective
175: Input Parameters:
176: + name_solver - name of a new user-defined compute-h module
177: - routine_create - routine to create method context
179: Level: developer
181: Notes:
182: MatMFFDRegister() may be called multiple times to add several user-defined solvers.
184: Sample usage:
185: .vb
186: MatMFFDRegister("my_h",MyHCreate);
187: .ve
189: Then, your solver can be chosen with the procedural interface via
190: $ MatMFFDSetType(mfctx,"my_h")
191: or at runtime via the option
192: $ -mat_mffd_type my_h
194: .keywords: MatMFFD, register
196: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
197: @*/
198: PetscErrorCode MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
199: {
203: PetscFunctionListAdd(&MatMFFDList,sname,function);
204: return(0);
205: }
207: /* ----------------------------------------------------------------------------------------*/
208: static PetscErrorCode MatDestroy_MFFD(Mat mat)
209: {
211: MatMFFD ctx = (MatMFFD)mat->data;
214: VecDestroy(&ctx->w);
215: VecDestroy(&ctx->drscale);
216: VecDestroy(&ctx->dlscale);
217: VecDestroy(&ctx->dshift);
218: VecDestroy(&ctx->dshiftw);
219: VecDestroy(&ctx->current_u);
220: if (ctx->current_f_allocated) {
221: VecDestroy(&ctx->current_f);
222: }
223: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
224: PetscHeaderDestroy(&ctx);
225: mat->data = 0;
227: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
228: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
229: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
230: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
231: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
232: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
233: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
234: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
235: return(0);
236: }
238: /*
239: MatMFFDView_MFFD - Views matrix-free parameters.
241: */
242: static PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
243: {
245: MatMFFD ctx = (MatMFFD)J->data;
246: PetscBool iascii, viewbase, viewfunction;
247: const char *prefix;
250: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
251: if (iascii) {
252: PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
253: PetscViewerASCIIPushTab(viewer);
254: PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
255: if (!((PetscObject)ctx)->type_name) {
256: PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
257: } else {
258: PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
259: }
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: /* w = u + ha */
360: if (ctx->drscale) {
361: VecPointwiseMult(ctx->drscale,a,U);
362: VecAYPX(U,h,w);
363: } else {
364: VecWAXPY(w,h,a,U);
365: }
367: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
368: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
369: (*ctx->func)(ctx->funcctx,U,F);
370: }
371: (*ctx->func)(ctx->funcctx,w,y);
373: VecAXPY(y,-1.0,F);
374: VecScale(y,1.0/h);
376: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
377: VecAXPBY(y,ctx->vshift,ctx->vscale,a);
378: }
379: if (ctx->dlscale) {
380: VecPointwiseMult(y,ctx->dlscale,y);
381: }
382: if (ctx->dshift) {
383: if (!ctx->dshiftw) {
384: VecDuplicate(y,&ctx->dshiftw);
385: }
386: VecPointwiseMult(ctx->dshift,a,ctx->dshiftw);
387: VecAXPY(y,1.0,ctx->dshiftw);
388: }
390: if (mat->nullsp) {MatNullSpaceRemove(mat->nullsp,y);}
392: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
393: return(0);
394: }
396: /*
397: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
399: y ~= (F(u + ha) - F(u))/h,
400: where F = nonlinear function, as set by SNESSetFunction()
401: u = current iterate
402: h = difference interval
403: */
404: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
405: {
406: MatMFFD ctx = (MatMFFD)mat->data;
407: PetscScalar h,*aa,*ww,v;
408: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
409: Vec w,U;
411: PetscInt i,rstart,rend;
414: if (!ctx->func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunction() first");
415: if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioni() first");
416: if (!ctx->funcisetbase) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Requires calling MatMFFDSetFunctioniBase() first");
417: w = ctx->w;
418: U = ctx->current_u;
419: (*ctx->func)(ctx->funcctx,U,a);
420: (*ctx->funcisetbase)(ctx->funcctx,U);
421: VecCopy(U,w);
423: VecGetOwnershipRange(a,&rstart,&rend);
424: VecGetArray(a,&aa);
425: for (i=rstart; i<rend; i++) {
426: VecGetArray(w,&ww);
427: h = ww[i-rstart];
428: if (h == 0.0) h = 1.0;
429: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
430: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
431: h *= epsilon;
433: ww[i-rstart] += h;
434: VecRestoreArray(w,&ww);
435: (*ctx->funci)(ctx->funcctx,i,w,&v);
436: aa[i-rstart] = (v - aa[i-rstart])/h;
438: /* possibly shift and scale result */
439: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
440: aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
441: }
443: VecGetArray(w,&ww);
444: ww[i-rstart] -= h;
445: VecRestoreArray(w,&ww);
446: }
447: VecRestoreArray(a,&aa);
448: return(0);
449: }
451: static PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
452: {
453: MatMFFD aij = (MatMFFD)mat->data;
457: if (ll && !aij->dlscale) {
458: VecDuplicate(ll,&aij->dlscale);
459: }
460: if (rr && !aij->drscale) {
461: VecDuplicate(rr,&aij->drscale);
462: }
463: if (ll) {
464: VecCopy(ll,aij->dlscale);
465: }
466: if (rr) {
467: VecCopy(rr,aij->drscale);
468: }
469: return(0);
470: }
472: static PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
473: {
474: MatMFFD aij = (MatMFFD)mat->data;
478: if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
479: if (!aij->dshift) {
480: VecDuplicate(ll,&aij->dshift);
481: }
482: VecAXPY(aij->dshift,1.0,ll);
483: return(0);
484: }
486: static PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
487: {
488: MatMFFD shell = (MatMFFD)Y->data;
491: shell->vshift += a;
492: return(0);
493: }
495: static PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
496: {
497: MatMFFD shell = (MatMFFD)Y->data;
500: shell->vscale *= a;
501: return(0);
502: }
504: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
505: {
507: MatMFFD ctx = (MatMFFD)J->data;
510: MatMFFDResetHHistory(J);
511: if (!ctx->current_u) {
512: VecDuplicate(U,&ctx->current_u);
513: VecLockPush(ctx->current_u);
514: }
515: VecLockPop(ctx->current_u);
516: VecCopy(U,ctx->current_u);
517: VecLockPush(ctx->current_u);
518: if (F) {
519: if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
520: ctx->current_f = F;
521: ctx->current_f_allocated = PETSC_FALSE;
522: } else if (!ctx->current_f_allocated) {
523: MatCreateVecs(J,NULL,&ctx->current_f);
525: ctx->current_f_allocated = PETSC_TRUE;
526: }
527: if (!ctx->w) {
528: VecDuplicate(ctx->current_u,&ctx->w);
529: }
530: J->assembled = PETSC_TRUE;
531: return(0);
532: }
534: typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
536: static PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
537: {
538: MatMFFD ctx = (MatMFFD)J->data;
541: ctx->checkh = fun;
542: ctx->checkhctx = ectx;
543: return(0);
544: }
546: /*@C
547: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
548: MatMFFD options in the database.
550: Collective on Mat
552: Input Parameter:
553: + A - the Mat context
554: - prefix - the prefix to prepend to all option names
556: Notes:
557: A hyphen (-) must NOT be given at the beginning of the prefix name.
558: The first character of all runtime options is AUTOMATICALLY the hyphen.
560: Level: advanced
562: .keywords: SNES, matrix-free, parameters
564: .seealso: MatSetFromOptions(), MatCreateSNESMF(), MatCreateMFFD()
565: @*/
566: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
568: {
569: MatMFFD mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;
575: PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
576: return(0);
577: }
579: static PetscErrorCode MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
580: {
581: MatMFFD mfctx = (MatMFFD)mat->data;
583: PetscBool flg;
584: char ftype[256];
589: PetscObjectOptionsBegin((PetscObject)mfctx);
590: PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
591: if (flg) {
592: MatMFFDSetType(mat,ftype);
593: }
595: PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
596: PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
598: flg = PETSC_FALSE;
599: PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
600: if (flg) {
601: MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
602: }
603: if (mfctx->ops->setfromoptions) {
604: (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
605: }
606: PetscOptionsEnd();
607: return(0);
608: }
610: static PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
611: {
612: MatMFFD ctx = (MatMFFD)mat->data;
615: ctx->recomputeperiod = period;
616: return(0);
617: }
619: static PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
620: {
621: MatMFFD ctx = (MatMFFD)mat->data;
624: ctx->func = func;
625: ctx->funcctx = funcctx;
626: return(0);
627: }
629: static PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
630: {
631: MatMFFD ctx = (MatMFFD)mat->data;
634: if (error != PETSC_DEFAULT) ctx->error_rel = error;
635: return(0);
636: }
638: static PetscErrorCode MatMissingDiagonal_MFFD(Mat A,PetscBool *missing,PetscInt *d)
639: {
641: *missing = PETSC_FALSE;
642: return(0);
643: }
645: /*MC
646: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
648: Level: advanced
650: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction(), MatMFFDSetType(),
651: MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
652: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
653: MatMFFDGetH(),
654: M*/
655: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
656: {
657: MatMFFD mfctx;
661: MatMFFDInitializePackage();
663: PetscHeaderCreate(mfctx,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),MatDestroy_MFFD,MatView_MFFD);
665: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
666: mfctx->recomputeperiod = 1;
667: mfctx->count = 0;
668: mfctx->currenth = 0.0;
669: mfctx->historyh = NULL;
670: mfctx->ncurrenth = 0;
671: mfctx->maxcurrenth = 0;
672: ((PetscObject)mfctx)->type_name = 0;
674: mfctx->vshift = 0.0;
675: mfctx->vscale = 1.0;
677: /*
678: Create the empty data structure to contain compute-h routines.
679: These will be filled in below from the command line options or
680: a later call with MatMFFDSetType() or if that is not called
681: then it will default in the first use of MatMult_MFFD()
682: */
683: mfctx->ops->compute = 0;
684: mfctx->ops->destroy = 0;
685: mfctx->ops->view = 0;
686: mfctx->ops->setfromoptions = 0;
687: mfctx->hctx = 0;
689: mfctx->func = 0;
690: mfctx->funcctx = 0;
691: mfctx->w = NULL;
693: A->data = mfctx;
695: A->ops->mult = MatMult_MFFD;
696: A->ops->destroy = MatDestroy_MFFD;
697: A->ops->view = MatView_MFFD;
698: A->ops->assemblyend = MatAssemblyEnd_MFFD;
699: A->ops->scale = MatScale_MFFD;
700: A->ops->shift = MatShift_MFFD;
701: A->ops->diagonalscale = MatDiagonalScale_MFFD;
702: A->ops->diagonalset = MatDiagonalSet_MFFD;
703: A->ops->setfromoptions = MatSetFromOptions_MFFD;
704: A->ops->missingdiagonal = MatMissingDiagonal_MFFD;
705: A->assembled = PETSC_TRUE;
707: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
708: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
709: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
710: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
711: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
712: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
713: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
714: PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
716: mfctx->mat = A;
718: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
719: return(0);
720: }
722: /*@
723: MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
725: Collective on Vec
727: Input Parameters:
728: + comm - MPI communicator
729: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
730: This value should be the same as the local size used in creating the
731: y vector for the matrix-vector product y = Ax.
732: . n - This value should be the same as the local size used in creating the
733: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
734: calculated if N is given) For square matrices n is almost always m.
735: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
736: - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
739: Output Parameter:
740: . J - the matrix-free matrix
742: Options Database Keys: call MatSetFromOptions() to trigger these
743: + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
744: - -mat_mffd_err - square root of estimated relative error in function evaluation
745: - -mat_mffd_period - how often h is recomputed, defaults to 1, everytime
748: Level: advanced
750: Notes:
751: The matrix-free matrix context merely contains the function pointers
752: and work space for performing finite difference approximations of
753: Jacobian-vector products, F'(u)*a,
755: The default code uses the following approach to compute h
757: .vb
758: F'(u)*a = [F(u+h*a) - F(u)]/h where
759: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
760: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
761: where
762: error_rel = square root of relative error in function evaluation
763: umin = minimum iterate parameter
764: .ve
766: You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
767: preconditioner matrix
769: The user can set the error_rel via MatMFFDSetFunctionError() and
770: umin via MatMFFDDSSetUmin(); see Users-Manual: ch_snes for details.
772: The user should call MatDestroy() when finished with the matrix-free
773: matrix context.
775: Options Database Keys:
776: + -mat_mffd_err <error_rel> - Sets error_rel
777: . -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
778: - -mat_mffd_check_positivity
780: .keywords: default, matrix-free, create, matrix
782: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
783: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
784: MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
786: @*/
787: PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
788: {
792: MatCreate(comm,J);
793: MatSetSizes(*J,m,n,M,N);
794: MatSetType(*J,MATMFFD);
795: MatSetUp(*J);
796: return(0);
797: }
799: /*@
800: MatMFFDGetH - Gets the last value that was used as the differencing
801: parameter.
803: Not Collective
805: Input Parameters:
806: . mat - the matrix obtained with MatCreateSNESMF()
808: Output Paramter:
809: . h - the differencing step size
811: Level: advanced
813: .keywords: SNES, matrix-free, parameters
815: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
816: @*/
817: PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
818: {
819: MatMFFD ctx = (MatMFFD)mat->data;
821: PetscBool match;
826: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
827: if (!match) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
829: *h = ctx->currenth;
830: return(0);
831: }
833: /*@C
834: MatMFFDSetFunction - Sets the function used in applying the matrix free.
836: Logically Collective on Mat
838: Input Parameters:
839: + mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
840: . func - the function to use
841: - funcctx - optional function context passed to function
843: Calling Sequence of func:
844: $ func (void *funcctx, Vec x, Vec f)
846: + funcctx - user provided context
847: . x - input vector
848: - f - computed output function
850: Level: advanced
852: Notes:
853: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
854: matrix inside your compute Jacobian routine
856: If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
858: .keywords: SNES, matrix-free, function
860: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
861: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
862: @*/
863: PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
864: {
869: PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
870: return(0);
871: }
873: /*@C
874: MatMFFDSetFunctioni - Sets the function for a single component
876: Logically Collective on Mat
878: Input Parameters:
879: + mat - the matrix free matrix created via MatCreateSNESMF()
880: - funci - the function to use
882: Level: advanced
884: Notes:
885: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
886: matrix inside your compute Jacobian routine.
887: This function is necessary to compute the diagonal of the matrix.
888: funci must not contain any MPI call as it is called inside a loop on the local portion of the vector.
890: .keywords: SNES, matrix-free, function
892: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
894: @*/
895: PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
896: {
901: PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
902: return(0);
903: }
905: /*@C
906: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
908: Logically Collective on Mat
910: Input Parameters:
911: + mat - the matrix free matrix created via MatCreateSNESMF()
912: - func - the function to use
914: Level: advanced
916: Notes:
917: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
918: matrix inside your compute Jacobian routine.
919: This function is necessary to compute the diagonal of the matrix.
922: .keywords: SNES, matrix-free, function
924: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
925: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction(), MatGetDiagonal()
926: @*/
927: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
928: {
933: PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
934: return(0);
935: }
937: /*@
938: MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime
940: Logically Collective on Mat
942: Input Parameters:
943: + mat - the matrix free matrix created via MatCreateSNESMF()
944: - period - 1 for everytime, 2 for every second etc
946: Options Database Keys:
947: + -mat_mffd_period <period>
949: Level: advanced
952: .keywords: SNES, matrix-free, parameters
954: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
955: MatMFFDSetHHistory(), MatMFFDResetHHistory()
956: @*/
957: PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
958: {
964: PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
965: return(0);
966: }
968: /*@
969: MatMFFDSetFunctionError - Sets the error_rel for the approximation of
970: matrix-vector products using finite differences.
972: Logically Collective on Mat
974: Input Parameters:
975: + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
976: - error_rel - relative error (should be set to the square root of
977: the relative error in the function evaluations)
979: Options Database Keys:
980: + -mat_mffd_err <error_rel> - Sets error_rel
982: Level: advanced
984: Notes:
985: The default matrix-free matrix-vector product routine computes
986: .vb
987: F'(u)*a = [F(u+h*a) - F(u)]/h where
988: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
989: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
990: .ve
992: .keywords: SNES, matrix-free, parameters
994: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
995: MatMFFDSetHHistory(), MatMFFDResetHHistory()
996: @*/
997: PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
998: {
1004: PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1005: return(0);
1006: }
1008: /*@
1009: MatMFFDSetHHistory - Sets an array to collect a history of the
1010: differencing values (h) computed for the matrix-free product.
1012: Logically Collective on Mat
1014: Input Parameters:
1015: + J - the matrix-free matrix context
1016: . histroy - space to hold the history
1017: - nhistory - number of entries in history, if more entries are generated than
1018: nhistory, then the later ones are discarded
1020: Level: advanced
1022: Notes:
1023: Use MatMFFDResetHHistory() to reset the history counter and collect
1024: a new batch of differencing parameters, h.
1026: .keywords: SNES, matrix-free, h history, differencing history
1028: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1029: MatMFFDResetHHistory(), MatMFFDSetFunctionError()
1031: @*/
1032: PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1033: {
1034: MatMFFD ctx = (MatMFFD)J->data;
1036: PetscBool match;
1042: PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1043: if (!match) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1044: ctx->historyh = history;
1045: ctx->maxcurrenth = nhistory;
1046: ctx->currenth = 0.;
1047: return(0);
1048: }
1050: /*@
1051: MatMFFDResetHHistory - Resets the counter to zero to begin
1052: collecting a new set of differencing histories.
1054: Logically Collective on Mat
1056: Input Parameters:
1057: . J - the matrix-free matrix context
1059: Level: advanced
1061: Notes:
1062: Use MatMFFDSetHHistory() to create the original history counter.
1064: .keywords: SNES, matrix-free, h history, differencing history
1066: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1067: MatMFFDSetHHistory(), MatMFFDSetFunctionError()
1069: @*/
1070: PetscErrorCode MatMFFDResetHHistory(Mat J)
1071: {
1076: PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1077: return(0);
1078: }
1080: /*@
1081: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1082: Jacobian are computed
1084: Logically Collective on Mat
1086: Input Parameters:
1087: + J - the MatMFFD matrix
1088: . U - the vector
1089: - F - (optional) vector that contains F(u) if it has been already computed
1091: Notes: This is rarely used directly
1093: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
1094: point during the first MatMult() after each call to MatMFFDSetBase().
1096: Level: advanced
1098: @*/
1099: PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
1100: {
1107: PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1108: return(0);
1109: }
1111: /*@C
1112: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1113: it to satisfy some criteria
1115: Logically Collective on Mat
1117: Input Parameters:
1118: + J - the MatMFFD matrix
1119: . fun - the function that checks h
1120: - ctx - any context needed by the function
1122: Options Database Keys:
1123: . -mat_mffd_check_positivity
1125: Level: advanced
1127: Notes: For example, MatMFFDCheckPositivity() insures that all entries
1128: of U + h*a are non-negative
1130: The function you provide is called after the default h has been computed and allows you to
1131: modify it.
1133: .seealso: MatMFFDCheckPositivity()
1134: @*/
1135: PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1136: {
1141: PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1142: return(0);
1143: }
1145: /*@
1146: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1147: zero, decreases h until this is satisfied.
1149: Logically Collective on Vec
1151: Input Parameters:
1152: + U - base vector that is added to
1153: . a - vector that is added
1154: . h - scaling factor on a
1155: - dummy - context variable (unused)
1157: Options Database Keys:
1158: . -mat_mffd_check_positivity
1160: Level: advanced
1162: Notes: This is rarely used directly, rather it is passed as an argument to
1163: MatMFFDSetCheckh()
1165: .seealso: MatMFFDSetCheckh()
1166: @*/
1167: PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1168: {
1169: PetscReal val, minval;
1170: PetscScalar *u_vec, *a_vec;
1172: PetscInt i,n;
1173: MPI_Comm comm;
1179: PetscObjectGetComm((PetscObject)U,&comm);
1180: VecGetArray(U,&u_vec);
1181: VecGetArray(a,&a_vec);
1182: VecGetLocalSize(U,&n);
1183: minval = PetscAbsScalar(*h*1.01);
1184: for (i=0; i<n; i++) {
1185: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1186: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1187: if (val < minval) minval = val;
1188: }
1189: }
1190: VecRestoreArray(U,&u_vec);
1191: VecRestoreArray(a,&a_vec);
1192: MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1193: if (val <= PetscAbsScalar(*h)) {
1194: PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1195: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1196: else *h = -0.99*val;
1197: }
1198: return(0);
1199: }