Actual source code: mffd.c
petsc-3.7.3 2016-08-01
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,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,NULL, "-log_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: }
209: /* ----------------------------------------------------------------------------------------*/
212: PetscErrorCode MatDestroy_MFFD(Mat mat)
213: {
215: MatMFFD ctx = (MatMFFD)mat->data;
218: VecDestroy(&ctx->w);
219: VecDestroy(&ctx->drscale);
220: VecDestroy(&ctx->dlscale);
221: VecDestroy(&ctx->dshift);
222: if (ctx->current_f_allocated) {
223: VecDestroy(&ctx->current_f);
224: }
225: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
226: PetscHeaderDestroy(&ctx);
227: mat->data = 0;
229: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetBase_C",NULL);
230: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioniBase_C",NULL);
231: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctioni_C",NULL);
232: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunction_C",NULL);
233: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetFunctionError_C",NULL);
234: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetCheckh_C",NULL);
235: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDSetPeriod_C",NULL);
236: PetscObjectComposeFunction((PetscObject)mat,"MatMFFDResetHHistory_C",NULL);
237: return(0);
238: }
242: /*
243: MatMFFDView_MFFD - Views matrix-free parameters.
245: */
246: PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
247: {
249: MatMFFD ctx = (MatMFFD)J->data;
250: PetscBool iascii, viewbase, viewfunction;
251: const char *prefix;
254: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
255: if (iascii) {
256: PetscViewerASCIIPrintf(viewer,"Matrix-free approximation:\n");
257: PetscViewerASCIIPushTab(viewer);
258: PetscViewerASCIIPrintf(viewer,"err=%g (relative error in function evaluation)\n",(double)ctx->error_rel);
259: if (!((PetscObject)ctx)->type_name) {
260: PetscViewerASCIIPrintf(viewer,"The compute h routine has not yet been set\n");
261: } else {
262: PetscViewerASCIIPrintf(viewer,"Using %s compute h routine\n",((PetscObject)ctx)->type_name);
263: }
264: if (ctx->ops->view) {
265: (*ctx->ops->view)(ctx,viewer);
266: }
267: PetscObjectGetOptionsPrefix((PetscObject)J, &prefix);
269: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_base", &viewbase);
270: if (viewbase) {
271: PetscViewerASCIIPrintf(viewer, "Base:\n");
272: VecView(ctx->current_u, viewer);
273: }
274: PetscOptionsHasName(((PetscObject)J)->options,prefix, "-mat_mffd_view_function", &viewfunction);
275: if (viewfunction) {
276: PetscViewerASCIIPrintf(viewer, "Function:\n");
277: VecView(ctx->current_f, viewer);
278: }
279: PetscViewerASCIIPopTab(viewer);
280: }
281: return(0);
282: }
286: /*
287: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
288: allows the user to indicate the beginning of a new linear solve by calling
289: MatAssemblyXXX() on the matrix free matrix. This then allows the
290: MatCreateMFFD_WP() to properly compute ||U|| only the first time
291: in the linear solver rather than every time.
293: This function is referenced directly from MatAssemblyEnd_SNESMF(), which may be in a different shared library hence
294: it must be labeled as PETSC_EXTERN
295: */
296: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
297: {
299: MatMFFD j = (MatMFFD)J->data;
302: MatMFFDResetHHistory(J);
303: j->vshift = 0.0;
304: j->vscale = 1.0;
305: return(0);
306: }
310: /*
311: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
313: y ~= (F(u + ha) - F(u))/h,
314: where F = nonlinear function, as set by SNESSetFunction()
315: u = current iterate
316: h = difference interval
317: */
318: PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
319: {
320: MatMFFD ctx = (MatMFFD)mat->data;
321: PetscScalar h;
322: Vec w,U,F;
324: PetscBool zeroa;
327: 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");
328: /* We log matrix-free matrix-vector products separately, so that we can
329: separate the performance monitoring from the cases that use conventional
330: storage. We may eventually modify event logging to associate events
331: with particular objects, hence alleviating the more general problem. */
332: PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);
334: w = ctx->w;
335: U = ctx->current_u;
336: F = ctx->current_f;
337: /*
338: Compute differencing parameter
339: */
340: if (!ctx->ops->compute) {
341: MatMFFDSetType(mat,MATMFFD_WP);
342: MatSetFromOptions(mat);
343: }
344: (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
345: if (zeroa) {
346: VecSet(y,0.0);
347: return(0);
348: }
350: if (mat->erroriffailure && PetscIsInfOrNanScalar(h)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed Nan differencing parameter h");
351: if (ctx->checkh) {
352: (*ctx->checkh)(ctx->checkhctx,U,a,&h);
353: }
355: /* keep a record of the current differencing parameter h */
356: ctx->currenth = h;
357: #if defined(PETSC_USE_COMPLEX)
358: PetscInfo2(mat,"Current differencing parameter: %g + %g i\n",(double)PetscRealPart(h),(double)PetscImaginaryPart(h));
359: #else
360: PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
361: #endif
362: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
363: ctx->historyh[ctx->ncurrenth] = h;
364: }
365: ctx->ncurrenth++;
367: /* w = u + ha */
368: if (ctx->drscale) {
369: VecPointwiseMult(ctx->drscale,a,U);
370: VecAYPX(U,h,w);
371: } else {
372: VecWAXPY(w,h,a,U);
373: }
375: /* compute func(U) as base for differencing; only needed first time in and not when provided by user */
376: if (ctx->ncurrenth == 1 && ctx->current_f_allocated) {
377: (*ctx->func)(ctx->funcctx,U,F);
378: }
379: (*ctx->func)(ctx->funcctx,w,y);
381: VecAXPY(y,-1.0,F);
382: VecScale(y,1.0/h);
384: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
385: VecAXPBY(y,ctx->vshift,ctx->vscale,a);
386: }
387: if (ctx->dlscale) {
388: VecPointwiseMult(y,ctx->dlscale,y);
389: }
390: if (ctx->dshift) {
391: VecPointwiseMult(ctx->dshift,a,U);
392: VecAXPY(y,1.0,U);
393: }
395: if (mat->nullsp) {MatNullSpaceRemove(mat->nullsp,y);}
397: PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);
398: return(0);
399: }
403: /*
404: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
406: y ~= (F(u + ha) - F(u))/h,
407: where F = nonlinear function, as set by SNESSetFunction()
408: u = current iterate
409: h = difference interval
410: */
411: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
412: {
413: MatMFFD ctx = (MatMFFD)mat->data;
414: PetscScalar h,*aa,*ww,v;
415: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
416: Vec w,U;
418: PetscInt i,rstart,rend;
421: if (!ctx->funci) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Requires calling MatMFFDSetFunctioni() first");
423: w = ctx->w;
424: U = ctx->current_u;
425: (*ctx->func)(ctx->funcctx,U,a);
426: (*ctx->funcisetbase)(ctx->funcctx,U);
427: VecCopy(U,w);
429: VecGetOwnershipRange(a,&rstart,&rend);
430: VecGetArray(a,&aa);
431: for (i=rstart; i<rend; i++) {
432: VecGetArray(w,&ww);
433: h = ww[i-rstart];
434: if (h == 0.0) h = 1.0;
435: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
436: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
437: h *= epsilon;
439: ww[i-rstart] += h;
440: VecRestoreArray(w,&ww);
441: (*ctx->funci)(ctx->funcctx,i,w,&v);
442: aa[i-rstart] = (v - aa[i-rstart])/h;
444: /* possibly shift and scale result */
445: if ((ctx->vshift != 0.0) || (ctx->vscale != 1.0)) {
446: aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
447: }
449: VecGetArray(w,&ww);
450: ww[i-rstart] -= h;
451: VecRestoreArray(w,&ww);
452: }
453: VecRestoreArray(a,&aa);
454: return(0);
455: }
459: PetscErrorCode MatDiagonalScale_MFFD(Mat mat,Vec ll,Vec rr)
460: {
461: MatMFFD aij = (MatMFFD)mat->data;
465: if (ll && !aij->dlscale) {
466: VecDuplicate(ll,&aij->dlscale);
467: }
468: if (rr && !aij->drscale) {
469: VecDuplicate(rr,&aij->drscale);
470: }
471: if (ll) {
472: VecCopy(ll,aij->dlscale);
473: }
474: if (rr) {
475: VecCopy(rr,aij->drscale);
476: }
477: return(0);
478: }
482: PetscErrorCode MatDiagonalSet_MFFD(Mat mat,Vec ll,InsertMode mode)
483: {
484: MatMFFD aij = (MatMFFD)mat->data;
488: if (mode == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No diagonal set with INSERT_VALUES");
489: if (!aij->dshift) {
490: VecDuplicate(ll,&aij->dshift);
491: }
492: VecAXPY(aij->dshift,1.0,ll);
493: return(0);
494: }
498: PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
499: {
500: MatMFFD shell = (MatMFFD)Y->data;
503: shell->vshift += a;
504: return(0);
505: }
509: PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
510: {
511: MatMFFD shell = (MatMFFD)Y->data;
514: shell->vscale *= a;
515: return(0);
516: }
520: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat J,Vec U,Vec F)
521: {
523: MatMFFD ctx = (MatMFFD)J->data;
526: MatMFFDResetHHistory(J);
528: ctx->current_u = U;
529: if (F) {
530: if (ctx->current_f_allocated) {VecDestroy(&ctx->current_f);}
531: ctx->current_f = F;
532: ctx->current_f_allocated = PETSC_FALSE;
533: } else if (!ctx->current_f_allocated) {
534: MatCreateVecs(J,NULL,&ctx->current_f);
536: ctx->current_f_allocated = PETSC_TRUE;
537: }
538: if (!ctx->w) {
539: VecDuplicate(ctx->current_u, &ctx->w);
540: }
541: J->assembled = PETSC_TRUE;
542: return(0);
543: }
545: typedef PetscErrorCode (*FCN3)(void*,Vec,Vec,PetscScalar*); /* force argument to next function to not be extern C*/
549: PetscErrorCode MatMFFDSetCheckh_MFFD(Mat J,FCN3 fun,void *ectx)
550: {
551: MatMFFD ctx = (MatMFFD)J->data;
554: ctx->checkh = fun;
555: ctx->checkhctx = ectx;
556: return(0);
557: }
561: /*@C
562: MatMFFDSetOptionsPrefix - Sets the prefix used for searching for all
563: MatMFFD options in the database.
565: Collective on Mat
567: Input Parameter:
568: + A - the Mat context
569: - prefix - the prefix to prepend to all option names
571: Notes:
572: A hyphen (-) must NOT be given at the beginning of the prefix name.
573: The first character of all runtime options is AUTOMATICALLY the hyphen.
575: Level: advanced
577: .keywords: SNES, matrix-free, parameters
579: .seealso: MatSetFromOptions(), MatCreateSNESMF()
580: @*/
581: PetscErrorCode MatMFFDSetOptionsPrefix(Mat mat,const char prefix[])
583: {
584: MatMFFD mfctx = mat ? (MatMFFD)mat->data : (MatMFFD)NULL;
590: PetscObjectSetOptionsPrefix((PetscObject)mfctx,prefix);
591: return(0);
592: }
596: PetscErrorCode MatSetFromOptions_MFFD(PetscOptionItems *PetscOptionsObject,Mat mat)
597: {
598: MatMFFD mfctx = (MatMFFD)mat->data;
600: PetscBool flg;
601: char ftype[256];
606: PetscObjectOptionsBegin((PetscObject)mfctx);
607: PetscOptionsFList("-mat_mffd_type","Matrix free type","MatMFFDSetType",MatMFFDList,((PetscObject)mfctx)->type_name,ftype,256,&flg);
608: if (flg) {
609: MatMFFDSetType(mat,ftype);
610: }
612: PetscOptionsReal("-mat_mffd_err","set sqrt relative error in function","MatMFFDSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
613: PetscOptionsInt("-mat_mffd_period","how often h is recomputed","MatMFFDSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
615: flg = PETSC_FALSE;
616: PetscOptionsBool("-mat_mffd_check_positivity","Insure that U + h*a is nonnegative","MatMFFDSetCheckh",flg,&flg,NULL);
617: if (flg) {
618: MatMFFDSetCheckh(mat,MatMFFDCheckPositivity,0);
619: }
620: if (mfctx->ops->setfromoptions) {
621: (*mfctx->ops->setfromoptions)(PetscOptionsObject,mfctx);
622: }
623: PetscOptionsEnd();
624: return(0);
625: }
629: PetscErrorCode MatMFFDSetPeriod_MFFD(Mat mat,PetscInt period)
630: {
631: MatMFFD ctx = (MatMFFD)mat->data;
635: ctx->recomputeperiod = period;
636: return(0);
637: }
641: PetscErrorCode MatMFFDSetFunction_MFFD(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
642: {
643: MatMFFD ctx = (MatMFFD)mat->data;
646: ctx->func = func;
647: ctx->funcctx = funcctx;
648: return(0);
649: }
653: PetscErrorCode MatMFFDSetFunctionError_MFFD(Mat mat,PetscReal error)
654: {
655: MatMFFD ctx = (MatMFFD)mat->data;
659: if (error != PETSC_DEFAULT) ctx->error_rel = error;
660: return(0);
661: }
665: static PetscErrorCode MatMissingDiagonal_MFFD(Mat A,PetscBool *missing,PetscInt *d)
666: {
668: *missing = PETSC_FALSE;
669: return(0);
670: }
672: /*MC
673: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
675: Level: advanced
677: .seealso: MatCreateMFFD(), MatCreateSNESMF(), MatMFFDSetFunction()
678: M*/
681: PETSC_EXTERN PetscErrorCode MatCreate_MFFD(Mat A)
682: {
683: MatMFFD mfctx;
687: MatMFFDInitializePackage();
689: PetscHeaderCreate(mfctx,MATMFFD_CLASSID,"MatMFFD","Matrix-free Finite Differencing","Mat",PetscObjectComm((PetscObject)A),MatDestroy_MFFD,MatView_MFFD);
691: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
692: mfctx->recomputeperiod = 1;
693: mfctx->count = 0;
694: mfctx->currenth = 0.0;
695: mfctx->historyh = NULL;
696: mfctx->ncurrenth = 0;
697: mfctx->maxcurrenth = 0;
698: ((PetscObject)mfctx)->type_name = 0;
700: mfctx->vshift = 0.0;
701: mfctx->vscale = 1.0;
703: /*
704: Create the empty data structure to contain compute-h routines.
705: These will be filled in below from the command line options or
706: a later call with MatMFFDSetType() or if that is not called
707: then it will default in the first use of MatMult_MFFD()
708: */
709: mfctx->ops->compute = 0;
710: mfctx->ops->destroy = 0;
711: mfctx->ops->view = 0;
712: mfctx->ops->setfromoptions = 0;
713: mfctx->hctx = 0;
715: mfctx->func = 0;
716: mfctx->funcctx = 0;
717: mfctx->w = NULL;
719: A->data = mfctx;
721: A->ops->mult = MatMult_MFFD;
722: A->ops->destroy = MatDestroy_MFFD;
723: A->ops->view = MatView_MFFD;
724: A->ops->assemblyend = MatAssemblyEnd_MFFD;
725: A->ops->getdiagonal = MatGetDiagonal_MFFD;
726: A->ops->scale = MatScale_MFFD;
727: A->ops->shift = MatShift_MFFD;
728: A->ops->diagonalscale = MatDiagonalScale_MFFD;
729: A->ops->diagonalset = MatDiagonalSet_MFFD;
730: A->ops->setfromoptions = MatSetFromOptions_MFFD;
731: A->ops->missingdiagonal = MatMissingDiagonal_MFFD;
732: A->assembled = PETSC_TRUE;
734: PetscLayoutSetUp(A->rmap);
735: PetscLayoutSetUp(A->cmap);
737: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C",MatMFFDSetBase_MFFD);
738: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioniBase_C",MatMFFDSetFunctioniBase_MFFD);
739: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctioni_C",MatMFFDSetFunctioni_MFFD);
740: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunction_C",MatMFFDSetFunction_MFFD);
741: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetCheckh_C",MatMFFDSetCheckh_MFFD);
742: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetPeriod_C",MatMFFDSetPeriod_MFFD);
743: PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetFunctionError_C",MatMFFDSetFunctionError_MFFD);
744: PetscObjectComposeFunction((PetscObject)A,"MatMFFDResetHHistory_C",MatMFFDResetHHistory_MFFD);
746: mfctx->mat = A;
748: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
749: return(0);
750: }
754: /*@
755: MatCreateMFFD - Creates a matrix-free matrix. See also MatCreateSNESMF()
757: Collective on Vec
759: Input Parameters:
760: + comm - MPI communicator
761: . m - number of local rows (or PETSC_DECIDE to have calculated if M is given)
762: This value should be the same as the local size used in creating the
763: y vector for the matrix-vector product y = Ax.
764: . n - This value should be the same as the local size used in creating the
765: x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
766: calculated if N is given) For square matrices n is almost always m.
767: . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given)
768: - N - number of global columns (or PETSC_DETERMINE to have calculated if n is given)
771: Output Parameter:
772: . J - the matrix-free matrix
774: Options Database Keys: call MatSetFromOptions() to trigger these
775: + -mat_mffd_type - wp or ds (see MATMFFD_WP or MATMFFD_DS)
776: - -mat_mffd_err - square root of estimated relative error in function evaluation
777: - -mat_mffd_period - how often h is recomputed, defaults to 1, everytime
780: Level: advanced
782: Notes:
783: The matrix-free matrix context merely contains the function pointers
784: and work space for performing finite difference approximations of
785: Jacobian-vector products, F'(u)*a,
787: The default code uses the following approach to compute h
789: .vb
790: F'(u)*a = [F(u+h*a) - F(u)]/h where
791: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
792: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
793: where
794: error_rel = square root of relative error in function evaluation
795: umin = minimum iterate parameter
796: .ve
798: You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
799: preconditioner matrix
801: The user can set the error_rel via MatMFFDSetFunctionError() and
802: umin via MatMFFDDSSetUmin(); see Users-Manual: ch_snes for details.
804: The user should call MatDestroy() when finished with the matrix-free
805: matrix context.
807: Options Database Keys:
808: + -mat_mffd_err <error_rel> - Sets error_rel
809: . -mat_mffd_unim <umin> - Sets umin (for default PETSc routine that computes h only)
810: - -mat_mffd_check_positivity
812: .keywords: default, matrix-free, create, matrix
814: .seealso: MatDestroy(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin(), MatMFFDSetFunction()
815: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateSNESMF(),
816: MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
818: @*/
819: PetscErrorCode MatCreateMFFD(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,Mat *J)
820: {
824: MatCreate(comm,J);
825: MatSetSizes(*J,m,n,M,N);
826: MatSetType(*J,MATMFFD);
827: MatSetUp(*J);
828: return(0);
829: }
834: /*@
835: MatMFFDGetH - Gets the last value that was used as the differencing
836: parameter.
838: Not Collective
840: Input Parameters:
841: . mat - the matrix obtained with MatCreateSNESMF()
843: Output Paramter:
844: . h - the differencing step size
846: Level: advanced
848: .keywords: SNES, matrix-free, parameters
850: .seealso: MatCreateSNESMF(),MatMFFDSetHHistory(), MatCreateMFFD(), MATMFFD, MatMFFDResetHHistory()
851: @*/
852: PetscErrorCode MatMFFDGetH(Mat mat,PetscScalar *h)
853: {
854: MatMFFD ctx = (MatMFFD)mat->data;
856: PetscBool match;
859: PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);
860: if (!match) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
862: *h = ctx->currenth;
863: return(0);
864: }
868: /*@C
869: MatMFFDSetFunction - Sets the function used in applying the matrix free.
871: Logically Collective on Mat
873: Input Parameters:
874: + mat - the matrix free matrix created via MatCreateSNESMF() or MatCreateMFFD()
875: . func - the function to use
876: - funcctx - optional function context passed to function
878: Calling Sequence of func:
879: $ func (void *funcctx, Vec x, Vec f)
881: + funcctx - user provided context
882: . x - input vector
883: - f - computed output function
885: Level: advanced
887: Notes:
888: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
889: matrix inside your compute Jacobian routine
891: If this is not set then it will use the function set with SNESSetFunction() if MatCreateSNESMF() was used.
893: .keywords: SNES, matrix-free, function
895: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD,
896: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
897: @*/
898: PetscErrorCode MatMFFDSetFunction(Mat mat,PetscErrorCode (*func)(void*,Vec,Vec),void *funcctx)
899: {
903: PetscTryMethod(mat,"MatMFFDSetFunction_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec),void*),(mat,func,funcctx));
904: return(0);
905: }
909: /*@C
910: MatMFFDSetFunctioni - Sets the function for a single component
912: Logically Collective on Mat
914: Input Parameters:
915: + mat - the matrix free matrix created via MatCreateSNESMF()
916: - funci - the function to use
918: Level: advanced
920: Notes:
921: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
922: matrix inside your compute Jacobian routine
925: .keywords: SNES, matrix-free, function
927: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
929: @*/
930: PetscErrorCode MatMFFDSetFunctioni(Mat mat,PetscErrorCode (*funci)(void*,PetscInt,Vec,PetscScalar*))
931: {
936: PetscTryMethod(mat,"MatMFFDSetFunctioni_C",(Mat,PetscErrorCode (*)(void*,PetscInt,Vec,PetscScalar*)),(mat,funci));
937: return(0);
938: }
943: /*@C
944: MatMFFDSetFunctioniBase - Sets the base vector for a single component function evaluation
946: Logically Collective on Mat
948: Input Parameters:
949: + mat - the matrix free matrix created via MatCreateSNESMF()
950: - func - the function to use
952: Level: advanced
954: Notes:
955: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
956: matrix inside your compute Jacobian routine
959: .keywords: SNES, matrix-free, function
961: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
962: MatMFFDSetHHistory(), MatMFFDResetHHistory(), SNESetFunction()
963: @*/
964: PetscErrorCode MatMFFDSetFunctioniBase(Mat mat,PetscErrorCode (*func)(void*,Vec))
965: {
970: PetscTryMethod(mat,"MatMFFDSetFunctioniBase_C",(Mat,PetscErrorCode (*)(void*,Vec)),(mat,func));
971: return(0);
972: }
976: /*@
977: MatMFFDSetPeriod - Sets how often h is recomputed, by default it is everytime
979: Logically Collective on Mat
981: Input Parameters:
982: + mat - the matrix free matrix created via MatCreateSNESMF()
983: - period - 1 for everytime, 2 for every second etc
985: Options Database Keys:
986: + -mat_mffd_period <period>
988: Level: advanced
991: .keywords: SNES, matrix-free, parameters
993: .seealso: MatCreateSNESMF(),MatMFFDGetH(),
994: MatMFFDSetHHistory(), MatMFFDResetHHistory()
995: @*/
996: PetscErrorCode MatMFFDSetPeriod(Mat mat,PetscInt period)
997: {
1001: PetscTryMethod(mat,"MatMFFDSetPeriod_C",(Mat,PetscInt),(mat,period));
1002: return(0);
1003: }
1007: /*@
1008: MatMFFDSetFunctionError - Sets the error_rel for the approximation of
1009: matrix-vector products using finite differences.
1011: Logically Collective on Mat
1013: Input Parameters:
1014: + mat - the matrix free matrix created via MatCreateMFFD() or MatCreateSNESMF()
1015: - error_rel - relative error (should be set to the square root of
1016: the relative error in the function evaluations)
1018: Options Database Keys:
1019: + -mat_mffd_err <error_rel> - Sets error_rel
1021: Level: advanced
1023: Notes:
1024: The default matrix-free matrix-vector product routine computes
1025: .vb
1026: F'(u)*a = [F(u+h*a) - F(u)]/h where
1027: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
1028: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
1029: .ve
1031: .keywords: SNES, matrix-free, parameters
1033: .seealso: MatCreateSNESMF(),MatMFFDGetH(), MatCreateMFFD(), MATMFFD
1034: MatMFFDSetHHistory(), MatMFFDResetHHistory()
1035: @*/
1036: PetscErrorCode MatMFFDSetFunctionError(Mat mat,PetscReal error)
1037: {
1041: PetscTryMethod(mat,"MatMFFDSetFunctionError_C",(Mat,PetscReal),(mat,error));
1042: return(0);
1043: }
1047: /*@
1048: MatMFFDSetHHistory - Sets an array to collect a history of the
1049: differencing values (h) computed for the matrix-free product.
1051: Logically Collective on Mat
1053: Input Parameters:
1054: + J - the matrix-free matrix context
1055: . histroy - space to hold the history
1056: - nhistory - number of entries in history, if more entries are generated than
1057: nhistory, then the later ones are discarded
1059: Level: advanced
1061: Notes:
1062: Use MatMFFDResetHHistory() to reset the history counter and collect
1063: a new batch of differencing parameters, h.
1065: .keywords: SNES, matrix-free, h history, differencing history
1067: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1068: MatMFFDResetHHistory(), MatMFFDSetFunctionError()
1070: @*/
1071: PetscErrorCode MatMFFDSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1072: {
1073: MatMFFD ctx = (MatMFFD)J->data;
1075: PetscBool match;
1078: PetscObjectTypeCompare((PetscObject)J,MATMFFD,&match);
1079: if (!match) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_WRONG,"Not a MFFD matrix");
1080: ctx->historyh = history;
1081: ctx->maxcurrenth = nhistory;
1082: ctx->currenth = 0.;
1083: return(0);
1084: }
1089: /*@
1090: MatMFFDResetHHistory - Resets the counter to zero to begin
1091: collecting a new set of differencing histories.
1093: Logically Collective on Mat
1095: Input Parameters:
1096: . J - the matrix-free matrix context
1098: Level: advanced
1100: Notes:
1101: Use MatMFFDSetHHistory() to create the original history counter.
1103: .keywords: SNES, matrix-free, h history, differencing history
1105: .seealso: MatMFFDGetH(), MatCreateSNESMF(),
1106: MatMFFDSetHHistory(), MatMFFDSetFunctionError()
1108: @*/
1109: PetscErrorCode MatMFFDResetHHistory(Mat J)
1110: {
1114: PetscTryMethod(J,"MatMFFDResetHHistory_C",(Mat),(J));
1115: return(0);
1116: }
1121: /*@
1122: MatMFFDSetBase - Sets the vector U at which matrix vector products of the
1123: Jacobian are computed
1125: Logically Collective on Mat
1127: Input Parameters:
1128: + J - the MatMFFD matrix
1129: . U - the vector
1130: - F - (optional) vector that contains F(u) if it has been already computed
1132: Notes: This is rarely used directly
1134: If F is provided then it is not recomputed. Otherwise the function is evaluated at the base
1135: point during the first MatMult() after each call to MatMFFDSetBase().
1137: Level: advanced
1139: @*/
1140: PetscErrorCode MatMFFDSetBase(Mat J,Vec U,Vec F)
1141: {
1148: PetscTryMethod(J,"MatMFFDSetBase_C",(Mat,Vec,Vec),(J,U,F));
1149: return(0);
1150: }
1154: /*@C
1155: MatMFFDSetCheckh - Sets a function that checks the computed h and adjusts
1156: it to satisfy some criteria
1158: Logically Collective on Mat
1160: Input Parameters:
1161: + J - the MatMFFD matrix
1162: . fun - the function that checks h
1163: - ctx - any context needed by the function
1165: Options Database Keys:
1166: . -mat_mffd_check_positivity
1168: Level: advanced
1170: Notes: For example, MatMFFDSetCheckPositivity() insures that all entries
1171: of U + h*a are non-negative
1173: .seealso: MatMFFDSetCheckPositivity()
1174: @*/
1175: PetscErrorCode MatMFFDSetCheckh(Mat J,PetscErrorCode (*fun)(void*,Vec,Vec,PetscScalar*),void *ctx)
1176: {
1181: PetscTryMethod(J,"MatMFFDSetCheckh_C",(Mat,PetscErrorCode (*)(void*,Vec,Vec,PetscScalar*),void*),(J,fun,ctx));
1182: return(0);
1183: }
1187: /*@
1188: MatMFFDCheckPositivity - Checks that all entries in U + h*a are positive or
1189: zero, decreases h until this is satisfied.
1191: Logically Collective on Vec
1193: Input Parameters:
1194: + U - base vector that is added to
1195: . a - vector that is added
1196: . h - scaling factor on a
1197: - dummy - context variable (unused)
1199: Options Database Keys:
1200: . -mat_mffd_check_positivity
1202: Level: advanced
1204: Notes: This is rarely used directly, rather it is passed as an argument to
1205: MatMFFDSetCheckh()
1207: .seealso: MatMFFDSetCheckh()
1208: @*/
1209: PetscErrorCode MatMFFDCheckPositivity(void *dummy,Vec U,Vec a,PetscScalar *h)
1210: {
1211: PetscReal val, minval;
1212: PetscScalar *u_vec, *a_vec;
1214: PetscInt i,n;
1215: MPI_Comm comm;
1218: PetscObjectGetComm((PetscObject)U,&comm);
1219: VecGetArray(U,&u_vec);
1220: VecGetArray(a,&a_vec);
1221: VecGetLocalSize(U,&n);
1222: minval = PetscAbsScalar(*h*1.01);
1223: for (i=0; i<n; i++) {
1224: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1225: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1226: if (val < minval) minval = val;
1227: }
1228: }
1229: VecRestoreArray(U,&u_vec);
1230: VecRestoreArray(a,&a_vec);
1231: MPIU_Allreduce(&minval,&val,1,MPIU_REAL,MPIU_MIN,comm);
1232: if (val <= PetscAbsScalar(*h)) {
1233: PetscInfo2(U,"Scaling back h from %g to %g\n",(double)PetscRealPart(*h),(double)(.99*val));
1234: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1235: else *h = -0.99*val;
1236: }
1237: return(0);
1238: }