Actual source code: snesmfj.c
1: /*$Id: snesmfj.c,v 1.131 2001/09/05 18:45:40 bsmith Exp $*/
3: #include src/mat/matimpl.h
4: #include src/snes/mf/snesmfj.h
6: PetscFList MatSNESMPetscFList = 0;
7: PetscTruth MatSNESMFRegisterAllCalled = PETSC_FALSE;
9: #undef __FUNCT__
11: /*@C
12: MatSNESMFSetType - Sets the method that is used to compute the
13: differencing parameter for finite differene matrix-free formulations.
15: Input Parameters:
16: + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMF()
17: or MatSetType(mat,MATMFFD);
18: - ftype - the type requested
20: Level: advanced
22: Notes:
23: For example, such routines can compute h for use in
24: Jacobian-vector products of the form
26: F(x+ha) - F(x)
27: F'(u)a ~= ----------------
28: h
30: .seealso: MatCreateSNESMF(), MatSNESMFRegisterDynamic)
31: @*/
32: int MatSNESMFSetType(Mat mat,MatSNESMFType ftype)
33: {
34: int ierr,(*r)(MatSNESMFCtx);
35: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
36: PetscTruth match;
37:
42: /* already set, so just return */
43: PetscTypeCompare((PetscObject)ctx,ftype,&match);
44: if (match) return(0);
46: /* destroy the old one if it exists */
47: if (ctx->ops->destroy) {
48: (*ctx->ops->destroy)(ctx);
49: }
51: /* Get the function pointers for the requrested method */
52: if (!MatSNESMFRegisterAllCalled) {MatSNESMFRegisterAll(PETSC_NULL);}
54: PetscFListFind(ctx->comm,MatSNESMPetscFList,ftype,(void (**)(void)) &r);
56: if (!r) SETERRQ(1,"Unknown MatSNESMF type given");
58: (*r)(ctx);
60: PetscObjectChangeTypeName((PetscObject)ctx,ftype);
62: return(0);
63: }
65: EXTERN_C_BEGIN
66: #undef __FUNCT__
68: int MatSNESMFSetFunctioniBase_FD(Mat mat,int (*func)(Vec,void *))
69: {
70: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
73: ctx->funcisetbase = func;
74: return(0);
75: }
76: EXTERN_C_END
78: EXTERN_C_BEGIN
79: #undef __FUNCT__
81: int MatSNESMFSetFunctioni_FD(Mat mat,int (*funci)(int,Vec,PetscScalar*,void *))
82: {
83: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
86: ctx->funci = funci;
87: return(0);
88: }
89: EXTERN_C_END
91: /*MC
92: MatSNESMFRegisterDynamic - Adds a method to the MatSNESMF registry.
94: Synopsis:
95: int MatSNESMFRegisterDynamic(char *name_solver,char *path,char *name_create,int (*routine_create)(MatSNESMF))
97: Not Collective
99: Input Parameters:
100: + name_solver - name of a new user-defined compute-h module
101: . path - path (either absolute or relative) the library containing this solver
102: . name_create - name of routine to create method context
103: - routine_create - routine to create method context
105: Level: developer
107: Notes:
108: MatSNESMFRegisterDynamic) may be called multiple times to add several user-defined solvers.
110: If dynamic libraries are used, then the fourth input argument (routine_create)
111: is ignored.
113: Sample usage:
114: .vb
115: MatSNESMFRegisterDynamic"my_h",/home/username/my_lib/lib/libO/solaris/mylib.a,
116: "MyHCreate",MyHCreate);
117: .ve
119: Then, your solver can be chosen with the procedural interface via
120: $ MatSNESMFSetType(mfctx,"my_h")
121: or at runtime via the option
122: $ -snes_mf_type my_h
124: .keywords: MatSNESMF, register
126: .seealso: MatSNESMFRegisterAll(), MatSNESMFRegisterDestroy()
127: M*/
129: #undef __FUNCT__
131: int MatSNESMFRegister(char *sname,char *path,char *name,int (*function)(MatSNESMFCtx))
132: {
134: char fullname[256];
137: PetscFListConcat(path,name,fullname);
138: PetscFListAdd(&MatSNESMPetscFList,sname,fullname,(void (*)(void))function);
139: return(0);
140: }
143: #undef __FUNCT__
145: /*@C
146: MatSNESMFRegisterDestroy - Frees the list of MatSNESMF methods that were
147: registered by MatSNESMFRegisterDynamic).
149: Not Collective
151: Level: developer
153: .keywords: MatSNESMF, register, destroy
155: .seealso: MatSNESMFRegisterDynamic), MatSNESMFRegisterAll()
156: @*/
157: int MatSNESMFRegisterDestroy(void)
158: {
162: if (MatSNESMPetscFList) {
163: PetscFListDestroy(&MatSNESMPetscFList);
164: MatSNESMPetscFList = 0;
165: }
166: MatSNESMFRegisterAllCalled = PETSC_FALSE;
167: return(0);
168: }
170: /* ----------------------------------------------------------------------------------------*/
171: #undef __FUNCT__
173: int MatDestroy_MFFD(Mat mat)
174: {
175: int ierr;
176: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
179: if (ctx->w != PETSC_NULL) {
180: VecDestroy(ctx->w);
181: }
182: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
183: if (ctx->sp) {MatNullSpaceDestroy(ctx->sp);}
184: PetscHeaderDestroy(ctx);
185: return(0);
186: }
188: #undef __FUNCT__
190: /*
191: MatSNESMFView_MFFD - Views matrix-free parameters.
193: */
194: int MatView_MFFD(Mat J,PetscViewer viewer)
195: {
196: int ierr;
197: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
198: PetscTruth isascii;
201: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
202: if (isascii) {
203: PetscViewerASCIIPrintf(viewer," SNES matrix-free approximation:n");
204: PetscViewerASCIIPrintf(viewer," err=%g (relative error in function evaluation)n",ctx->error_rel);
205: if (!ctx->type_name) {
206: PetscViewerASCIIPrintf(viewer," The compute h routine has not yet been setn");
207: } else {
208: PetscViewerASCIIPrintf(viewer," Using %s compute h routinen",ctx->type_name);
209: }
210: if (ctx->ops->view) {
211: (*ctx->ops->view)(ctx,viewer);
212: }
213: } else {
214: SETERRQ1(1,"Viewer type %s not supported for SNES matrix free matrix",((PetscObject)viewer)->type_name);
215: }
216: return(0);
217: }
219: #undef __FUNCT__
221: /*
222: MatSNESMFAssemblyEnd_Private - Resets the ctx->ncurrenth to zero. This
223: allows the user to indicate the beginning of a new linear solve by calling
224: MatAssemblyXXX() on the matrix free matrix. This then allows the
225: MatSNESMFCreate_WP() to properly compute ||U|| only the first time
226: in the linear solver rather than every time.
227: */
228: int MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
229: {
230: int ierr;
231: MatSNESMFCtx j = (MatSNESMFCtx)J->data;
234: MatSNESMFResetHHistory(J);
235: if (j->usesnes) {
236: SNESGetSolution(j->snes,&j->current_u);
237: SNESGetFunction(j->snes,&j->current_f,PETSC_NULL,PETSC_NULL);
238: if (j->w == PETSC_NULL) {
239: VecDuplicate(j->current_u, &j->w);
240: }
241: }
242: j->vshift = 0.0;
243: j->vscale = 1.0;
244: return(0);
245: }
247: #undef __FUNCT__
249: /*
250: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
252: y ~= (F(u + ha) - F(u))/h,
253: where F = nonlinear function, as set by SNESSetFunction()
254: u = current iterate
255: h = difference interval
256: */
257: int MatMult_MFFD(Mat mat,Vec a,Vec y)
258: {
259: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
260: SNES snes;
261: PetscScalar h,mone = -1.0;
262: Vec w,U,F;
263: int ierr,(*eval_fct)(SNES,Vec,Vec)=0;
266: /* We log matrix-free matrix-vector products separately, so that we can
267: separate the performance monitoring from the cases that use conventional
268: storage. We may eventually modify event logging to associate events
269: with particular objects, hence alleviating the more general problem. */
270: PetscLogEventBegin(MAT_MultMatrixFree,a,y,0,0);
272: snes = ctx->snes;
273: w = ctx->w;
274: U = ctx->current_u;
276: /*
277: Compute differencing parameter
278: */
279: if (!ctx->ops->compute) {
280: MatSNESMFSetType(mat,MATSNESMF_DEFAULT);
281: MatSNESMFSetFromOptions(mat);
282: }
283: (*ctx->ops->compute)(ctx,U,a,&h);
285: if (ctx->checkh) {
286: (*ctx->checkh)(U,a,&h,ctx->checkhctx);
287: }
289: /* keep a record of the current differencing parameter h */
290: ctx->currenth = h;
291: #if defined(PETSC_USE_COMPLEX)
292: PetscLogInfo(mat,"MatMult_MFFD:Current differencing parameter: %g + %g in",PetscRealPart(h),PetscImaginaryPart(h));
293: #else
294: PetscLogInfo(mat,"MatMult_MFFD:Current differencing parameter: %15.12en",h);
295: #endif
296: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
297: ctx->historyh[ctx->ncurrenth] = h;
298: }
299: ctx->ncurrenth++;
301: /* w = u + ha */
302: VecWAXPY(&h,a,U,w);
304: if (ctx->usesnes) {
305: eval_fct = SNESComputeFunction;
306: F = ctx->current_f;
307: if (!F) SETERRQ(1,"You must call MatAssembly() even on matrix-free matrices");
308: (*eval_fct)(snes,w,y);
309: } else {
310: F = ctx->funcvec;
311: /* compute func(U) as base for differencing */
312: if (ctx->ncurrenth == 1) {
313: (*ctx->func)(snes,U,F,ctx->funcctx);
314: }
315: (*ctx->func)(snes,w,y,ctx->funcctx);
316: }
318: VecAXPY(&mone,F,y);
319: h = 1.0/h;
320: VecScale(&h,y);
323: if (ctx->vshift != 0.0 && ctx->vscale != 1.0) {
324: VecAXPBY(&ctx->vshift,&ctx->vscale,a,y);
325: } else if (ctx->vscale != 1.0) {
326: VecScale(&ctx->vscale,y);
327: } else if (ctx->vshift != 0.0) {
328: VecAXPY(&ctx->vshift,a,y);
329: }
331: if (ctx->sp) {MatNullSpaceRemove(ctx->sp,y,PETSC_NULL);}
333: PetscLogEventEnd(MAT_MultMatrixFree,a,y,0,0);
334: return(0);
335: }
337: #undef __FUNCT__
339: /*
340: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
342: y ~= (F(u + ha) - F(u))/h,
343: where F = nonlinear function, as set by SNESSetFunction()
344: u = current iterate
345: h = difference interval
346: */
347: int MatGetDiagonal_MFFD(Mat mat,Vec a)
348: {
349: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
350: PetscScalar h,*aa,*ww,v;
351: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
352: Vec w,U;
353: int i,ierr,rstart,rend;
356: if (!ctx->funci) {
357: SETERRQ(1,"Requirers calling MatSNESMFSetFunctioni() first");
358: }
360: w = ctx->w;
361: U = ctx->current_u;
362: (*ctx->func)(0,U,a,ctx->funcctx);
363: (*ctx->funcisetbase)(U,ctx->funcctx);
364: VecCopy(U,w);
366: VecGetOwnershipRange(a,&rstart,&rend);
367: VecGetArray(a,&aa);
368: for (i=rstart; i<rend; i++) {
369: VecGetArray(w,&ww);
370: h = ww[i-rstart];
371: if (h == 0.0) h = 1.0;
372: #if !defined(PETSC_USE_COMPLEX)
373: if (h < umin && h >= 0.0) h = umin;
374: else if (h < 0.0 && h > -umin) h = -umin;
375: #else
376: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
377: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
378: #endif
379: h *= epsilon;
380:
381: ww[i-rstart] += h;
382: VecRestoreArray(w,&ww);
383: ierr = (*ctx->funci)(i,w,&v,ctx->funcctx);
384: aa[i-rstart] = (v - aa[i-rstart])/h;
386: /* possibly shift and scale result */
387: aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
389: VecGetArray(w,&ww);
390: ww[i-rstart] -= h;
391: VecRestoreArray(w,&ww);
392: }
393: VecRestoreArray(a,&aa);
394: return(0);
395: }
397: #undef __FUNCT__
399: int MatShift_MFFD(PetscScalar *a,Mat Y)
400: {
401: MatSNESMFCtx shell = (MatSNESMFCtx)Y->data;
403: shell->vshift += *a;
404: return(0);
405: }
407: #undef __FUNCT__
409: int MatScale_MFFD(PetscScalar *a,Mat Y)
410: {
411: MatSNESMFCtx shell = (MatSNESMFCtx)Y->data;
413: shell->vscale *= *a;
414: return(0);
415: }
418: #undef __FUNCT__
420: /*@C
421: MatCreateSNESMF - Creates a matrix-free matrix context for use with
422: a SNES solver. This matrix can be used as the Jacobian argument for
423: the routine SNESSetJacobian().
425: Collective on SNES and Vec
427: Input Parameters:
428: + snes - the SNES context
429: - x - vector where SNES solution is to be stored.
431: Output Parameter:
432: . J - the matrix-free matrix
434: Level: advanced
436: Notes:
437: The matrix-free matrix context merely contains the function pointers
438: and work space for performing finite difference approximations of
439: Jacobian-vector products, F'(u)*a,
441: The default code uses the following approach to compute h
443: .vb
444: F'(u)*a = [F(u+h*a) - F(u)]/h where
445: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
446: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
447: where
448: error_rel = square root of relative error in function evaluation
449: umin = minimum iterate parameter
450: .ve
452: The user can set the error_rel via MatSNESMFSetFunctionError() and
453: umin via MatSNESMFDefaultSetUmin(); see the nonlinear solvers chapter
454: of the users manual for details.
456: The user should call MatDestroy() when finished with the matrix-free
457: matrix context.
459: Options Database Keys:
460: + -snes_mf_err <error_rel> - Sets error_rel
461: . -snes_mf_unim <umin> - Sets umin (for default PETSc routine that computes h only)
462: - -snes_mf_ksp_monitor - KSP monitor routine that prints differencing h
464: .keywords: SNES, default, matrix-free, create, matrix
466: .seealso: MatDestroy(), MatSNESMFSetFunctionError(), MatSNESMFDefaultSetUmin()
467: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(), MatCreateMF(),
468: MatSNESMFGetH(),MatSNESMFKSPMonitor(), MatSNESMFRegisterDynamic), MatSNESMFComputeJacobian()
469:
470: @*/
471: int MatCreateSNESMF(SNES snes,Vec x,Mat *J)
472: {
473: MatSNESMFCtx mfctx;
474: int ierr;
477: MatCreateMF(x,J);
479: mfctx = (MatSNESMFCtx)(*J)->data;
480: mfctx->snes = snes;
481: mfctx->usesnes = PETSC_TRUE;
482: PetscLogObjectParent(snes,*J);
483: return(0);
484: }
486: EXTERN_C_BEGIN
487: #undef __FUNCT__
489: int MatSNESMFSetBase_FD(Mat J,Vec U)
490: {
491: int ierr;
492: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
495: MatSNESMFResetHHistory(J);
496: ctx->current_u = U;
497: ctx->usesnes = PETSC_FALSE;
498: if (ctx->w == PETSC_NULL) {
499: VecDuplicate(ctx->current_u, &ctx->w);
500: }
501: return(0);
502: }
503: EXTERN_C_END
505: EXTERN_C_BEGIN
506: #undef __FUNCT__
508: int MatSNESMFSetCheckh_FD(Mat J,int (*fun)(Vec,Vec,PetscScalar*,void*),void*ectx)
509: {
510: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
513: ctx->checkh = fun;
514: ctx->checkhctx = ectx;
515: return(0);
516: }
517: EXTERN_C_END
519: #undef __FUNCT__
521: /*@
522: MatSNESMFSetFromOptions - Sets the MatSNESMF options from the command line
523: parameter.
525: Collective on Mat
527: Input Parameters:
528: . mat - the matrix obtained with MatCreateSNESMF()
530: Options Database Keys:
531: + -snes_mf_type - <default,wp>
532: - -snes_mf_err - square root of estimated relative error in function evaluation
533: - -snes_mf_period - how often h is recomputed, defaults to 1, everytime
535: Level: advanced
537: .keywords: SNES, matrix-free, parameters
539: .seealso: MatCreateSNESMF(),MatSNESMFSetHHistory(),
540: MatSNESMFResetHHistory(), MatSNESMFKSPMonitor()
541: @*/
542: int MatSNESMFSetFromOptions(Mat mat)
543: {
544: MatSNESMFCtx mfctx = (MatSNESMFCtx)mat->data;
545: int ierr;
546: PetscTruth flg;
547: char ftype[256];
550: if (!MatSNESMFRegisterAllCalled) {MatSNESMFRegisterAll(PETSC_NULL);}
551:
552: PetscOptionsBegin(mfctx->comm,mfctx->prefix,"Set matrix free computation parameters","MatSNESMF");
553: PetscOptionsList("-snes_mf_type","Matrix free type","MatSNESMFSetType",MatSNESMPetscFList,mfctx->type_name,ftype,256,&flg);
554: if (flg) {
555: MatSNESMFSetType(mat,ftype);
556: }
558: PetscOptionsReal("-snes_mf_err","set sqrt relative error in function","MatSNESMFSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
559: PetscOptionsInt("-snes_mf_period","how often h is recomputed","MatSNESMFSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
560: if (mfctx->snes) {
561: PetscOptionsName("-snes_mf_ksp_monitor","Monitor matrix-free parameters","MatSNESMFKSPMonitor",&flg);
562: if (flg) {
563: SLES sles;
564: KSP ksp;
565: SNESGetSLES(mfctx->snes,&sles);
566: SLESGetKSP(sles,&ksp);
567: KSPSetMonitor(ksp,MatSNESMFKSPMonitor,PETSC_NULL,0);
568: }
569: }
570: PetscOptionsName("-snes_mf_check_positivity","Insure that U + h*a is nonnegative","MatSNESMFSetCheckh",&flg);
571: if (flg) {
572: MatSNESMFSetCheckh(mat,MatSNESMFCheckPositivity,0);
573: }
574: if (mfctx->ops->setfromoptions) {
575: (*mfctx->ops->setfromoptions)(mfctx);
576: }
577: PetscOptionsEnd();
578: return(0);
579: }
581: #undef __FUNCT__
583: EXTERN_C_BEGIN
584: int MatCreate_MFFD(Mat A)
585: {
586: MatSNESMFCtx mfctx;
587: int ierr;
590: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
591: SNESInitializePackage(PETSC_NULL);
592: #endif
594: PetscHeaderCreate(mfctx,_p_MatSNESMFCtx,struct _MFOps,MATSNESMFCTX_COOKIE,0,"SNESMF",A->comm,MatDestroy_MFFD,MatView_MFFD);
595: PetscLogObjectCreate(mfctx);
596: mfctx->sp = 0;
597: mfctx->snes = 0;
598: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
599: mfctx->recomputeperiod = 1;
600: mfctx->count = 0;
601: mfctx->currenth = 0.0;
602: mfctx->historyh = PETSC_NULL;
603: mfctx->ncurrenth = 0;
604: mfctx->maxcurrenth = 0;
605: mfctx->type_name = 0;
606: mfctx->usesnes = PETSC_FALSE;
608: mfctx->vshift = 0.0;
609: mfctx->vscale = 1.0;
611: /*
612: Create the empty data structure to contain compute-h routines.
613: These will be filled in below from the command line options or
614: a later call with MatSNESMFSetType() or if that is not called
615: then it will default in the first use of MatMult_MFFD()
616: */
617: mfctx->ops->compute = 0;
618: mfctx->ops->destroy = 0;
619: mfctx->ops->view = 0;
620: mfctx->ops->setfromoptions = 0;
621: mfctx->hctx = 0;
623: mfctx->func = 0;
624: mfctx->funcctx = 0;
625: mfctx->funcvec = 0;
626: mfctx->w = PETSC_NULL;
628: A->data = mfctx;
630: A->ops->mult = MatMult_MFFD;
631: A->ops->destroy = MatDestroy_MFFD;
632: A->ops->view = MatView_MFFD;
633: A->ops->assemblyend = MatAssemblyEnd_MFFD;
634: A->ops->getdiagonal = MatGetDiagonal_MFFD;
635: A->ops->scale = MatScale_MFFD;
636: A->ops->shift = MatShift_MFFD;
637: A->ops->setfromoptions = MatSNESMFSetFromOptions;
639: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetBase_C","MatSNESMFSetBase_FD",MatSNESMFSetBase_FD);
640: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetFunctioniBase_C","MatSNESMFSetFunctioniBase_FD",MatSNESMFSetFunctioniBase_FD);
641: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetFunctioni_C","MatSNESMFSetFunctioni_FD",MatSNESMFSetFunctioni_FD);
642: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetCheckh_C","MatSNESMFSetCheckh_FD",MatSNESMFSetCheckh_FD);
643: mfctx->mat = A;
645: return(0);
646: }
648: EXTERN_C_END
650: #undef __FUNCT__
652: /*@C
653: MatCreateMF - Creates a matrix-free matrix. See also MatCreateSNESMF()
655: Collective on Vec
657: Input Parameters:
658: . x - vector that defines layout of the vectors and matrices
660: Output Parameter:
661: . J - the matrix-free matrix
663: Level: advanced
665: Notes:
666: The matrix-free matrix context merely contains the function pointers
667: and work space for performing finite difference approximations of
668: Jacobian-vector products, F'(u)*a,
670: The default code uses the following approach to compute h
672: .vb
673: F'(u)*a = [F(u+h*a) - F(u)]/h where
674: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
675: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
676: where
677: error_rel = square root of relative error in function evaluation
678: umin = minimum iterate parameter
679: .ve
681: The user can set the error_rel via MatSNESMFSetFunctionError() and
682: umin via MatSNESMFDefaultSetUmin(); see the nonlinear solvers chapter
683: of the users manual for details.
685: The user should call MatDestroy() when finished with the matrix-free
686: matrix context.
688: Options Database Keys:
689: + -snes_mf_err <error_rel> - Sets error_rel
690: . -snes_mf_unim <umin> - Sets umin (for default PETSc routine that computes h only)
691: . -snes_mf_ksp_monitor - KSP monitor routine that prints differencing h
692: - -snes_mf_check_positivity
694: .keywords: default, matrix-free, create, matrix
696: .seealso: MatDestroy(), MatSNESMFSetFunctionError(), MatSNESMFDefaultSetUmin()
697: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(), MatCreateSNESMF(),
698: MatSNESMFGetH(),MatSNESMFKSPMonitor(), MatSNESMFRegisterDynamic),, MatSNESMFComputeJacobian()
699:
700: @*/
701: int MatCreateMF(Vec x,Mat *J)
702: {
703: MPI_Comm comm;
704: int n,nloc,ierr;
707: PetscObjectGetComm((PetscObject)x,&comm);
708: VecGetSize(x,&n);
709: VecGetLocalSize(x,&nloc);
710: MatCreate(comm,nloc,nloc,n,n,J);
711: MatRegister(MATMFFD,0,"MatCreate_MFFD",MatCreate_MFFD);
712: MatSetType(*J,MATMFFD);
713: return(0);
714: }
717: #undef __FUNCT__
719: /*@
720: MatSNESMFGetH - Gets the last value that was used as the differencing
721: parameter.
723: Not Collective
725: Input Parameters:
726: . mat - the matrix obtained with MatCreateSNESMF()
728: Output Paramter:
729: . h - the differencing step size
731: Level: advanced
733: .keywords: SNES, matrix-free, parameters
735: .seealso: MatCreateSNESMF(),MatSNESMFSetHHistory(),
736: MatSNESMFResetHHistory(),MatSNESMFKSPMonitor()
737: @*/
738: int MatSNESMFGetH(Mat mat,PetscScalar *h)
739: {
740: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
743: *h = ctx->currenth;
744: return(0);
745: }
747: #undef __FUNCT__
749: /*
750: MatSNESMFKSPMonitor - A KSP monitor for use with the default PETSc
751: SNES matrix free routines. Prints the differencing parameter used at
752: each step.
753: */
754: int MatSNESMFKSPMonitor(KSP ksp,int n,PetscReal rnorm,void *dummy)
755: {
756: PC pc;
757: MatSNESMFCtx ctx;
758: int ierr;
759: Mat mat;
760: MPI_Comm comm;
761: PetscTruth nonzeroinitialguess;
764: PetscObjectGetComm((PetscObject)ksp,&comm);
765: KSPGetPC(ksp,&pc);
766: KSPGetInitialGuessNonzero(ksp,&nonzeroinitialguess);
767: PCGetOperators(pc,&mat,PETSC_NULL,PETSC_NULL);
768: ctx = (MatSNESMFCtx)mat->data;
770: if (n > 0 || nonzeroinitialguess) {
771: #if defined(PETSC_USE_COMPLEX)
772: PetscPrintf(comm,"%d KSP Residual norm %14.12e h %g + %g in",n,rnorm,
773: PetscRealPart(ctx->currenth),PetscImaginaryPart(ctx->currenth));
774: #else
775: PetscPrintf(comm,"%d KSP Residual norm %14.12e h %g n",n,rnorm,ctx->currenth);
776: #endif
777: } else {
778: PetscPrintf(comm,"%d KSP Residual norm %14.12en",n,rnorm);
779: }
780: return(0);
781: }
783: #undef __FUNCT__
785: /*@C
786: MatSNESMFSetFunction - Sets the function used in applying the matrix free.
788: Collective on Mat
790: Input Parameters:
791: + mat - the matrix free matrix created via MatCreateSNESMF()
792: . v - workspace vector
793: . func - the function to use
794: - funcctx - optional function context passed to function
796: Level: advanced
798: Notes:
799: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
800: matrix inside your compute Jacobian routine
802: If this is not set then it will use the function set with SNESSetFunction()
804: .keywords: SNES, matrix-free, function
806: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
807: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
808: MatSNESMFKSPMonitor(), SNESetFunction()
809: @*/
810: int MatSNESMFSetFunction(Mat mat,Vec v,int (*func)(SNES,Vec,Vec,void *),void *funcctx)
811: {
812: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
815: ctx->func = func;
816: ctx->funcctx = funcctx;
817: ctx->funcvec = v;
818: return(0);
819: }
821: #undef __FUNCT__
823: /*@C
824: MatSNESMFSetFunctioni - Sets the function for a single component
826: Collective on Mat
828: Input Parameters:
829: + mat - the matrix free matrix created via MatCreateSNESMF()
830: - funci - the function to use
832: Level: advanced
834: Notes:
835: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
836: matrix inside your compute Jacobian routine
839: .keywords: SNES, matrix-free, function
841: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
842: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
843: MatSNESMFKSPMonitor(), SNESetFunction()
844: @*/
845: int MatSNESMFSetFunctioni(Mat mat,int (*funci)(int,Vec,PetscScalar*,void *))
846: {
847: int ierr,(*f)(Mat,int (*)(int,Vec,PetscScalar*,void *));
851: PetscObjectQueryFunction((PetscObject)mat,"MatSNESMFSetFunctioni_C",(void (**)(void))&f);
852: if (f) {
853: (*f)(mat,funci);
854: }
855: return(0);
856: }
859: #undef __FUNCT__
861: /*@C
862: MatSNESMFSetFunctioniBase - Sets the base vector for a single component function evaluation
864: Collective on Mat
866: Input Parameters:
867: + mat - the matrix free matrix created via MatCreateSNESMF()
868: - func - the function to use
870: Level: advanced
872: Notes:
873: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
874: matrix inside your compute Jacobian routine
877: .keywords: SNES, matrix-free, function
879: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
880: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
881: MatSNESMFKSPMonitor(), SNESetFunction()
882: @*/
883: int MatSNESMFSetFunctioniBase(Mat mat,int (*func)(Vec,void *))
884: {
885: int ierr,(*f)(Mat,int (*)(Vec,void *));
889: PetscObjectQueryFunction((PetscObject)mat,"MatSNESMFSetFunctioniBase_C",(void (**)(void))&f);
890: if (f) {
891: (*f)(mat,func);
892: }
893: return(0);
894: }
897: #undef __FUNCT__
899: /*@
900: MatSNESMFSetPeriod - Sets how often h is recomputed, by default it is everytime
902: Collective on Mat
904: Input Parameters:
905: + mat - the matrix free matrix created via MatCreateSNESMF()
906: - period - 1 for everytime, 2 for every second etc
908: Options Database Keys:
909: + -snes_mf_period <period>
911: Level: advanced
914: .keywords: SNES, matrix-free, parameters
916: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
917: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
918: MatSNESMFKSPMonitor()
919: @*/
920: int MatSNESMFSetPeriod(Mat mat,int period)
921: {
922: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
925: ctx->recomputeperiod = period;
926: return(0);
927: }
929: #undef __FUNCT__
931: /*@
932: MatSNESMFSetFunctionError - Sets the error_rel for the approximation of
933: matrix-vector products using finite differences.
935: Collective on Mat
937: Input Parameters:
938: + mat - the matrix free matrix created via MatCreateSNESMF()
939: - error_rel - relative error (should be set to the square root of
940: the relative error in the function evaluations)
942: Options Database Keys:
943: + -snes_mf_err <error_rel> - Sets error_rel
945: Level: advanced
947: Notes:
948: The default matrix-free matrix-vector product routine computes
949: .vb
950: F'(u)*a = [F(u+h*a) - F(u)]/h where
951: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
952: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
953: .ve
955: .keywords: SNES, matrix-free, parameters
957: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
958: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
959: MatSNESMFKSPMonitor()
960: @*/
961: int MatSNESMFSetFunctionError(Mat mat,PetscReal error)
962: {
963: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
966: if (error != PETSC_DEFAULT) ctx->error_rel = error;
967: return(0);
968: }
970: #undef __FUNCT__
972: /*@
973: MatSNESMFAddNullSpace - Provides a null space that an operator is
974: supposed to have. Since roundoff will create a small component in
975: the null space, if you know the null space you may have it
976: automatically removed.
978: Collective on Mat
980: Input Parameters:
981: + J - the matrix-free matrix context
982: - nullsp - object created with MatNullSpaceCreate()
984: Level: advanced
986: .keywords: SNES, matrix-free, null space
988: .seealso: MatNullSpaceCreate(), MatSNESMFGetH(), MatCreateSNESMF(),
989: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
990: MatSNESMFKSPMonitor(), MatSNESMFErrorRel()
991: @*/
992: int MatSNESMFAddNullSpace(Mat J,MatNullSpace nullsp)
993: {
994: int ierr;
995: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
996: MPI_Comm comm;
999: PetscObjectGetComm((PetscObject)J,&comm);
1001: ctx->sp = nullsp;
1002: ierr = PetscObjectReference((PetscObject)nullsp);
1003: return(0);
1004: }
1006: #undef __FUNCT__
1008: /*@
1009: MatSNESMFSetHHistory - Sets an array to collect a history of the
1010: differencing values (h) computed for the matrix-free product.
1012: Collective on Mat
1014: Input Parameters:
1015: + J - the matrix-free matrix context
1016: . histroy - space to hold the history
1017: - nhistory - number of entries in history, if more entries are generated than
1018: nhistory, then the later ones are discarded
1020: Level: advanced
1022: Notes:
1023: Use MatSNESMFResetHHistory() to reset the history counter and collect
1024: a new batch of differencing parameters, h.
1026: .keywords: SNES, matrix-free, h history, differencing history
1028: .seealso: MatSNESMFGetH(), MatCreateSNESMF(),
1029: MatSNESMFResetHHistory(),
1030: MatSNESMFKSPMonitor(), MatSNESMFSetFunctionError()
1032: @*/
1033: int MatSNESMFSetHHistory(Mat J,PetscScalar *history,int nhistory)
1034: {
1035: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
1038: ctx->historyh = history;
1039: ctx->maxcurrenth = nhistory;
1040: ctx->currenth = 0;
1041: return(0);
1042: }
1044: #undef __FUNCT__
1046: /*@
1047: MatSNESMFResetHHistory - Resets the counter to zero to begin
1048: collecting a new set of differencing histories.
1050: Collective on Mat
1052: Input Parameters:
1053: . J - the matrix-free matrix context
1055: Level: advanced
1057: Notes:
1058: Use MatSNESMFSetHHistory() to create the original history counter.
1060: .keywords: SNES, matrix-free, h history, differencing history
1062: .seealso: MatSNESMFGetH(), MatCreateSNESMF(),
1063: MatSNESMFSetHHistory(),
1064: MatSNESMFKSPMonitor(), MatSNESMFSetFunctionError()
1066: @*/
1067: int MatSNESMFResetHHistory(Mat J)
1068: {
1069: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
1072: ctx->ncurrenth = 0;
1073: return(0);
1074: }
1076: #undef __FUNCT__
1078: int MatSNESMFComputeJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
1079: {
1082: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
1083: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
1084: return(0);
1085: }
1087: #undef __FUNCT__
1089: /*@
1090: MatSNESMFSetBase - Sets the vector U at which matrix vector products of the
1091: Jacobian are computed
1093: Collective on Mat
1095: Input Parameters:
1096: + J - the MatSNESMF matrix
1097: - U - the vector
1099: Notes: This is rarely used directly
1101: Level: advanced
1103: @*/
1104: int MatSNESMFSetBase(Mat J,Vec U)
1105: {
1106: int ierr,(*f)(Mat,Vec);
1111: PetscObjectQueryFunction((PetscObject)J,"MatSNESMFSetBase_C",(void (**)(void))&f);
1112: if (f) {
1113: (*f)(J,U);
1114: }
1115: return(0);
1116: }
1118: #undef __FUNCT__
1120: /*@C
1121: MatSNESMFSetCheckh - Sets a function that checks the computed h and adjusts
1122: it to satisfy some criteria
1124: Collective on Mat
1126: Input Parameters:
1127: + J - the MatSNESMF matrix
1128: . fun - the function that checks h
1129: - ctx - any context needed by the function
1131: Options Database Keys:
1132: . -snes_mf_check_positivity
1134: Level: advanced
1136: Notes: For example, MatSNESMFSetCheckPositivity() insures that all entries
1137: of U + h*a are non-negative
1139: .seealso: MatSNESMFSetCheckPositivity()
1140: @*/
1141: int MatSNESMFSetCheckh(Mat J,int (*fun)(Vec,Vec,PetscScalar*,void*),void* ctx)
1142: {
1143: int ierr,(*f)(Mat,int (*)(Vec,Vec,PetscScalar*,void*),void*);
1147: PetscObjectQueryFunction((PetscObject)J,"MatSNESMFSetCheckh_C",(void (**)(void))&f);
1148: if (f) {
1149: (*f)(J,fun,ctx);
1150: }
1151: return(0);
1152: }
1154: #undef __FUNCT__
1156: /*@
1157: MatSNESMFCheckPositivity - Checks that all entries in U + h*a are positive or
1158: zero, decreases h until this is satisfied.
1160: Collective on Vec
1162: Input Parameters:
1163: + U - base vector that is added to
1164: . a - vector that is added
1165: . h - scaling factor on a
1166: - dummy - context variable (unused)
1168: Options Database Keys:
1169: . -snes_mf_check_positivity
1171: Level: advanced
1173: Notes: This is rarely used directly, rather it is passed as an argument to
1174: MatSNESMFSetCheckh()
1176: .seealso: MatSNESMFSetCheckh()
1177: @*/
1178: int MatSNESMFCheckPositivity(Vec U,Vec a,PetscScalar *h,void *dummy)
1179: {
1180: PetscReal val, minval;
1181: PetscScalar *u_vec, *a_vec;
1182: int ierr, i, size;
1183: MPI_Comm comm;
1186: PetscObjectGetComm((PetscObject)U,&comm);
1187: VecGetArray(U,&u_vec);
1188: VecGetArray(a,&a_vec);
1189: VecGetLocalSize(U,&size);
1190: minval = PetscAbsScalar(*h*1.01);
1191: for(i=0;i<size;i++) {
1192: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1193: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1194: if (val < minval) minval = val;
1195: }
1196: }
1197: VecRestoreArray(U,&u_vec);
1198: VecRestoreArray(a,&a_vec);
1199: PetscGlobalMin(&minval,&val,comm);
1200: if (val <= PetscAbsScalar(*h)) {
1201: PetscLogInfo(U,"MatSNESMFCheckPositivity: Scaling back h from %g to %gn",PetscRealPart(*h),.99*val);
1202: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1203: else *h = -0.99*val;
1204: }
1205: return(0);
1206: }