Actual source code: mffd.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/matimpl.h>
3: #include <../src/mat/impls/mffd/mffdimpl.h> /*I "petscmat.h" I*/
5: PetscFList MatMFFDList = 0;
6: PetscBool MatMFFDRegisterAllCalled = PETSC_FALSE;
8: PetscClassId MATMFFD_CLASSID;
9: PetscLogEvent MATMFFD_Mult;
11: static PetscBool MatMFFDPackageInitialized = PETSC_FALSE;
14: /*@C
15: MatMFFDFinalizePackage - This function destroys everything in the MatMFFD package. It is
16: called from PetscFinalize().
18: Level: developer
20: .keywords: Petsc, destroy, package
21: .seealso: PetscFinalize()
22: @*/
23: PetscErrorCode MatMFFDFinalizePackage(void)
24: {
26: MatMFFDPackageInitialized = PETSC_FALSE;
27: MatMFFDRegisterAllCalled = PETSC_FALSE;
28: MatMFFDList = PETSC_NULL;
29: return(0);
30: }
34: /*@C
35: MatMFFDInitializePackage - This function initializes everything in the MatMFFD package. It is called
36: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to MatCreate_MFFD()
37: when using static libraries.
39: Input Parameter:
40: . path - The dynamic library path, or PETSC_NULL
42: Level: developer
44: .keywords: Vec, initialize, package
45: .seealso: PetscInitialize()
46: @*/
47: PetscErrorCode MatMFFDInitializePackage(const char path[])
48: {
49: char logList[256];
50: char *className;
51: PetscBool opt;
52: PetscErrorCode ierr;
55: if (MatMFFDPackageInitialized) return(0);
56: MatMFFDPackageInitialized = PETSC_TRUE;
57: /* Register Classes */
58: PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);
59: /* Register Constructors */
60: MatMFFDRegisterAll(path);
61: /* Register Events */
62: PetscLogEventRegister("MatMult MF", MATMFFD_CLASSID,&MATMFFD_Mult);
64: /* Process info exclusions */
65: PetscOptionsGetString(PETSC_NULL, "-info_exclude", logList, 256, &opt);
66: if (opt) {
67: PetscStrstr(logList, "matmffd", &className);
68: if (className) {
69: PetscInfoDeactivateClass(MATMFFD_CLASSID);
70: }
71: }
72: /* Process summary exclusions */
73: PetscOptionsGetString(PETSC_NULL, "-log_summary_exclude", logList, 256, &opt);
74: if (opt) {
75: PetscStrstr(logList, "matmffd", &className);
76: if (className) {
77: PetscLogEventDeactivateClass(MATMFFD_CLASSID);
78: }
79: }
80: PetscRegisterFinalize(MatMFFDFinalizePackage);
81: return(0);
82: }
86: /*@C
87: MatMFFDSetType - Sets the method that is used to compute the
88: differencing parameter for finite differene matrix-free formulations.
90: Input Parameters:
91: + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
92: or MatSetType(mat,MATMFFD);
93: - ftype - the type requested, either MATMFFD_WP or MATMFFD_DS
95: Level: advanced
97: Notes:
98: For example, such routines can compute h for use in
99: Jacobian-vector products of the form
101: F(x+ha) - F(x)
102: F'(u)a ~= ----------------
103: h
105: .seealso: MatCreateSNESMF(), MatMFFDRegisterDynamic(), MatMFFDSetFunction()
106: @*/
107: PetscErrorCode MatMFFDSetType(Mat mat,const MatMFFDType ftype)
108: {
109: PetscErrorCode ierr,(*r)(MatMFFD);
110: MatMFFD ctx = (MatMFFD)mat->data;
111: PetscBool match;
112:
117: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
118: if (!match) return(0);
120: /* already set, so just return */
121: PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
122: if (match) return(0);
124: /* destroy the old one if it exists */
125: if (ctx->ops->destroy) {
126: (*ctx->ops->destroy)(ctx);
127: }
129: PetscFListFind(MatMFFDList,((PetscObject)ctx)->comm,ftype,PETSC_TRUE,(void (**)(void)) &r);
130: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
131: (*r)(ctx);
132: PetscObjectChangeTypeName((PetscObject)ctx,ftype);
133: return(0);
134: }
136: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
137: EXTERN_C_BEGIN
140: PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
141: {
142: MatMFFD ctx = (MatMFFD)mat->data;
145: ctx->funcisetbase = func;
146: return(0);
147: }
148: EXTERN_C_END
150: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
151: EXTERN_C_BEGIN
154: PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
155: {
156: MatMFFD ctx = (MatMFFD)mat->data;
159: ctx->funci = funci;
160: return(0);
161: }
162: EXTERN_C_END
164: EXTERN_C_BEGIN
167: PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
168: {
169: MatMFFD ctx = (MatMFFD)J->data;
172: ctx->ncurrenth = 0;
173: return(0);
174: }
175: EXTERN_C_END
179: PetscErrorCode MatMFFDRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(MatMFFD))
180: {
182: char fullname[PETSC_MAX_PATH_LEN];
185: PetscFListConcat(path,name,fullname);
186: PetscFListAdd(&MatMFFDList,sname,fullname,(void (*)(void))function);
187: return(0);
188: }
193: /*@C
194: MatMFFDRegisterDestroy - Frees the list of MatMFFD methods that were
195: registered by MatMFFDRegisterDynamic).
197: Not Collective
199: Level: developer
201: .keywords: MatMFFD, register, destroy
203: .seealso: MatMFFDRegisterDynamic), MatMFFDRegisterAll()
204: @*/
205: PetscErrorCode MatMFFDRegisterDestroy(void)
206: {
210: PetscFListDestroy(&MatMFFDList);
211: MatMFFDRegisterAllCalled = PETSC_FALSE;
212: return(0);
213: }
215: EXTERN_C_BEGIN
218: PetscErrorCode MatMFFDAddNullSpace_MFFD(Mat J,MatNullSpace nullsp)
219: {
221: MatMFFD ctx = (MatMFFD)J->data;
224: PetscObjectReference((PetscObject)nullsp);
225: if (ctx->sp) { MatNullSpaceDestroy(&ctx->sp); }
226: ctx->sp = nullsp;
227: return(0);
228: }
229: EXTERN_C_END
231: /* ----------------------------------------------------------------------------------------*/
234: PetscErrorCode MatDestroy_MFFD(Mat mat)
235: {
237: MatMFFD ctx = (MatMFFD)mat->data;
240: VecDestroy(&ctx->w);
241: VecDestroy(&ctx->drscale);
242: VecDestroy(&ctx->dlscale);
243: VecDestroy(&ctx->dshift);
244: if (ctx->current_f_allocated) {
245: VecDestroy(&ctx->current_f);
246: }
247: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
248: MatNullSpaceDestroy(&ctx->sp);
249: PetscHeaderDestroy(&ctx);
250: mat->data = 0;
252: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C","",PETSC_NULL);
253: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C","",PETSC_NULL);
254: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C","",PETSC_NULL);
255: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C","",PETSC_NULL);
256: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C","",PETSC_NULL);
257: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C","",PETSC_NULL);
258: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C","",PETSC_NULL);
259: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C","",PETSC_NULL);
260: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDAddNullSpace_C","",PETSC_NULL);
262: return(0);
263: }
267: /*
268: MatMFFDView_MFFD - Views matrix-free parameters.
270: */
271: PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
272: {
274: MatMFFD ctx = (MatMFFD)J->data;
275: PetscBool iascii, viewbase, viewfunction;
276: const char* prefix;
279: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
280: if (iascii) {
281: PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
282: PetscViewerASCIIPushTab(viewer);
283: PetscViewerASCIIPrintf(viewer,"err=%G (relative error in function evaluation)\n",ctx->error_rel);
284: if (!((PetscObject)ctx)->type_name) {
285: PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
286: } else {
287: PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
288: }
289: if (ctx->ops->view) {
290: (*ctx->ops->view)(ctx,viewer);
291: }
292: PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);
293:
294: PetscOptionsHasName(prefix, "-mat_mffd_view_base", &viewbase);
295: if(viewbase) {
296: PetscViewerASCIIPrintf(viewer, "Base:\n");
297: VecView(ctx->current_u, viewer);
298: }
299: PetscOptionsHasName(prefix, "-mat_mffd_view_function", &viewfunction);
300: if(viewfunction) {
301: PetscViewerASCIIPrintf(viewer, "Function:\n");
302: VecView(ctx->current_f, viewer);
303: }
304: PetscViewerASCIIPopTab(viewer);
305: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for matrix-free matrix",((PetscObject)viewer)->type_name);
306: return(0);
307: }
311: /*
312: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
313: allows the user to indicate the beginning of a new linear solve by calling
314: MatAssemblyXXX() on the matrix free matrix. This then allows the
315: MatCreateMFFD_WP() to properly compute ||U|| only the first time
316: in the linear solver rather than every time.
317: */
318: PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
319: {
321: MatMFFD j = (MatMFFD)J->data;
324: MatMFFDResetHHistory(J);
325: j->vshift = 0.0;
326: j->vscale = 1.0;
327: return(0);
328: }
332: /*
333: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
335: y ~= (F(u + ha) - F(u))/h,
336: where F = nonlinear function, as set by SNESSetFunction()
337: u = current iterate
338: h = difference interval
339: */
340: PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
341: {
342: MatMFFD ctx = (MatMFFD)mat->data;
343: PetscScalar h;
344: Vec w,U,F;
346: PetscBool zeroa;
349: if (!ctx->current_u) SETERRQ(((PetscObject)mat)->comm,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");
350: /* We log matrix-free matrix-vector products separately, so that we can
351: separate the performance monitoring from the cases that use conventional
352: storage. We may eventually modify event logging to associate events
353: with particular objects, hence alleviating the more general problem. */
354: PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);
356: w = ctx->w;
357: U = ctx->current_u;
358: F = ctx->current_f;
359: /*
360: Compute differencing parameter
361: */
362: if (!ctx->ops->compute) {
363: MatMFFDSetType(mat,MATMFFD_WP);
364: MatSetFromOptions(mat);
365: }
366: (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
367: if (zeroa) {
368: VecSet(y,0.0);
369: return(0);
370: }
372: if (PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
373: if (ctx->checkh) {
374: (*ctx->checkh)(ctx->checkhctx,U,a,&h);
375: }
377: /* keep a record of the current differencing parameter h */
378: ctx->currenth = h;
379: #if defined(PETSC_USE_COMPLEX)
380: PetscInfo2(mat,"Current differencing parameter: %G + %G i\n",PetscRealPart(h),PetscImaginaryPart(h));
381: #else
382: PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
383: #endif
384: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
385: ctx->historyh[ctx->ncurrenth] = h;
386: }
387: ctx->ncurrenth++;
389: /* w = u + ha */
390: if (ctx->drscale) {
391: VecPointwiseMult(ctx->drscale,a,U);
392: VecAYPX(U,h,w);
393: } else {
394: VecWAXPY(w,h,a,U);
395: }
397: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
398: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
399: (*ctx->func)(ctx->funcctx,U,F);
400: }
401: (*ctx->func)(ctx->funcctx,w,y);
403: VecAXPY(y,-1.0,F);
404: VecScale(y,1.0/h);
406: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
407: VecAXPBY(y,ctx->vshift,ctx->vscale,a);
408: }
409: if (ctx->dlscale) {
410: VecPointwiseMult(y,ctx->dlscale,y);
411: }
412: if (ctx->dshift) {
413: VecPointwiseMult(ctx->dshift,a,U);
414: VecAXPY(y,1.0,U);
415: }
417: if (ctx->sp) {MatNullSpaceRemove(ctx->sp,y,PETSC_NULL);}
419: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
420: return(0);
421: }
425: /*
426: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
428: y ~= (F(u + ha) - F(u))/h,
429: where F = nonlinear function, as set by SNESSetFunction()
430: u = current iterate
431: h = difference interval
432: */
433: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
434: {
435: MatMFFD ctx = (MatMFFD)mat->data;
436: PetscScalar h,*aa,*ww,v;
437: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
438: Vec w,U;
440: PetscInt i,rstart,rend;
443: if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Requires calling MatMFFDSetFunctioni() first");
445: w = ctx->w;
446: U = ctx->current_u;
447: (*ctx->func)(ctx->funcctx,U,a);
448: (*ctx->funcisetbase)(ctx->funcctx,U);
449: VecCopy(U,w);
451: VecGetOwnershipRange(a,&rstart,&rend);
452: VecGetArray(a,&aa);
453: for (i=rstart; i<rend; i++) {
454: VecGetArray(w,&ww);
455: h = ww[i-rstart];
456: if (h == 0.0) h = 1.0;
457: #if !defined(PETSC_USE_COMPLEX)
458: if (h < umin && h >= 0.0) h = umin;
459: else if (h < 0.0 && h > -umin) h = -umin;
460: #else
461: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
462: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
463: #endif
464: h *= epsilon;
465:
466: ww[i-rstart] += h;
467: VecRestoreArray(w,&ww);
468: (*ctx->funci)(ctx->funcctx,i,w,&v);
469: aa[i-rstart] = (v - aa[i-rstart])/h;
471: /* possibly shift and scale result */
472: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
473: aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
474: }
476: VecGetArray(w,&ww);
477: ww[i-rstart] -= h;
478: VecRestoreArray(w,&ww);
479: }
480: VecRestoreArray(a,&aa);
481: return(0);
482: }
486: PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
487: {
488: MatMFFD aij = (MatMFFD)mat->data;
492: if (ll && !aij->dlscale) {
493: VecDuplicate(ll,&aij->dlscale);
494: }
495: if (rr && !aij->drscale) {
496: VecDuplicate(rr,&aij->drscale);
497: }
498: if (ll) {
499: VecCopy(ll,aij->dlscale);
500: }
501: if (rr) {
502: VecCopy(rr,aij->drscale);
503: }
504: return(0);
505: }
509: PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
510: {
511: MatMFFD aij = (MatMFFD)mat->data;
515: if (mode == INSERT_VALUES) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
516: if (!aij->dshift) {
517: VecDuplicate(ll,&aij->dshift);
518: }
519: VecAXPY(aij->dshift,1.0,ll);
520: return(0);
521: }
525: PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
526: {
527: MatMFFD shell = (MatMFFD)Y->data;
529: shell->vshift += a;
530: return(0);
531: }
535: PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
536: {
537: MatMFFD shell = (MatMFFD)Y->data;
539: shell->vscale *= a;
540: return(0);
541: }
543: EXTERN_C_BEGIN
546: PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
547: {
549: MatMFFD ctx = (MatMFFD)J->data;
552: MatMFFDResetHHistory(J);
553: ctx->current_u = U;
554: if (F) {
555: if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
556: ctx->current_f = F;
557: ctx->current_f_allocated = PETSC_FALSE;
558: } else if (!ctx->current_f_allocated) {
559: VecDuplicate(ctx->current_u, &ctx->current_f);
560: ctx->current_f_allocated = PETSC_TRUE;
561: }
562: if (!ctx->w) {
563: VecDuplicate(ctx->current_u, &ctx->w);
564: }
565: J->assembled = PETSC_TRUE;
566: return(0);
567: }
568: EXTERN_C_END
569: typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
570: EXTERN_C_BEGIN
573: PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void*ectx)
574: {
575: MatMFFD ctx = (MatMFFD)J->data;
578: ctx->checkh = fun;
579: ctx->checkhctx = ectx;
580: return(0);
581: }
582: EXTERN_C_END
586: /*@C
587: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
588: MatMFFD options in the database.
590: Collective on Mat
592: Input Parameter:
593: + A - the Mat context
594: - prefix - the prefix to prepend to all option names
596: Notes:
597: A hyphen (-) must NOT be given at the beginning of the prefix name.
598: The first character of all runtime options is AUTOMATICALLY the hyphen.
600: Level: advanced
602: .keywords: SNES, matrix-free, parameters
604: .seealso: MatSetFromOptions(), MatCreateSNESMF()
605: @*/
606: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
608: {
609: MatMFFD mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)PETSC_NULL;
614: PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
615: return(0);
616: }
620: PetscErrorCode MatSetFromOptions_MFFD(Mat mat)
621: {
622: MatMFFD mfctx = (MatMFFD)mat->data;
624: PetscBool flg;
625: char ftype[256];
630: PetscObjectOptionsBegin((PetscObject)mfctx);
631: PetscOptionsList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
632: if (flg) {
633: MatMFFDSetType(mat,ftype);
634: }
636: PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
637: PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
639: flg = PETSC_FALSE;
640: PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,PETSC_NULL);
641: if (flg) {
642: MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
643: }
644: if (mfctx->ops->setfromoptions) {
645: (*mfctx->ops->setfromoptions)(mfctx);
646: }
647: PetscOptionsEnd();
648: return(0);
649: }
651: EXTERN_C_BEGIN
654: PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
655: {
656: MatMFFD ctx = (MatMFFD)mat->data;
660: ctx->recomputeperiod = period;
661: return(0);
662: }
663: EXTERN_C_END
665: EXTERN_C_BEGIN
668: PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
669: {
670: MatMFFD ctx = (MatMFFD)mat->data;
673: ctx->func = func;
674: ctx->funcctx = funcctx;
675: return(0);
676: }
677: EXTERN_C_END
679: EXTERN_C_BEGIN
682: PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
683: {
684: MatMFFD ctx = (MatMFFD)mat->data;
688: if (error != PETSC_DEFAULT) ctx->error_rel = error;
689: return(0);
690: }
691: EXTERN_C_END
693: /*MC
694: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
696: Level: advanced
698: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction()
699: M*/
700: EXTERN_C_BEGIN
703: PetscErrorCode MatCreate_MFFD(Mat A)
704: {
705: MatMFFD mfctx;
706: PetscErrorCode ierr;
709: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
710: MatMFFDInitializePackage(PETSC_NULL);
711: #endif
713: PetscHeaderCreate(mfctx,_p_MatMFFD,struct _MFOps,MATMFFD_CLASSID,0,"MatMFFD","Matrix-free Finite Differencing","Mat",((PetscObject)A)->comm,MatDestroy_MFFD,MatView_MFFD);
714: mfctx->sp = 0;
715: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
716: mfctx->recomputeperiod = 1;
717: mfctx->count = 0;
718: mfctx->currenth = 0.0;
719: mfctx->historyh = PETSC_NULL;
720: mfctx->ncurrenth = 0;
721: mfctx->maxcurrenth = 0;
722: ((PetscObject)mfctx)->type_name = 0;
724: mfctx->vshift = 0.0;
725: mfctx->vscale = 1.0;
727: /*
728: Create the empty data structure to contain compute-h routines.
729: These will be filled in below from the command line options or
730: a later call with MatMFFDSetType() or if that is not called
731: then it will default in the first use of MatMult_MFFD()
732: */
733: mfctx->ops->compute = 0;
734: mfctx->ops->destroy = 0;
735: mfctx->ops->view = 0;
736: mfctx->ops->setfromoptions = 0;
737: mfctx->hctx = 0;
739: mfctx->func = 0;
740: mfctx->funcctx = 0;
741: mfctx->w = PETSC_NULL;
743: A->data = mfctx;
745: A->ops->mult = MatMult_MFFD;
746: A->ops->destroy = MatDestroy_MFFD;
747: A->ops->view = MatView_MFFD;
748: A->ops->assemblyend = MatAssemblyEnd_MFFD;
749: A->ops->getdiagonal = MatGetDiagonal_MFFD;
750: A->ops->scale = MatScale_MFFD;
751: A->ops->shift = MatShift_MFFD;
752: A->ops->diagonalscale = MatDiagonalScale_MFFD;
753: A->ops->diagonalset = MatDiagonalSet_MFFD;
754: A->ops->setfromoptions = MatSetFromOptions_MFFD;
755: A->assembled = PETSC_TRUE;
757: PetscLayoutSetUp(A->rmap);
758: PetscLayoutSetUp(A->cmap);
760: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMFFDSetBase_C","MatMFFDSetBase_MFFD",MatMFFDSetBase_MFFD);
761: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMFFDSetFunctioniBase_C","MatMFFDSetFunctioniBase_MFFD",MatMFFDSetFunctioniBase_MFFD);
762: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMFFDSetFunctioni_C","MatMFFDSetFunctioni_MFFD",MatMFFDSetFunctioni_MFFD);
763: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMFFDSetFunction_C","MatMFFDSetFunction_MFFD",MatMFFDSetFunction_MFFD);
764: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMFFDSetCheckh_C","MatMFFDSetCheckh_MFFD",MatMFFDSetCheckh_MFFD);
765: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMFFDSetPeriod_C","MatMFFDSetPeriod_MFFD",MatMFFDSetPeriod_MFFD);
766: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMFFDSetFunctionError_C","MatMFFDSetFunctionError_MFFD",MatMFFDSetFunctionError_MFFD);
767: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMFFDResetHHistory_C","MatMFFDResetHHistory_MFFD",MatMFFDResetHHistory_MFFD);
768: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMFFDAddNullSpace_C","MatMFFDAddNullSpace_MFFD",MatMFFDAddNullSpace_MFFD);
769: mfctx->mat = A;
770: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
771: return(0);
772: }
773: EXTERN_C_END
777: /*@
778: MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
780: Collective on Vec
782: Input Parameters:
783: + comm - MPI communicator
784: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
785: This value should be the same as the local size used in creating the
786: y vector for the matrix-vector product y = Ax.
787: . n - This value should be the same as the local size used in creating the
788: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
789: calculated if N is given) For square matrices n is almost always m.
790: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
791: - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
794: Output Parameter:
795: . J - the matrix-free matrix
797: Options Database Keys: call MatSetFromOptions() to trigger these
798: + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
799: - -mat_mffd_err - square root of estimated relative error in function evaluation
800: - -mat_mffd_period - how often h is recomputed, defaults to 1, everytime
803: Level: advanced
805: Notes:
806: The matrix-free matrix context merely contains the function pointers
807: and work space for performing finite difference approximations of
808: Jacobian-vector products, F'(u)*a,
810: The default code uses the following approach to compute h
812: .vb
813: F'(u)*a = [F(u+h*a) - F(u)]/h where
814: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
815: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
816: where
817: error_rel = square root of relative error in function evaluation
818: umin = minimum iterate parameter
819: .ve
821: The user can set the error_rel via MatMFFDSetFunctionError() and
822: umin via MatMFFDDSSetUmin(); see the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A> for details.
824: The user should call MatDestroy() when finished with the matrix-free
825: matrix context.
827: Options Database Keys:
828: + -mat_mffd_err <error_rel> - Sets error_rel
829: . -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
830: - -mat_mffd_check_positivity
832: .keywords: default, matrix-free, create, matrix
834: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
835: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
836: MatMFFDGetH(), MatMFFDRegisterDynamic), MatMFFDComputeJacobian()
837:
838: @*/
839: PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
840: {
844: MatCreate(comm,J);
845: MatSetSizes(*J,m,n,M,N);
846: MatSetType(*J,MATMFFD);
847: MatSetUp(*J);
848: return(0);
849: }
854: /*@
855: MatMFFDGetH - Gets the last value that was used as the differencing
856: parameter.
858: Not Collective
860: Input Parameters:
861: . mat - the matrix obtained with MatCreateSNESMF()
863: Output Paramter:
864: . h - the differencing step size
866: Level: advanced
868: .keywords: SNES, matrix-free, parameters
870: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
871: @*/
872: PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
873: {
874: MatMFFD ctx = (MatMFFD)mat->data;
876: PetscBool match;
879: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
880: if (!match) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
882: *h = ctx->currenth;
883: return(0);
884: }
888: /*@C
889: MatMFFDSetFunction - Sets the function used in applying the matrix free.
891: Logically Collective on Mat
893: Input Parameters:
894: + mat - the matrix free matrix created via MatCreateSNESMF()
895: . func - the function to use
896: - funcctx - optional function context passed to function
898: Level: advanced
900: Notes:
901: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
902: matrix inside your compute Jacobian routine
904: If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
906: .keywords: SNES, matrix-free, function
908: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
909: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
910: @*/
911: PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
912: {
916: PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
917: return(0);
918: }
922: /*@C
923: MatMFFDSetFunctioni - Sets the function for a single component
925: Logically Collective on Mat
927: Input Parameters:
928: + mat - the matrix free matrix created via MatCreateSNESMF()
929: - funci - the function to use
931: Level: advanced
933: Notes:
934: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
935: matrix inside your compute Jacobian routine
938: .keywords: SNES, matrix-free, function
940: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
942: @*/
943: PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
944: {
949: PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
950: return(0);
951: }
956: /*@C
957: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
959: Logically Collective on Mat
961: Input Parameters:
962: + mat - the matrix free matrix created via MatCreateSNESMF()
963: - func - the function to use
965: Level: advanced
967: Notes:
968: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
969: matrix inside your compute Jacobian routine
972: .keywords: SNES, matrix-free, function
974: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
975: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
976: @*/
977: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
978: {
983: PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
984: return(0);
985: }
989: /*@
990: MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime
992: Logically Collective on Mat
994: Input Parameters:
995: + mat - the matrix free matrix created via MatCreateSNESMF()
996: - period - 1 for everytime, 2 for every second etc
998: Options Database Keys:
999: + -mat_mffd_period <period>
1001: Level: advanced
1004: .keywords: SNES, matrix-free, parameters
1006: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
1007: MatMFFDSetHHistory(), MatMFFDResetHHistory()
1008: @*/
1009: PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
1010: {
1014: PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
1015: return(0);
1016: }
1020: /*@
1021: MatMFFDSetFunctionError - Sets the error_rel for the approximation of
1022: matrix-vector products using finite differences.
1024: Logically Collective on Mat
1026: Input Parameters:
1027: + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
1028: - error_rel - relative error (should be set to the square root of
1029: the relative error in the function evaluations)
1031: Options Database Keys:
1032: + -mat_mffd_err <error_rel> - Sets error_rel
1034: Level: advanced
1036: Notes:
1037: The default matrix-free matrix-vector product routine computes
1038: .vb
1039: F'(u)*a = [F(u+h*a) - F(u)]/h where
1040: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
1041: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
1042: .ve
1044: .keywords: SNES, matrix-free, parameters
1046: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
1047: MatMFFDSetHHistory(), MatMFFDResetHHistory()
1048: @*/
1049: PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
1050: {
1054: PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1055: return(0);
1056: }
1060: /*@
1061: MatMFFDAddNullSpace - Provides a null space that an operator is
1062: supposed to have. Since roundoff will create a small component in
1063: the null space, if you know the null space you may have it
1064: automatically removed.
1066: Logically Collective on Mat
1068: Input Parameters:
1069: + J - the matrix-free matrix context
1070: - nullsp - object created with MatNullSpaceCreate()
1072: Level: advanced
1074: .keywords: SNES, matrix-free, null space
1076: .seealso: MatNullSpaceCreate(), MatMFFDGetH(), MatCreateSNESMF(), MatCreateMFFD(), MATMFFD
1077: MatMFFDSetHHistory(), MatMFFDResetHHistory()
1078: @*/
1079: PetscErrorCode MatMFFDAddNullSpace(Mat J,MatNullSpace nullsp)
1080: {
1084: PetscTryMethod(J,"MatMFFDAddNullSpace_C",(Mat,MatNullSpace),(J,nullsp));
1085: return(0);
1086: }
1090: /*@
1091: MatMFFDSetHHistory - Sets an array to collect a history of the
1092: differencing values (h) computed for the matrix-free product.
1094: Logically Collective on Mat
1096: Input Parameters:
1097: + J - the matrix-free matrix context
1098: . histroy - space to hold the history
1099: - nhistory - number of entries in history, if more entries are generated than
1100: nhistory, then the later ones are discarded
1102: Level: advanced
1104: Notes:
1105: Use MatMFFDResetHHistory() to reset the history counter and collect
1106: a new batch of differencing parameters, h.
1108: .keywords: SNES, matrix-free, h history, differencing history
1110: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1111: MatMFFDResetHHistory(), MatMFFDSetFunctionError()
1113: @*/
1114: PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1115: {
1116: MatMFFD ctx = (MatMFFD)J->data;
1118: PetscBool match;
1121: PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1122: if (!match) SETERRQ(((PetscObject)J)->comm,PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1123: ctx->historyh = history;
1124: ctx->maxcurrenth = nhistory;
1125: ctx->currenth = 0.;
1126: return(0);
1127: }
1132: /*@
1133: MatMFFDResetHHistory - Resets the counter to zero to begin
1134: collecting a new set of differencing histories.
1136: Logically Collective on Mat
1138: Input Parameters:
1139: . J - the matrix-free matrix context
1141: Level: advanced
1143: Notes:
1144: Use MatMFFDSetHHistory() to create the original history counter.
1146: .keywords: SNES, matrix-free, h history, differencing history
1148: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1149: MatMFFDSetHHistory(), MatMFFDSetFunctionError()
1151: @*/
1152: PetscErrorCode MatMFFDResetHHistory(Mat J)
1153: {
1157: PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1158: return(0);
1159: }
1164: /*@
1165: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1166: Jacobian are computed
1168: Logically Collective on Mat
1170: Input Parameters:
1171: + J - the MatMFFD matrix
1172: . U - the vector
1173: - F - (optional) vector that contains F(u) if it has been already computed
1175: Notes: This is rarely used directly
1177: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
1178: point during the first MatMult() after each call to MatMFFDSetBase().
1180: Level: advanced
1182: @*/
1183: PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
1184: {
1191: PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1192: return(0);
1193: }
1197: /*@C
1198: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1199: it to satisfy some criteria
1201: Logically Collective on Mat
1203: Input Parameters:
1204: + J - the MatMFFD matrix
1205: . fun - the function that checks h
1206: - ctx - any context needed by the function
1208: Options Database Keys:
1209: . -mat_mffd_check_positivity
1211: Level: advanced
1213: Notes: For example, MatMFFDSetCheckPositivity() insures that all entries
1214: of U + h*a are non-negative
1216: .seealso: MatMFFDSetCheckPositivity()
1217: @*/
1218: PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void* ctx)
1219: {
1224: PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1225: return(0);
1226: }
1230: /*@
1231: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1232: zero, decreases h until this is satisfied.
1234: Logically Collective on Vec
1236: Input Parameters:
1237: + U - base vector that is added to
1238: . a - vector that is added
1239: . h - scaling factor on a
1240: - dummy - context variable (unused)
1242: Options Database Keys:
1243: . -mat_mffd_check_positivity
1245: Level: advanced
1247: Notes: This is rarely used directly, rather it is passed as an argument to
1248: MatMFFDSetCheckh()
1250: .seealso: MatMFFDSetCheckh()
1251: @*/
1252: PetscErrorCode MatMFFDCheckPositivity(void* dummy,Vec U,Vec a,PetscScalar *h)
1253: {
1254: PetscReal val, minval;
1255: PetscScalar *u_vec, *a_vec;
1257: PetscInt i,n;
1258: MPI_Comm comm;
1261: PetscObjectGetComm((PetscObject)U,&comm);
1262: VecGetArray(U,&u_vec);
1263: VecGetArray(a,&a_vec);
1264: VecGetLocalSize(U,&n);
1265: minval = PetscAbsScalar(*h*1.01);
1266: for(i=0;i<n;i++) {
1267: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1268: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1269: if (val < minval) minval = val;
1270: }
1271: }
1272: VecRestoreArray(U,&u_vec);
1273: VecRestoreArray(a,&a_vec);
1274: MPI_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1275: if (val <= PetscAbsScalar(*h)) {
1276: PetscInfo2(U,"Scaling back h from %G to %G\n",PetscRealPart(*h),.99*val);
1277: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1278: else *h = -0.99*val;
1279: }
1280: return(0);
1281: }