Actual source code: mffd.c
petsc-3.5.4 2015-05-23
2: #include <petsc-private/matimpl.h>
3: #include <../src/mat/impls/mffd/mffdimpl.h> /*I "petscmat.h" I*/
5: PetscFunctionList 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: {
28: PetscFunctionListDestroy(&MatMFFDList);
29: MatMFFDPackageInitialized = PETSC_FALSE;
30: MatMFFDRegisterAllCalled = PETSC_FALSE;
31: return(0);
32: }
36: /*@C
37: MatMFFDInitializePackage - This function initializes everything in the MatMFFD package. It is called
38: from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to MatCreate_MFFD()
39: when using static libraries.
41: Level: developer
43: .keywords: Vec, initialize, package
44: .seealso: PetscInitialize()
45: @*/
46: PetscErrorCode MatMFFDInitializePackage(void)
47: {
48: char logList[256];
49: char *className;
50: PetscBool opt;
54: if (MatMFFDPackageInitialized) return(0);
55: MatMFFDPackageInitialized = PETSC_TRUE;
56: /* Register Classes */
57: PetscClassIdRegister("MatMFFD",&MATMFFD_CLASSID);
58: /* Register Constructors */
59: MatMFFDRegisterAll();
60: /* Register Events */
61: PetscLogEventRegister("MatMult MF", MATMFFD_CLASSID,&MATMFFD_Mult);
63: /* Process info exclusions */
64: PetscOptionsGetString(NULL, "-info_exclude", logList, 256, &opt);
65: if (opt) {
66: PetscStrstr(logList, "matmffd", &className);
67: if (className) {
68: PetscInfoDeactivateClass(MATMFFD_CLASSID);
69: }
70: }
71: /* Process summary exclusions */
72: PetscOptionsGetString(NULL, "-log_summary_exclude", logList, 256, &opt);
73: if (opt) {
74: PetscStrstr(logList, "matmffd", &className);
75: if (className) {
76: PetscLogEventDeactivateClass(MATMFFD_CLASSID);
77: }
78: }
79: PetscRegisterFinalize(MatMFFDFinalizePackage);
80: return(0);
81: }
85: /*@C
86: MatMFFDSetType - Sets the method that is used to compute the
87: differencing parameter for finite differene matrix-free formulations.
89: Input Parameters:
90: + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
91: or MatSetType(mat,MATMFFD);
92: - ftype - the type requested, either MATMFFD_WP or MATMFFD_DS
94: Level: advanced
96: Notes:
97: For example, such routines can compute h for use in
98: Jacobian-vector products of the form
100: F(x+ha) - F(x)
101: F'(u)a ~= ----------------
102: h
104: .seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction()
105: @*/
106: PetscErrorCode MatMFFDSetType(Mat mat,MatMFFDType ftype)
107: {
108: PetscErrorCode ierr,(*r)(MatMFFD);
109: MatMFFD ctx = (MatMFFD)mat->data;
110: PetscBool match;
116: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
117: if (!match) return(0);
119: /* already set, so just return */
120: PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);
121: if (match) return(0);
123: /* destroy the old one if it exists */
124: if (ctx->ops->destroy) {
125: (*ctx->ops->destroy)(ctx);
126: }
128: PetscFunctionListFind(MatMFFDList,ftype,&r);
129: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
130: (*r)(ctx);
131: PetscObjectChangeTypeName((PetscObject)ctx,ftype);
132: return(0);
133: }
135: typedef PetscErrorCode (*FCN1)(void*,Vec); /* force argument to next function to not be extern C*/
138: PetscErrorCode MatMFFDSetFunctioniBase_MFFD(Mat mat,FCN1 func)
139: {
140: MatMFFD ctx = (MatMFFD)mat->data;
143: ctx->funcisetbase = func;
144: return(0);
145: }
147: typedef PetscErrorCode (*FCN2)(void*,PetscInt,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
150: PetscErrorCode MatMFFDSetFunctioni_MFFD(Mat mat,FCN2 funci)
151: {
152: MatMFFD ctx = (MatMFFD)mat->data;
155: ctx->funci = funci;
156: return(0);
157: }
161: PetscErrorCode MatMFFDResetHHistory_MFFD(Mat J)
162: {
163: MatMFFD ctx = (MatMFFD)J->data;
166: ctx->ncurrenth = 0;
167: return(0);
168: }
172: /*@C
173: MatMFFDRegister - Adds a method to the MatMFFD registry.
175: Not Collective
177: Input Parameters:
178: + name_solver - name of a new user-defined compute-h module
179: - routine_create - routine to create method context
181: Level: developer
183: Notes:
184: MatMFFDRegister() may be called multiple times to add several user-defined solvers.
186: Sample usage:
187: .vb
188: MatMFFDRegister("my_h",MyHCreate);
189: .ve
191: Then, your solver can be chosen with the procedural interface via
192: $ MatMFFDSetType(mfctx,"my_h")
193: or at runtime via the option
194: $ -mat_mffd_type my_h
196: .keywords: MatMFFD, register
198: .seealso: MatMFFDRegisterAll(), MatMFFDRegisterDestroy()
199: @*/
200: PetscErrorCode MatMFFDRegister(const char sname[],PetscErrorCode (*function)(MatMFFD))
201: {
205: PetscFunctionListAdd(&MatMFFDList,sname,function);
206: return(0);
207: }
211: PetscErrorCode MatMFFDAddNullSpace_MFFD(Mat J,MatNullSpace nullsp)
212: {
214: MatMFFD ctx = (MatMFFD)J->data;
217: PetscObjectReference((PetscObject)nullsp);
218: if (ctx->sp) { MatNullSpaceDestroy(&ctx->sp); }
219: ctx->sp = nullsp;
220: return(0);
221: }
223: /* ----------------------------------------------------------------------------------------*/
226: PetscErrorCode MatDestroy_MFFD(Mat mat)
227: {
229: MatMFFD ctx = (MatMFFD)mat->data;
232: VecDestroy(&ctx->w);
233: VecDestroy(&ctx->drscale);
234: VecDestroy(&ctx->dlscale);
235: VecDestroy(&ctx->dshift);
236: if (ctx->current_f_allocated) {
237: VecDestroy(&ctx->current_f);
238: }
239: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
240: MatNullSpaceDestroy(&ctx->sp);
241: PetscHeaderDestroy(&ctx);
242: mat->data = 0;
244: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
245: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
246: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
247: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
248: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
249: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
250: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
251: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
252: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDAddNullSpace_C",NULL);
253: return(0);
254: }
258: /*
259: MatMFFDView_MFFD - Views matrix-free parameters.
261: */
262: PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
263: {
265: MatMFFD ctx = (MatMFFD)J->data;
266: PetscBool iascii, viewbase, viewfunction;
267: const char *prefix;
270: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
271: if (iascii) {
272: PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
273: PetscViewerASCIIPushTab(viewer);
274: PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
275: if (!((PetscObject)ctx)->type_name) {
276: PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
277: } else {
278: PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
279: }
280: if (ctx->ops->view) {
281: (*ctx->ops->view)(ctx,viewer);
282: }
283: PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);
285: PetscOptionsHasName(prefix, "-mat_mffd_view_base", &viewbase);
286: if (viewbase) {
287: PetscViewerASCIIPrintf(viewer, "Base:\n");
288: VecView(ctx->current_u, viewer);
289: }
290: PetscOptionsHasName(prefix, "-mat_mffd_view_function", &viewfunction);
291: if (viewfunction) {
292: PetscViewerASCIIPrintf(viewer, "Function:\n");
293: VecView(ctx->current_f, viewer);
294: }
295: PetscViewerASCIIPopTab(viewer);
296: }
297: return(0);
298: }
302: /*
303: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
304: allows the user to indicate the beginning of a new linear solve by calling
305: MatAssemblyXXX() on the matrix free matrix. This then allows the
306: MatCreateMFFD_WP() to properly compute ||U|| only the first time
307: in the linear solver rather than every time.
309: This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library.
310: */
311: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
312: {
314: MatMFFD j = (MatMFFD)J->data;
317: MatMFFDResetHHistory(J);
318: j->vshift = 0.0;
319: j->vscale = 1.0;
320: return(0);
321: }
325: /*
326: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
328: y ~= (F(u + ha) - F(u))/h,
329: where F = nonlinear function, as set by SNESSetFunction()
330: u = current iterate
331: h = difference interval
332: */
333: PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
334: {
335: MatMFFD ctx = (MatMFFD)mat->data;
336: PetscScalar h;
337: Vec w,U,F;
339: PetscBool zeroa;
342: 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");
343: /* We log matrix-free matrix-vector products separately, so that we can
344: separate the performance monitoring from the cases that use conventional
345: storage. We may eventually modify event logging to associate events
346: with particular objects, hence alleviating the more general problem. */
347: PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);
349: w = ctx->w;
350: U = ctx->current_u;
351: F = ctx->current_f;
352: /*
353: Compute differencing parameter
354: */
355: if (!ctx->ops->compute) {
356: MatMFFDSetType(mat,MATMFFD_WP);
357: MatSetFromOptions(mat);
358: }
359: (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
360: if (zeroa) {
361: VecSet(y,0.0);
362: return(0);
363: }
365: if (PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
366: if (ctx->checkh) {
367: (*ctx->checkh)(ctx->checkhctx,U,a,&h);
368: }
370: /* keep a record of the current differencing parameter h */
371: ctx->currenth = h;
372: #if defined(PETSC_USE_COMPLEX)
373: PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
374: #else
375: PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
376: #endif
377: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
378: ctx->historyh[ctx->ncurrenth] = h;
379: }
380: ctx->ncurrenth++;
382: /* w = u + ha */
383: if (ctx->drscale) {
384: VecPointwiseMult(ctx->drscale,a,U);
385: VecAYPX(U,h,w);
386: } else {
387: VecWAXPY(w,h,a,U);
388: }
390: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
391: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
392: (*ctx->func)(ctx->funcctx,U,F);
393: }
394: (*ctx->func)(ctx->funcctx,w,y);
396: VecAXPY(y,-1.0,F);
397: VecScale(y,1.0/h);
399: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
400: VecAXPBY(y,ctx->vshift,ctx->vscale,a);
401: }
402: if (ctx->dlscale) {
403: VecPointwiseMult(y,ctx->dlscale,y);
404: }
405: if (ctx->dshift) {
406: VecPointwiseMult(ctx->dshift,a,U);
407: VecAXPY(y,1.0,U);
408: }
410: if (ctx->sp) {MatNullSpaceRemove(ctx->sp,y);}
412: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
413: return(0);
414: }
418: /*
419: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
421: y ~= (F(u + ha) - F(u))/h,
422: where F = nonlinear function, as set by SNESSetFunction()
423: u = current iterate
424: h = difference interval
425: */
426: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
427: {
428: MatMFFD ctx = (MatMFFD)mat->data;
429: PetscScalar h,*aa,*ww,v;
430: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
431: Vec w,U;
433: PetscInt i,rstart,rend;
436: if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Requires calling MatMFFDSetFunctioni() first");
438: w = ctx->w;
439: U = ctx->current_u;
440: (*ctx->func)(ctx->funcctx,U,a);
441: (*ctx->funcisetbase)(ctx->funcctx,U);
442: VecCopy(U,w);
444: VecGetOwnershipRange(a,&rstart,&rend);
445: VecGetArray(a,&aa);
446: for (i=rstart; i<rend; i++) {
447: VecGetArray(w,&ww);
448: h = ww[i-rstart];
449: if (h == 0.0) h = 1.0;
450: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
451: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
452: h *= epsilon;
454: ww[i-rstart] += h;
455: VecRestoreArray(w,&ww);
456: (*ctx->funci)(ctx->funcctx,i,w,&v);
457: aa[i-rstart] = (v - aa[i-rstart])/h;
459: /* possibly shift and scale result */
460: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
461: aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
462: }
464: VecGetArray(w,&ww);
465: ww[i-rstart] -= h;
466: VecRestoreArray(w,&ww);
467: }
468: VecRestoreArray(a,&aa);
469: return(0);
470: }
474: PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
475: {
476: MatMFFD aij = (MatMFFD)mat->data;
480: if (ll && !aij->dlscale) {
481: VecDuplicate(ll,&aij->dlscale);
482: }
483: if (rr && !aij->drscale) {
484: VecDuplicate(rr,&aij->drscale);
485: }
486: if (ll) {
487: VecCopy(ll,aij->dlscale);
488: }
489: if (rr) {
490: VecCopy(rr,aij->drscale);
491: }
492: return(0);
493: }
497: PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
498: {
499: MatMFFD aij = (MatMFFD)mat->data;
503: if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
504: if (!aij->dshift) {
505: VecDuplicate(ll,&aij->dshift);
506: }
507: VecAXPY(aij->dshift,1.0,ll);
508: return(0);
509: }
513: PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
514: {
515: MatMFFD shell = (MatMFFD)Y->data;
518: shell->vshift += a;
519: return(0);
520: }
524: PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
525: {
526: MatMFFD shell = (MatMFFD)Y->data;
529: shell->vscale *= a;
530: return(0);
531: }
535: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
536: {
538: MatMFFD ctx = (MatMFFD)J->data;
541: MatMFFDResetHHistory(J);
543: ctx->current_u = U;
544: if (F) {
545: if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
546: ctx->current_f = F;
547: ctx->current_f_allocated = PETSC_FALSE;
548: } else if (!ctx->current_f_allocated) {
549: MatGetVecs(J,NULL,&ctx->current_f);
551: ctx->current_f_allocated = PETSC_TRUE;
552: }
553: if (!ctx->w) {
554: VecDuplicate(ctx->current_u, &ctx->w);
555: }
556: J->assembled = PETSC_TRUE;
557: return(0);
558: }
560: typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
564: PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
565: {
566: MatMFFD ctx = (MatMFFD)J->data;
569: ctx->checkh = fun;
570: ctx->checkhctx = ectx;
571: return(0);
572: }
576: /*@C
577: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
578: MatMFFD options in the database.
580: Collective on Mat
582: Input Parameter:
583: + A - the Mat context
584: - prefix - the prefix to prepend to all option names
586: Notes:
587: A hyphen (-) must NOT be given at the beginning of the prefix name.
588: The first character of all runtime options is AUTOMATICALLY the hyphen.
590: Level: advanced
592: .keywords: SNES, matrix-free, parameters
594: .seealso: MatSetFromOptions(), MatCreateSNESMF()
595: @*/
596: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
598: {
599: MatMFFD mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;
605: PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
606: return(0);
607: }
611: PetscErrorCode MatSetFromOptions_MFFD(Mat mat)
612: {
613: MatMFFD mfctx = (MatMFFD)mat->data;
615: PetscBool flg;
616: char ftype[256];
621: PetscObjectOptionsBegin((PetscObject)mfctx);
622: PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
623: if (flg) {
624: MatMFFDSetType(mat,ftype);
625: }
627: PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
628: PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
630: flg = PETSC_FALSE;
631: PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
632: if (flg) {
633: MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
634: }
635: if (mfctx->ops->setfromoptions) {
636: (*mfctx->ops->setfromoptions)(mfctx);
637: }
638: PetscOptionsEnd();
639: return(0);
640: }
644: PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
645: {
646: MatMFFD ctx = (MatMFFD)mat->data;
650: ctx->recomputeperiod = period;
651: return(0);
652: }
656: PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
657: {
658: MatMFFD ctx = (MatMFFD)mat->data;
661: ctx->func = func;
662: ctx->funcctx = funcctx;
663: return(0);
664: }
668: PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
669: {
670: MatMFFD ctx = (MatMFFD)mat->data;
674: if (error != PETSC_DEFAULT) ctx->error_rel = error;
675: return(0);
676: }
678: /*MC
679: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
681: Level: advanced
683: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction()
684: M*/
687: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
688: {
689: MatMFFD mfctx;
693: MatMFFDInitializePackage();
695: PetscHeaderCreate(mfctx,_p_MatMFFD,struct _MFOps,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),MatDestroy_MFFD,MatView_MFFD);
697: mfctx->sp = 0;
698: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
699: mfctx->recomputeperiod = 1;
700: mfctx->count = 0;
701: mfctx->currenth = 0.0;
702: mfctx->historyh = NULL;
703: mfctx->ncurrenth = 0;
704: mfctx->maxcurrenth = 0;
705: ((PetscObject)mfctx)->type_name = 0;
707: mfctx->vshift = 0.0;
708: mfctx->vscale = 1.0;
710: /*
711: Create the empty data structure to contain compute-h routines.
712: These will be filled in below from the command line options or
713: a later call with MatMFFDSetType() or if that is not called
714: then it will default in the first use of MatMult_MFFD()
715: */
716: mfctx->ops->compute = 0;
717: mfctx->ops->destroy = 0;
718: mfctx->ops->view = 0;
719: mfctx->ops->setfromoptions = 0;
720: mfctx->hctx = 0;
722: mfctx->func = 0;
723: mfctx->funcctx = 0;
724: mfctx->w = NULL;
726: A->data = mfctx;
728: A->ops->mult = MatMult_MFFD;
729: A->ops->destroy = MatDestroy_MFFD;
730: A->ops->view = MatView_MFFD;
731: A->ops->assemblyend = MatAssemblyEnd_MFFD;
732: A->ops->getdiagonal = MatGetDiagonal_MFFD;
733: A->ops->scale = MatScale_MFFD;
734: A->ops->shift = MatShift_MFFD;
735: A->ops->diagonalscale = MatDiagonalScale_MFFD;
736: A->ops->diagonalset = MatDiagonalSet_MFFD;
737: A->ops->setfromoptions = MatSetFromOptions_MFFD;
738: A->assembled = PETSC_TRUE;
740: PetscLayoutSetUp(A->rmap);
741: PetscLayoutSetUp(A->cmap);
743: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
744: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
745: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
746: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
747: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
748: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
749: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
750: PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
751: PetscObjectComposeFunction((PetscObject)A,"MatMFFDAddNullSpace_C",MatMFFDAddNullSpace_MFFD);
753: mfctx->mat = A;
755: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
756: return(0);
757: }
761: /*@
762: MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
764: Collective on Vec
766: Input Parameters:
767: + comm - MPI communicator
768: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
769: This value should be the same as the local size used in creating the
770: y vector for the matrix-vector product y = Ax.
771: . n - This value should be the same as the local size used in creating the
772: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
773: calculated if N is given) For square matrices n is almost always m.
774: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
775: - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
778: Output Parameter:
779: . J - the matrix-free matrix
781: Options Database Keys: call MatSetFromOptions() to trigger these
782: + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
783: - -mat_mffd_err - square root of estimated relative error in function evaluation
784: - -mat_mffd_period - how often h is recomputed, defaults to 1, everytime
787: Level: advanced
789: Notes:
790: The matrix-free matrix context merely contains the function pointers
791: and work space for performing finite difference approximations of
792: Jacobian-vector products, F'(u)*a,
794: The default code uses the following approach to compute h
796: .vb
797: F'(u)*a = [F(u+h*a) - F(u)]/h where
798: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
799: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
800: where
801: error_rel = square root of relative error in function evaluation
802: umin = minimum iterate parameter
803: .ve
805: You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
806: preconditioner matrix
808: The user can set the error_rel via MatMFFDSetFunctionError() and
809: umin via MatMFFDDSSetUmin(); see Users-Manual: ch_snes for details.
811: The user should call MatDestroy() when finished with the matrix-free
812: matrix context.
814: Options Database Keys:
815: + -mat_mffd_err <error_rel> - Sets error_rel
816: . -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
817: - -mat_mffd_check_positivity
819: .keywords: default, matrix-free, create, matrix
821: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
822: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
823: MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
825: @*/
826: PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
827: {
831: MatCreate(comm,J);
832: MatSetSizes(*J,m,n,M,N);
833: MatSetType(*J,MATMFFD);
834: MatSetUp(*J);
835: return(0);
836: }
841: /*@
842: MatMFFDGetH - Gets the last value that was used as the differencing
843: parameter.
845: Not Collective
847: Input Parameters:
848: . mat - the matrix obtained with MatCreateSNESMF()
850: Output Paramter:
851: . h - the differencing step size
853: Level: advanced
855: .keywords: SNES, matrix-free, parameters
857: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
858: @*/
859: PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
860: {
861: MatMFFD ctx = (MatMFFD)mat->data;
863: PetscBool match;
866: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
867: if (!match) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
869: *h = ctx->currenth;
870: return(0);
871: }
875: /*@C
876: MatMFFDSetFunction - Sets the function used in applying the matrix free.
878: Logically Collective on Mat
880: Input Parameters:
881: + mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
882: . func - the function to use
883: - funcctx - optional function context passed to function
885: Calling Sequence of func:
886: $ func (void *funcctx, Vec x, Vec f)
888: + funcctx - user provided context
889: . x - input vector
890: - f - computed output function
892: Level: advanced
894: Notes:
895: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
896: matrix inside your compute Jacobian routine
898: If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
900: .keywords: SNES, matrix-free, function
902: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
903: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
904: @*/
905: PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
906: {
910: PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
911: return(0);
912: }
916: /*@C
917: MatMFFDSetFunctioni - Sets the function for a single component
919: Logically Collective on Mat
921: Input Parameters:
922: + mat - the matrix free matrix created via MatCreateSNESMF()
923: - funci - the function to use
925: Level: advanced
927: Notes:
928: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
929: matrix inside your compute Jacobian routine
932: .keywords: SNES, matrix-free, function
934: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
936: @*/
937: PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
938: {
943: PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
944: return(0);
945: }
950: /*@C
951: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
953: Logically Collective on Mat
955: Input Parameters:
956: + mat - the matrix free matrix created via MatCreateSNESMF()
957: - func - the function to use
959: Level: advanced
961: Notes:
962: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
963: matrix inside your compute Jacobian routine
966: .keywords: SNES, matrix-free, function
968: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
969: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
970: @*/
971: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
972: {
977: PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
978: return(0);
979: }
983: /*@
984: MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime
986: Logically Collective on Mat
988: Input Parameters:
989: + mat - the matrix free matrix created via MatCreateSNESMF()
990: - period - 1 for everytime, 2 for every second etc
992: Options Database Keys:
993: + -mat_mffd_period <period>
995: Level: advanced
998: .keywords: SNES, matrix-free, parameters
1000: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
1001: MatMFFDSetHHistory(), MatMFFDResetHHistory()
1002: @*/
1003: PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
1004: {
1008: PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
1009: return(0);
1010: }
1014: /*@
1015: MatMFFDSetFunctionError - Sets the error_rel for the approximation of
1016: matrix-vector products using finite differences.
1018: Logically Collective on Mat
1020: Input Parameters:
1021: + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
1022: - error_rel - relative error (should be set to the square root of
1023: the relative error in the function evaluations)
1025: Options Database Keys:
1026: + -mat_mffd_err <error_rel> - Sets error_rel
1028: Level: advanced
1030: Notes:
1031: The default matrix-free matrix-vector product routine computes
1032: .vb
1033: F'(u)*a = [F(u+h*a) - F(u)]/h where
1034: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
1035: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
1036: .ve
1038: .keywords: SNES, matrix-free, parameters
1040: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
1041: MatMFFDSetHHistory(), MatMFFDResetHHistory()
1042: @*/
1043: PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
1044: {
1048: PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1049: return(0);
1050: }
1054: /*@
1055: MatMFFDAddNullSpace - Provides a null space that an operator is
1056: supposed to have. Since roundoff will create a small component in
1057: the null space, if you know the null space you may have it
1058: automatically removed.
1060: Logically Collective on Mat
1062: Input Parameters:
1063: + J - the matrix-free matrix context
1064: - nullsp - object created with MatNullSpaceCreate()
1066: Level: advanced
1068: .keywords: SNES, matrix-free, null space
1070: .seealso: MatNullSpaceCreate(), MatMFFDGetH(), MatCreateSNESMF(), MatCreateMFFD(), MATMFFD
1071: MatMFFDSetHHistory(), MatMFFDResetHHistory()
1072: @*/
1073: PetscErrorCode MatMFFDAddNullSpace(Mat J,MatNullSpace nullsp)
1074: {
1078: PetscTryMethod(J,"MatMFFDAddNullSpace_C",(Mat,MatNullSpace),(J,nullsp));
1079: return(0);
1080: }
1084: /*@
1085: MatMFFDSetHHistory - Sets an array to collect a history of the
1086: differencing values (h) computed for the matrix-free product.
1088: Logically Collective on Mat
1090: Input Parameters:
1091: + J - the matrix-free matrix context
1092: . histroy - space to hold the history
1093: - nhistory - number of entries in history, if more entries are generated than
1094: nhistory, then the later ones are discarded
1096: Level: advanced
1098: Notes:
1099: Use MatMFFDResetHHistory() to reset the history counter and collect
1100: a new batch of differencing parameters, h.
1102: .keywords: SNES, matrix-free, h history, differencing history
1104: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1105: MatMFFDResetHHistory(), MatMFFDSetFunctionError()
1107: @*/
1108: PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1109: {
1110: MatMFFD ctx = (MatMFFD)J->data;
1112: PetscBool match;
1115: PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1116: if (!match) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1117: ctx->historyh = history;
1118: ctx->maxcurrenth = nhistory;
1119: ctx->currenth = 0.;
1120: return(0);
1121: }
1126: /*@
1127: MatMFFDResetHHistory - Resets the counter to zero to begin
1128: collecting a new set of differencing histories.
1130: Logically Collective on Mat
1132: Input Parameters:
1133: . J - the matrix-free matrix context
1135: Level: advanced
1137: Notes:
1138: Use MatMFFDSetHHistory() to create the original history counter.
1140: .keywords: SNES, matrix-free, h history, differencing history
1142: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1143: MatMFFDSetHHistory(), MatMFFDSetFunctionError()
1145: @*/
1146: PetscErrorCode MatMFFDResetHHistory(Mat J)
1147: {
1151: PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1152: return(0);
1153: }
1158: /*@
1159: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1160: Jacobian are computed
1162: Logically Collective on Mat
1164: Input Parameters:
1165: + J - the MatMFFD matrix
1166: . U - the vector
1167: - F - (optional) vector that contains F(u) if it has been already computed
1169: Notes: This is rarely used directly
1171: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
1172: point during the first MatMult() after each call to MatMFFDSetBase().
1174: Level: advanced
1176: @*/
1177: PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
1178: {
1185: PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1186: return(0);
1187: }
1191: /*@C
1192: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1193: it to satisfy some criteria
1195: Logically Collective on Mat
1197: Input Parameters:
1198: + J - the MatMFFD matrix
1199: . fun - the function that checks h
1200: - ctx - any context needed by the function
1202: Options Database Keys:
1203: . -mat_mffd_check_positivity
1205: Level: advanced
1207: Notes: For example, MatMFFDSetCheckPositivity() insures that all entries
1208: of U + h*a are non-negative
1210: .seealso: MatMFFDSetCheckPositivity()
1211: @*/
1212: PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1213: {
1218: PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1219: return(0);
1220: }
1224: /*@
1225: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1226: zero, decreases h until this is satisfied.
1228: Logically Collective on Vec
1230: Input Parameters:
1231: + U - base vector that is added to
1232: . a - vector that is added
1233: . h - scaling factor on a
1234: - dummy - context variable (unused)
1236: Options Database Keys:
1237: . -mat_mffd_check_positivity
1239: Level: advanced
1241: Notes: This is rarely used directly, rather it is passed as an argument to
1242: MatMFFDSetCheckh()
1244: .seealso: MatMFFDSetCheckh()
1245: @*/
1246: PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1247: {
1248: PetscReal val, minval;
1249: PetscScalar *u_vec, *a_vec;
1251: PetscInt i,n;
1252: MPI_Comm comm;
1255: PetscObjectGetComm((PetscObject)U,&comm);
1256: VecGetArray(U,&u_vec);
1257: VecGetArray(a,&a_vec);
1258: VecGetLocalSize(U,&n);
1259: minval = PetscAbsScalar(*h*1.01);
1260: for (i=0; i<n; i++) {
1261: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1262: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1263: if (val < minval) minval = val;
1264: }
1265: }
1266: VecRestoreArray(U,&u_vec);
1267: VecRestoreArray(a,&a_vec);
1268: MPI_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1269: if (val <= PetscAbsScalar(*h)) {
1270: PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1271: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1272: else *h = -0.99*val;
1273: }
1274: return(0);
1275: }