Actual source code: snesmfj.c
1: #define PETSCSNES_DLL
3: #include include/private/matimpl.h
4: #include src/snes/mf/snesmfj.h
6: PetscFList MatSNESMPetscFList = 0;
7: PetscTruth MatSNESMFRegisterAllCalled = PETSC_FALSE;
9: PetscCookie MATSNESMFCTX_COOKIE = 0;
10: PetscEvent MATSNESMF_Mult = 0;
14: /*@C
15: MatSNESMFSetType - Sets the method that is used to compute the
16: differencing parameter for finite differene matrix-free formulations.
18: Input Parameters:
19: + mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMF()
20: or MatSetType(mat,MATMFFD);
21: - ftype - the type requested, either MATSNESMF_WP or MATSNESMF_DS
23: Level: advanced
25: Notes:
26: For example, such routines can compute h for use in
27: Jacobian-vector products of the form
29: F(x+ha) - F(x)
30: F'(u)a ~= ----------------
31: h
33: .seealso: MatCreateSNESMF(), MatSNESMFRegisterDynamic)
34: @*/
35: PetscErrorCode MatSNESMFSetType(Mat mat,const MatSNESMFType ftype)
36: {
37: PetscErrorCode ierr,(*r)(MatSNESMFCtx);
38: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
39: PetscTruth match;
40:
45: /* already set, so just return */
46: PetscTypeCompare((PetscObject)ctx,ftype,&match);
47: if (match) return(0);
49: /* destroy the old one if it exists */
50: if (ctx->ops->destroy) {
51: (*ctx->ops->destroy)(ctx);
52: }
54: /* Get the function pointers for the requrested method */
55: if (!MatSNESMFRegisterAllCalled) {MatSNESMFRegisterAll(PETSC_NULL);}
56: PetscFListFind(ctx->comm,MatSNESMPetscFList,ftype,(void (**)(void)) &r);
57: if (!r) SETERRQ1(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatSNESMF type %s given",ftype);
58: (*r)(ctx);
59: PetscObjectChangeTypeName((PetscObject)ctx,ftype);
60: return(0);
61: }
67: PetscErrorCode MatSNESMFSetFunctioniBase_FD(Mat mat,FCN1 func)
68: {
69: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
72: ctx->funcisetbase = func;
73: return(0);
74: }
81: PetscErrorCode MatSNESMFSetFunctioni_FD(Mat mat,FCN2 funci)
82: {
83: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
86: ctx->funci = funci;
87: return(0);
88: }
94: PetscErrorCode MatSNESMFRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(MatSNESMFCtx))
95: {
97: char fullname[PETSC_MAX_PATH_LEN];
100: PetscFListConcat(path,name,fullname);
101: PetscFListAdd(&MatSNESMPetscFList,sname,fullname,(void (*)(void))function);
102: return(0);
103: }
108: /*@C
109: MatSNESMFRegisterDestroy - Frees the list of MatSNESMF methods that were
110: registered by MatSNESMFRegisterDynamic).
112: Not Collective
114: Level: developer
116: .keywords: MatSNESMF, register, destroy
118: .seealso: MatSNESMFRegisterDynamic), MatSNESMFRegisterAll()
119: @*/
120: PetscErrorCode MatSNESMFRegisterDestroy(void)
121: {
125: if (MatSNESMPetscFList) {
126: PetscFListDestroy(&MatSNESMPetscFList);
127: MatSNESMPetscFList = 0;
128: }
129: MatSNESMFRegisterAllCalled = PETSC_FALSE;
130: return(0);
131: }
133: /* ----------------------------------------------------------------------------------------*/
136: PetscErrorCode MatDestroy_MFFD(Mat mat)
137: {
139: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
142: if (ctx->w) {
143: VecDestroy(ctx->w);
144: }
145: if (ctx->ops->destroy) {(*ctx->ops->destroy)(ctx);}
146: if (ctx->sp) {MatNullSpaceDestroy(ctx->sp);}
147: PetscHeaderDestroy(ctx);
149: PetscObjectComposeFunction((PetscObject)mat,"MatSNESMFSetBase_C","",PETSC_NULL);
150: PetscObjectComposeFunction((PetscObject)mat,"MatSNESMFSetFunctioniBase_C","",PETSC_NULL);
151: PetscObjectComposeFunction((PetscObject)mat,"MatSNESMFSetFunctioni_C","",PETSC_NULL);
152: PetscObjectComposeFunction((PetscObject)mat,"MatSNESMFSetCheckh_C","",PETSC_NULL);
154: return(0);
155: }
159: /*
160: MatSNESMFView_MFFD - Views matrix-free parameters.
162: */
163: PetscErrorCode MatView_MFFD(Mat J,PetscViewer viewer)
164: {
166: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
167: PetscTruth iascii;
170: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
171: if (iascii) {
172: PetscViewerASCIIPrintf(viewer," SNES matrix-free approximation:\n");
173: PetscViewerASCIIPrintf(viewer," err=%G (relative error in function evaluation)\n",ctx->error_rel);
174: if (!ctx->type_name) {
175: PetscViewerASCIIPrintf(viewer," The compute h routine has not yet been set\n");
176: } else {
177: PetscViewerASCIIPrintf(viewer," Using %s compute h routine\n",ctx->type_name);
178: }
179: if (ctx->ops->view) {
180: (*ctx->ops->view)(ctx,viewer);
181: }
182: } else {
183: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES matrix free matrix",((PetscObject)viewer)->type_name);
184: }
185: return(0);
186: }
190: /*
191: MatAssemblyEnd_MFFD - Resets the ctx->ncurrenth to zero. This
192: allows the user to indicate the beginning of a new linear solve by calling
193: MatAssemblyXXX() on the matrix free matrix. This then allows the
194: MatSNESMFCreate_WP() to properly compute ||U|| only the first time
195: in the linear solver rather than every time.
196: */
197: PetscErrorCode MatAssemblyEnd_MFFD(Mat J,MatAssemblyType mt)
198: {
200: MatSNESMFCtx j = (MatSNESMFCtx)J->data;
203: MatSNESMFResetHHistory(J);
204: if (j->usesnes) {
205: SNESGetSolution(j->snes,&j->current_u);
206: SNESGetFunction(j->snes,&j->current_f,PETSC_NULL,PETSC_NULL);
207: if (!j->w) {
208: VecDuplicate(j->current_u, &j->w);
209: }
210: }
211: j->vshift = 0.0;
212: j->vscale = 1.0;
213: return(0);
214: }
218: /*
219: MatMult_MFFD - Default matrix-free form for Jacobian-vector product, y = F'(u)*a:
221: y ~= (F(u + ha) - F(u))/h,
222: where F = nonlinear function, as set by SNESSetFunction()
223: u = current iterate
224: h = difference interval
225: */
226: PetscErrorCode MatMult_MFFD(Mat mat,Vec a,Vec y)
227: {
228: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
229: SNES snes;
230: PetscScalar h;
231: Vec w,U,F;
232: PetscErrorCode ierr,(*eval_fct)(SNES,Vec,Vec)=0;
233: PetscTruth zeroa;
236: /* We log matrix-free matrix-vector products separately, so that we can
237: separate the performance monitoring from the cases that use conventional
238: storage. We may eventually modify event logging to associate events
239: with particular objects, hence alleviating the more general problem. */
242: snes = ctx->snes;
243: w = ctx->w;
244: U = ctx->current_u;
246: /*
247: Compute differencing parameter
248: */
249: if (!ctx->ops->compute) {
250: MatSNESMFSetType(mat,MATSNESMF_WP);
251: MatSNESMFSetFromOptions(mat);
252: }
253: (*ctx->ops->compute)(ctx,U,a,&h,&zeroa);
254: if (zeroa) {
255: VecSet(y,0.0);
256: return(0);
257: }
259: if (ctx->checkh) {
260: (*ctx->checkh)(U,a,&h,ctx->checkhctx);
261: }
263: /* keep a record of the current differencing parameter h */
264: ctx->currenth = h;
265: #if defined(PETSC_USE_COMPLEX)
266: PetscInfo2(mat,"Current differencing parameter: %G + %G i\n",PetscRealPart(h),PetscImaginaryPart(h));
267: #else
268: PetscInfo1(mat,"Current differencing parameter: %15.12e\n",h);
269: #endif
270: if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
271: ctx->historyh[ctx->ncurrenth] = h;
272: }
273: ctx->ncurrenth++;
275: /* w = u + ha */
276: VecWAXPY(w,h,a,U);
278: if (ctx->usesnes) {
279: eval_fct = SNESComputeFunction;
280: F = ctx->current_f;
281: if (!F) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"You must call MatAssembly() even on matrix-free matrices");
282: (*eval_fct)(snes,w,y);
283: } else {
284: F = ctx->funcvec;
285: /* compute func(U) as base for differencing */
286: if (ctx->ncurrenth == 1) {
287: (*ctx->func)(snes,U,F,ctx->funcctx);
288: }
289: (*ctx->func)(snes,w,y,ctx->funcctx);
290: }
292: VecAXPY(y,-1.0,F);
293: VecScale(y,1.0/h);
295: VecAXPBY(y,ctx->vshift,ctx->vscale,a);
297: if (ctx->sp) {MatNullSpaceRemove(ctx->sp,y,PETSC_NULL);}
300: return(0);
301: }
305: /*
306: MatGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix
308: y ~= (F(u + ha) - F(u))/h,
309: where F = nonlinear function, as set by SNESSetFunction()
310: u = current iterate
311: h = difference interval
312: */
313: PetscErrorCode MatGetDiagonal_MFFD(Mat mat,Vec a)
314: {
315: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
316: PetscScalar h,*aa,*ww,v;
317: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
318: Vec w,U;
320: PetscInt i,rstart,rend;
323: if (!ctx->funci) {
324: SETERRQ(PETSC_ERR_ORDER,"Requires calling MatSNESMFSetFunctioni() first");
325: }
327: w = ctx->w;
328: U = ctx->current_u;
329: (*ctx->func)(0,U,a,ctx->funcctx);
330: (*ctx->funcisetbase)(U,ctx->funcctx);
331: VecCopy(U,w);
333: VecGetOwnershipRange(a,&rstart,&rend);
334: VecGetArray(a,&aa);
335: for (i=rstart; i<rend; i++) {
336: VecGetArray(w,&ww);
337: h = ww[i-rstart];
338: if (h == 0.0) h = 1.0;
339: #if !defined(PETSC_USE_COMPLEX)
340: if (h < umin && h >= 0.0) h = umin;
341: else if (h < 0.0 && h > -umin) h = -umin;
342: #else
343: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
344: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
345: #endif
346: h *= epsilon;
347:
348: ww[i-rstart] += h;
349: VecRestoreArray(w,&ww);
350: (*ctx->funci)(i,w,&v,ctx->funcctx);
351: aa[i-rstart] = (v - aa[i-rstart])/h;
353: /* possibly shift and scale result */
354: aa[i - rstart] = ctx->vshift + ctx->vscale*aa[i-rstart];
356: VecGetArray(w,&ww);
357: ww[i-rstart] -= h;
358: VecRestoreArray(w,&ww);
359: }
360: VecRestoreArray(a,&aa);
361: return(0);
362: }
366: PetscErrorCode MatShift_MFFD(Mat Y,PetscScalar a)
367: {
368: MatSNESMFCtx shell = (MatSNESMFCtx)Y->data;
370: shell->vshift += a;
371: return(0);
372: }
376: PetscErrorCode MatScale_MFFD(Mat Y,PetscScalar a)
377: {
378: MatSNESMFCtx shell = (MatSNESMFCtx)Y->data;
380: shell->vscale *= a;
381: return(0);
382: }
387: /*@
388: MatCreateSNESMF - Creates a matrix-free matrix context for use with
389: a SNES solver. This matrix can be used as the Jacobian argument for
390: the routine SNESSetJacobian().
392: Collective on SNES and Vec
394: Input Parameters:
395: + snes - the SNES context
396: - x - vector where SNES solution is to be stored.
398: Output Parameter:
399: . J - the matrix-free matrix
401: Level: advanced
403: Notes:
404: The matrix-free matrix context merely contains the function pointers
405: and work space for performing finite difference approximations of
406: Jacobian-vector products, F'(u)*a,
408: The default code uses the following approach to compute h
410: .vb
411: F'(u)*a = [F(u+h*a) - F(u)]/h where
412: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
413: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
414: where
415: error_rel = square root of relative error in function evaluation
416: umin = minimum iterate parameter
417: .ve
418: (see MATSNESMF_WP or MATSNESMF_DS)
419:
420: The user can set the error_rel via MatSNESMFSetFunctionError() and
421: umin via MatSNESMFDefaultSetUmin(); see the nonlinear solvers chapter
422: of the users manual for details.
424: The user should call MatDestroy() when finished with the matrix-free
425: matrix context.
427: Options Database Keys:
428: + -snes_mf_err <error_rel> - Sets error_rel
429: + -snes_mf_type - wp or ds (see MATSNESMF_WP or MATSNESMF_DS)
430: . -snes_mf_unim <umin> - Sets umin (for default PETSc routine that computes h only)
431: - -snes_mf_ksp_monitor - KSP monitor routine that prints differencing h
433: .keywords: SNES, default, matrix-free, create, matrix
435: .seealso: MatDestroy(), MatSNESMFSetFunctionError(), MatSNESMFDefaultSetUmin()
436: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(), MatCreateMF(),
437: MatSNESMFGetH(),MatSNESMFKSPMonitor(), MatSNESMFRegisterDynamic), MatSNESMFComputeJacobian()
438:
439: @*/
440: PetscErrorCode MatCreateSNESMF(SNES snes,Vec x,Mat *J)
441: {
442: MatSNESMFCtx mfctx;
446: MatCreateMF(x,J);
448: mfctx = (MatSNESMFCtx)(*J)->data;
449: mfctx->snes = snes;
450: mfctx->usesnes = PETSC_TRUE;
451: PetscLogObjectParent(snes,*J);
452: return(0);
453: }
458: PetscErrorCode MatSNESMFSetBase_FD(Mat J,Vec U)
459: {
461: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
464: MatSNESMFResetHHistory(J);
465: ctx->current_u = U;
466: ctx->usesnes = PETSC_FALSE;
467: if (!ctx->w) {
468: VecDuplicate(ctx->current_u, &ctx->w);
469: }
470: J->assembled = PETSC_TRUE;
471: return(0);
472: }
478: PetscErrorCode MatSNESMFSetCheckh_FD(Mat J,FCN3 fun,void*ectx)
479: {
480: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
483: ctx->checkh = fun;
484: ctx->checkhctx = ectx;
485: return(0);
486: }
491: /*@
492: MatSNESMFSetFromOptions - Sets the MatSNESMF options from the command line
493: parameter.
495: Collective on Mat
497: Input Parameters:
498: . mat - the matrix obtained with MatCreateSNESMF()
500: Options Database Keys:
501: + -snes_mf_type - wp or ds (see MATSNESMF_WP or MATSNESMF_DS)
502: - -snes_mf_err - square root of estimated relative error in function evaluation
503: - -snes_mf_period - how often h is recomputed, defaults to 1, everytime
505: Level: advanced
507: .keywords: SNES, matrix-free, parameters
509: .seealso: MatCreateSNESMF(),MatSNESMFSetHHistory(),
510: MatSNESMFResetHHistory(), MatSNESMFKSPMonitor()
511: @*/
512: PetscErrorCode MatSNESMFSetFromOptions(Mat mat)
513: {
514: MatSNESMFCtx mfctx = (MatSNESMFCtx)mat->data;
516: PetscTruth flg;
517: char ftype[256];
520: if (!MatSNESMFRegisterAllCalled) {MatSNESMFRegisterAll(PETSC_NULL);}
521:
522: PetscOptionsBegin(mfctx->comm,mfctx->prefix,"Set matrix free computation parameters","MatSNESMF");
523: PetscOptionsList("-snes_mf_type","Matrix free type","MatSNESMFSetType",MatSNESMPetscFList,mfctx->type_name,ftype,256,&flg);
524: if (flg) {
525: MatSNESMFSetType(mat,ftype);
526: }
528: PetscOptionsReal("-snes_mf_err","set sqrt relative error in function","MatSNESMFSetFunctionError",mfctx->error_rel,&mfctx->error_rel,0);
529: PetscOptionsInt("-snes_mf_period","how often h is recomputed","MatSNESMFSetPeriod",mfctx->recomputeperiod,&mfctx->recomputeperiod,0);
530: if (mfctx->snes) {
531: PetscOptionsName("-snes_mf_ksp_monitor","Monitor matrix-free parameters","MatSNESMFKSPMonitor",&flg);
532: if (flg) {
533: KSP ksp;
534: SNESGetKSP(mfctx->snes,&ksp);
535: KSPMonitorSet(ksp,MatSNESMFKSPMonitor,PETSC_NULL,0);
536: }
537: }
538: PetscOptionsName("-snes_mf_check_positivity","Insure that U + h*a is nonnegative","MatSNESMFSetCheckh",&flg);
539: if (flg) {
540: MatSNESMFSetCheckh(mat,MatSNESMFCheckPositivity,0);
541: }
542: if (mfctx->ops->setfromoptions) {
543: (*mfctx->ops->setfromoptions)(mfctx);
544: }
545: PetscOptionsEnd();
546: return(0);
547: }
549: /*MC
550: MATMFFD - MATMFFD = "mffd" - A matrix free matrix type.
552: Level: advanced
554: .seealso: MatCreateMF(), MatCreateSNESMF()
555: M*/
559: PetscErrorCode MatCreate_MFFD(Mat A)
560: {
561: MatSNESMFCtx mfctx;
565: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
566: SNESInitializePackage(PETSC_NULL);
567: #endif
569: PetscHeaderCreate(mfctx,_p_MatSNESMFCtx,struct _MFOps,MATSNESMFCTX_COOKIE,0,"SNESMF",A->comm,MatDestroy_MFFD,MatView_MFFD);
570: mfctx->sp = 0;
571: mfctx->snes = 0;
572: mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON;
573: mfctx->recomputeperiod = 1;
574: mfctx->count = 0;
575: mfctx->currenth = 0.0;
576: mfctx->historyh = PETSC_NULL;
577: mfctx->ncurrenth = 0;
578: mfctx->maxcurrenth = 0;
579: mfctx->type_name = 0;
580: mfctx->usesnes = PETSC_FALSE;
582: mfctx->vshift = 0.0;
583: mfctx->vscale = 1.0;
585: /*
586: Create the empty data structure to contain compute-h routines.
587: These will be filled in below from the command line options or
588: a later call with MatSNESMFSetType() or if that is not called
589: then it will default in the first use of MatMult_MFFD()
590: */
591: mfctx->ops->compute = 0;
592: mfctx->ops->destroy = 0;
593: mfctx->ops->view = 0;
594: mfctx->ops->setfromoptions = 0;
595: mfctx->hctx = 0;
597: mfctx->func = 0;
598: mfctx->funcctx = 0;
599: mfctx->funcvec = 0;
600: mfctx->w = PETSC_NULL;
602: A->data = mfctx;
604: A->ops->mult = MatMult_MFFD;
605: A->ops->destroy = MatDestroy_MFFD;
606: A->ops->view = MatView_MFFD;
607: A->ops->assemblyend = MatAssemblyEnd_MFFD;
608: A->ops->getdiagonal = MatGetDiagonal_MFFD;
609: A->ops->scale = MatScale_MFFD;
610: A->ops->shift = MatShift_MFFD;
611: A->ops->setfromoptions = MatSNESMFSetFromOptions;
612: A->assembled = PETSC_TRUE;
614: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetBase_C","MatSNESMFSetBase_FD",MatSNESMFSetBase_FD);
615: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetFunctioniBase_C","MatSNESMFSetFunctioniBase_FD",MatSNESMFSetFunctioniBase_FD);
616: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetFunctioni_C","MatSNESMFSetFunctioni_FD",MatSNESMFSetFunctioni_FD);
617: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSNESMFSetCheckh_C","MatSNESMFSetCheckh_FD",MatSNESMFSetCheckh_FD);
618: mfctx->mat = A;
619: PetscObjectChangeTypeName((PetscObject)A,MATMFFD);
620: return(0);
621: }
626: /*@
627: MatCreateMF - Creates a matrix-free matrix. See also MatCreateSNESMF()
629: Collective on Vec
631: Input Parameters:
632: . x - vector that defines layout of the vectors and matrices
634: Output Parameter:
635: . J - the matrix-free matrix
637: Level: advanced
639: Notes:
640: The matrix-free matrix context merely contains the function pointers
641: and work space for performing finite difference approximations of
642: Jacobian-vector products, F'(u)*a,
644: The default code uses the following approach to compute h
646: .vb
647: F'(u)*a = [F(u+h*a) - F(u)]/h where
648: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
649: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise
650: where
651: error_rel = square root of relative error in function evaluation
652: umin = minimum iterate parameter
653: .ve
655: The user can set the error_rel via MatSNESMFSetFunctionError() and
656: umin via MatSNESMFDefaultSetUmin(); see the nonlinear solvers chapter
657: of the users manual for details.
659: The user should call MatDestroy() when finished with the matrix-free
660: matrix context.
662: Options Database Keys:
663: + -snes_mf_err <error_rel> - Sets error_rel
664: . -snes_mf_unim <umin> - Sets umin (for default PETSc routine that computes h only)
665: . -snes_mf_ksp_monitor - KSP monitor routine that prints differencing h
666: - -snes_mf_check_positivity
668: .keywords: default, matrix-free, create, matrix
670: .seealso: MatDestroy(), MatSNESMFSetFunctionError(), MatSNESMFDefaultSetUmin(), MatSNESMFSetFunction()
671: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(), MatCreateSNESMF(),
672: MatSNESMFGetH(),MatSNESMFKSPMonitor(), MatSNESMFRegisterDynamic),, MatSNESMFComputeJacobian()
673:
674: @*/
675: PetscErrorCode MatCreateMF(Vec x,Mat *J)
676: {
677: MPI_Comm comm;
679: PetscInt n,nloc;
682: PetscObjectGetComm((PetscObject)x,&comm);
683: VecGetSize(x,&n);
684: VecGetLocalSize(x,&nloc);
685: MatCreate(comm,J);
686: MatSetSizes(*J,nloc,nloc,n,n);
687: MatRegisterDynamic(MATMFFD,0,"MatCreate_MFFD",MatCreate_MFFD);
688: MatSetType(*J,MATMFFD);
689: return(0);
690: }
695: /*@
696: MatSNESMFGetH - Gets the last value that was used as the differencing
697: parameter.
699: Not Collective
701: Input Parameters:
702: . mat - the matrix obtained with MatCreateSNESMF()
704: Output Paramter:
705: . h - the differencing step size
707: Level: advanced
709: .keywords: SNES, matrix-free, parameters
711: .seealso: MatCreateSNESMF(),MatSNESMFSetHHistory(),
712: MatSNESMFResetHHistory(),MatSNESMFKSPMonitor()
713: @*/
714: PetscErrorCode MatSNESMFGetH(Mat mat,PetscScalar *h)
715: {
716: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
719: *h = ctx->currenth;
720: return(0);
721: }
725: /*@C
726: MatSNESMFKSPMonitor - A KSP monitor for use with the PETSc
727: SNES matrix free routines. Prints the differencing parameter used at
728: each step.
730: Collective on KSP
732: Input Parameters:
733: + ksp - the Krylov solver object
734: . n - the current iteration
735: . rnorm - the current residual norm (may be preconditioned or not depending on solver and options
736: _ dummy - unused argument
738: Options Database:
739: . -snes_mf_ksp_monitor - turn this monitor on
741: Notes: Use KSPMonitorSet(ksp,MatSNESMFKSPMonitor,PETSC_NULL,PETSC_NULL);
743: .seealso: KSPMonitorSet()
745: @*/
746: PetscErrorCode MatSNESMFKSPMonitor(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
747: {
748: MatSNESMFCtx ctx;
750: Mat mat;
751: MPI_Comm comm;
752: PetscTruth nonzeroinitialguess;
755: PetscObjectGetComm((PetscObject)ksp,&comm);
756: KSPGetInitialGuessNonzero(ksp,&nonzeroinitialguess);
757: KSPGetOperators(ksp,&mat,PETSC_NULL,PETSC_NULL);
758: ctx = (MatSNESMFCtx)mat->data;
760: if (n > 0 || nonzeroinitialguess) {
761: #if defined(PETSC_USE_COMPLEX)
762: PetscPrintf(comm,"%D KSP Residual norm %14.12e h %G + %G i\n",n,rnorm,
763: PetscRealPart(ctx->currenth),PetscImaginaryPart(ctx->currenth));
764: #else
765: PetscPrintf(comm,"%D KSP Residual norm %14.12e h %G \n",n,rnorm,ctx->currenth);
766: #endif
767: } else {
768: PetscPrintf(comm,"%D KSP Residual norm %14.12e\n",n,rnorm);
769: }
770: return(0);
771: }
775: /*@C
776: MatSNESMFSetFunction - Sets the function used in applying the matrix free.
778: Collective on Mat
780: Input Parameters:
781: + mat - the matrix free matrix created via MatCreateSNESMF()
782: . v - workspace vector
783: . func - the function to use
784: - funcctx - optional function context passed to function
786: Level: advanced
788: Notes:
789: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
790: matrix inside your compute Jacobian routine
792: If this is not set then it will use the function set with SNESSetFunction()
794: .keywords: SNES, matrix-free, function
796: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
797: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
798: MatSNESMFKSPMonitor(), SNESetFunction()
799: @*/
800: PetscErrorCode MatSNESMFSetFunction(Mat mat,Vec v,PetscErrorCode (*func)(SNES,Vec,Vec,void *),void *funcctx)
801: {
802: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
805: ctx->func = func;
806: ctx->funcctx = funcctx;
807: ctx->funcvec = v;
808: return(0);
809: }
813: /*@C
814: MatSNESMFSetFunctioni - Sets the function for a single component
816: Collective on Mat
818: Input Parameters:
819: + mat - the matrix free matrix created via MatCreateSNESMF()
820: - funci - the function to use
822: Level: advanced
824: Notes:
825: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
826: matrix inside your compute Jacobian routine
829: .keywords: SNES, matrix-free, function
831: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
832: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
833: MatSNESMFKSPMonitor(), SNESetFunction()
834: @*/
835: PetscErrorCode MatSNESMFSetFunctioni(Mat mat,PetscErrorCode (*funci)(PetscInt,Vec,PetscScalar*,void *))
836: {
837: PetscErrorCode ierr,(*f)(Mat,PetscErrorCode (*)(PetscInt,Vec,PetscScalar*,void *));
841: PetscObjectQueryFunction((PetscObject)mat,"MatSNESMFSetFunctioni_C",(void (**)(void))&f);
842: if (f) {
843: (*f)(mat,funci);
844: }
845: return(0);
846: }
851: /*@C
852: MatSNESMFSetFunctioniBase - Sets the base vector for a single component function evaluation
854: Collective on Mat
856: Input Parameters:
857: + mat - the matrix free matrix created via MatCreateSNESMF()
858: - func - the function to use
860: Level: advanced
862: Notes:
863: If you use this you MUST call MatAssemblyBegin()/MatAssemblyEnd() on the matrix free
864: matrix inside your compute Jacobian routine
867: .keywords: SNES, matrix-free, function
869: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
870: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
871: MatSNESMFKSPMonitor(), SNESetFunction()
872: @*/
873: PetscErrorCode MatSNESMFSetFunctioniBase(Mat mat,PetscErrorCode (*func)(Vec,void *))
874: {
875: PetscErrorCode ierr,(*f)(Mat,PetscErrorCode (*)(Vec,void *));
879: PetscObjectQueryFunction((PetscObject)mat,"MatSNESMFSetFunctioniBase_C",(void (**)(void))&f);
880: if (f) {
881: (*f)(mat,func);
882: }
883: return(0);
884: }
889: /*@
890: MatSNESMFSetPeriod - Sets how often h is recomputed, by default it is everytime
892: Collective on Mat
894: Input Parameters:
895: + mat - the matrix free matrix created via MatCreateSNESMF()
896: - period - 1 for everytime, 2 for every second etc
898: Options Database Keys:
899: + -snes_mf_period <period>
901: Level: advanced
904: .keywords: SNES, matrix-free, parameters
906: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
907: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
908: MatSNESMFKSPMonitor()
909: @*/
910: PetscErrorCode MatSNESMFSetPeriod(Mat mat,PetscInt period)
911: {
912: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
915: ctx->recomputeperiod = period;
916: return(0);
917: }
921: /*@
922: MatSNESMFSetFunctionError - Sets the error_rel for the approximation of
923: matrix-vector products using finite differences.
925: Collective on Mat
927: Input Parameters:
928: + mat - the matrix free matrix created via MatCreateSNESMF()
929: - error_rel - relative error (should be set to the square root of
930: the relative error in the function evaluations)
932: Options Database Keys:
933: + -snes_mf_err <error_rel> - Sets error_rel
935: Level: advanced
937: Notes:
938: The default matrix-free matrix-vector product routine computes
939: .vb
940: F'(u)*a = [F(u+h*a) - F(u)]/h where
941: h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1}
942: = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else
943: .ve
945: .keywords: SNES, matrix-free, parameters
947: .seealso: MatCreateSNESMF(),MatSNESMFGetH(),
948: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
949: MatSNESMFKSPMonitor()
950: @*/
951: PetscErrorCode MatSNESMFSetFunctionError(Mat mat,PetscReal error)
952: {
953: MatSNESMFCtx ctx = (MatSNESMFCtx)mat->data;
956: if (error != PETSC_DEFAULT) ctx->error_rel = error;
957: return(0);
958: }
962: /*@
963: MatSNESMFAddNullSpace - Provides a null space that an operator is
964: supposed to have. Since roundoff will create a small component in
965: the null space, if you know the null space you may have it
966: automatically removed.
968: Collective on Mat
970: Input Parameters:
971: + J - the matrix-free matrix context
972: - nullsp - object created with MatNullSpaceCreate()
974: Level: advanced
976: .keywords: SNES, matrix-free, null space
978: .seealso: MatNullSpaceCreate(), MatSNESMFGetH(), MatCreateSNESMF(),
979: MatSNESMFSetHHistory(), MatSNESMFResetHHistory(),
980: MatSNESMFKSPMonitor(), MatSNESMFErrorRel()
981: @*/
982: PetscErrorCode MatSNESMFAddNullSpace(Mat J,MatNullSpace nullsp)
983: {
985: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
986: MPI_Comm comm;
989: PetscObjectGetComm((PetscObject)J,&comm);
991: ctx->sp = nullsp;
992: PetscObjectReference((PetscObject)nullsp);
993: return(0);
994: }
998: /*@
999: MatSNESMFSetHHistory - Sets an array to collect a history of the
1000: differencing values (h) computed for the matrix-free product.
1002: Collective on Mat
1004: Input Parameters:
1005: + J - the matrix-free matrix context
1006: . histroy - space to hold the history
1007: - nhistory - number of entries in history, if more entries are generated than
1008: nhistory, then the later ones are discarded
1010: Level: advanced
1012: Notes:
1013: Use MatSNESMFResetHHistory() to reset the history counter and collect
1014: a new batch of differencing parameters, h.
1016: .keywords: SNES, matrix-free, h history, differencing history
1018: .seealso: MatSNESMFGetH(), MatCreateSNESMF(),
1019: MatSNESMFResetHHistory(),
1020: MatSNESMFKSPMonitor(), MatSNESMFSetFunctionError()
1022: @*/
1023: PetscErrorCode MatSNESMFSetHHistory(Mat J,PetscScalar history[],PetscInt nhistory)
1024: {
1025: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
1028: ctx->historyh = history;
1029: ctx->maxcurrenth = nhistory;
1030: ctx->currenth = 0;
1031: return(0);
1032: }
1036: /*@
1037: MatSNESMFResetHHistory - Resets the counter to zero to begin
1038: collecting a new set of differencing histories.
1040: Collective on Mat
1042: Input Parameters:
1043: . J - the matrix-free matrix context
1045: Level: advanced
1047: Notes:
1048: Use MatSNESMFSetHHistory() to create the original history counter.
1050: .keywords: SNES, matrix-free, h history, differencing history
1052: .seealso: MatSNESMFGetH(), MatCreateSNESMF(),
1053: MatSNESMFSetHHistory(),
1054: MatSNESMFKSPMonitor(), MatSNESMFSetFunctionError()
1056: @*/
1057: PetscErrorCode MatSNESMFResetHHistory(Mat J)
1058: {
1059: MatSNESMFCtx ctx = (MatSNESMFCtx)J->data;
1062: ctx->ncurrenth = 0;
1063: return(0);
1064: }
1068: /*@
1069: MatSNESMFComputeJacobian - Tells the matrix-free Jacobian object the new location at which
1070: Jacobian matrix vector products will be computed at, i.e. J(x) * a.
1072: Collective on SNES
1074: Input Parameters:
1075: + snes - the nonlinear solver context
1076: . x - the point at which the Jacobian vector products will be performed
1077: . jac - the matrix-free Jacobian object
1078: . B - either the same as jac or another matrix type (ignored)
1079: . flag - not relevent for matrix-free form
1080: - dummy - the user context (ignored)
1082: Level: developer
1084: Notes:
1085: This can be passed into SNESSetJacobian() when using a completely matrix-free solver,
1086: that is the B matrix is also the same matrix operator. This is used when you select
1087: -snes_mf but rarely used directly by users.
1089: .seealso: MatSNESMFGetH(), MatCreateSNESMF(),
1090: MatSNESMFSetHHistory(),
1091: MatSNESMFKSPMonitor(), MatSNESMFSetFunctionError(), MatSNESMFCreate(), SNESSetJacobian()
1093: @*/
1094: PetscErrorCode MatSNESMFComputeJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
1095: {
1098: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
1099: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
1100: return(0);
1101: }
1105: /*@
1106: MatSNESMFSetBase - Sets the vector U at which matrix vector products of the
1107: Jacobian are computed
1109: Collective on Mat
1111: Input Parameters:
1112: + J - the MatSNESMF matrix
1113: - U - the vector
1115: Notes: This is rarely used directly
1117: Level: advanced
1119: @*/
1120: PetscErrorCode MatSNESMFSetBase(Mat J,Vec U)
1121: {
1122: PetscErrorCode ierr,(*f)(Mat,Vec);
1127: PetscObjectQueryFunction((PetscObject)J,"MatSNESMFSetBase_C",(void (**)(void))&f);
1128: if (f) {
1129: (*f)(J,U);
1130: }
1131: return(0);
1132: }
1136: /*@C
1137: MatSNESMFSetCheckh - Sets a function that checks the computed h and adjusts
1138: it to satisfy some criteria
1140: Collective on Mat
1142: Input Parameters:
1143: + J - the MatSNESMF matrix
1144: . fun - the function that checks h
1145: - ctx - any context needed by the function
1147: Options Database Keys:
1148: . -snes_mf_check_positivity
1150: Level: advanced
1152: Notes: For example, MatSNESMFSetCheckPositivity() insures that all entries
1153: of U + h*a are non-negative
1155: .seealso: MatSNESMFSetCheckPositivity()
1156: @*/
1157: PetscErrorCode MatSNESMFSetCheckh(Mat J,PetscErrorCode (*fun)(Vec,Vec,PetscScalar*,void*),void* ctx)
1158: {
1159: PetscErrorCode ierr,(*f)(Mat,PetscErrorCode (*)(Vec,Vec,PetscScalar*,void*),void*);
1163: PetscObjectQueryFunction((PetscObject)J,"MatSNESMFSetCheckh_C",(void (**)(void))&f);
1164: if (f) {
1165: (*f)(J,fun,ctx);
1166: }
1167: return(0);
1168: }
1172: /*@
1173: MatSNESMFCheckPositivity - Checks that all entries in U + h*a are positive or
1174: zero, decreases h until this is satisfied.
1176: Collective on Vec
1178: Input Parameters:
1179: + U - base vector that is added to
1180: . a - vector that is added
1181: . h - scaling factor on a
1182: - dummy - context variable (unused)
1184: Options Database Keys:
1185: . -snes_mf_check_positivity
1187: Level: advanced
1189: Notes: This is rarely used directly, rather it is passed as an argument to
1190: MatSNESMFSetCheckh()
1192: .seealso: MatSNESMFSetCheckh()
1193: @*/
1194: PetscErrorCode MatSNESMFCheckPositivity(Vec U,Vec a,PetscScalar *h,void *dummy)
1195: {
1196: PetscReal val, minval;
1197: PetscScalar *u_vec, *a_vec;
1199: PetscInt i,n;
1200: MPI_Comm comm;
1203: PetscObjectGetComm((PetscObject)U,&comm);
1204: VecGetArray(U,&u_vec);
1205: VecGetArray(a,&a_vec);
1206: VecGetLocalSize(U,&n);
1207: minval = PetscAbsScalar(*h*1.01);
1208: for(i=0;i<n;i++) {
1209: if (PetscRealPart(u_vec[i] + *h*a_vec[i]) <= 0.0) {
1210: val = PetscAbsScalar(u_vec[i]/a_vec[i]);
1211: if (val < minval) minval = val;
1212: }
1213: }
1214: VecRestoreArray(U,&u_vec);
1215: VecRestoreArray(a,&a_vec);
1216: PetscGlobalMin(&minval,&val,comm);
1217: if (val <= PetscAbsScalar(*h)) {
1218: PetscInfo2(U,"Scaling back h from %G to %G\n",PetscRealPart(*h),.99*val);
1219: if (PetscRealPart(*h) > 0.0) *h = 0.99*val;
1220: else *h = -0.99*val;
1221: }
1222: return(0);
1223: }