xref: /petsc/src/snes/mf/snesmfj.c (revision 9a6cb015ac3094bf0593987cefb2928a6c0323c2)
1*9a6cb015SBarry Smith 
2a5eb4965SSatish Balay #ifdef PETSC_RCS_HEADER
3*9a6cb015SBarry Smith static char vcid[] = "$Id: snesmfj.c,v 1.72 1998/11/05 04:28:55 bsmith Exp bsmith $";
481e6777dSBarry Smith #endif
581e6777dSBarry Smith 
6*9a6cb015SBarry Smith #include "src/snes/snesimpl.h"
7*9a6cb015SBarry Smith #include "src/snes/mf/snesmfj.h"   /*I  "snes.h"   I*/
881e6777dSBarry Smith 
9*9a6cb015SBarry Smith FList MatSNESFDMFList              = 0;
10*9a6cb015SBarry Smith int   MatSNESFDMFRegisterAllCalled = 0;
11a4d4d686SBarry Smith 
1239e2f89bSBarry Smith 
135615d1e5SSatish Balay #undef __FUNC__
14a4d4d686SBarry Smith #define __FUNC__ "MatSNESFDMFSetType"
15*9a6cb015SBarry Smith /*@
16*9a6cb015SBarry Smith       MatSNESFDMFSetType - Sets the method that is used to compute the h in the
17*9a6cb015SBarry Smith             finite difference matrix free formulation.
18*9a6cb015SBarry Smith 
19*9a6cb015SBarry Smith    Input Parameters:
20*9a6cb015SBarry Smith +     mat - the matrix free matrix created via MatCreateSNESFDMF()
21*9a6cb015SBarry Smith -     ftype - the type requested
22*9a6cb015SBarry Smith 
23*9a6cb015SBarry Smith .seealso: MatCreateSNESFDMF(), MatSNESFDMFRegister()
24*9a6cb015SBarry Smith @*/
25*9a6cb015SBarry Smith int MatSNESFDMFSetType(Mat mat,char *ftype)
26b9fa9cd0SBarry Smith {
27*9a6cb015SBarry Smith   int            ierr, (*r)(MatSNESFDMFCtx);
28a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
29a4d4d686SBarry Smith 
30a4d4d686SBarry Smith   PetscFunctionBegin;
31a4d4d686SBarry Smith   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
32a4d4d686SBarry Smith 
33*9a6cb015SBarry Smith   /* already set, so just return */
34*9a6cb015SBarry Smith   if (!PetscStrcmp(ctx->type_name,ftype)) PetscFunctionReturn(0);
35a4d4d686SBarry Smith 
36*9a6cb015SBarry Smith   /* destroy the old one if it exists */
37*9a6cb015SBarry Smith   if (ctx->ops->destroy) {
38*9a6cb015SBarry Smith     ierr = (*ctx->ops->destroy)(ctx);CHKERRQ(ierr);
39*9a6cb015SBarry Smith   }
40*9a6cb015SBarry Smith 
41*9a6cb015SBarry Smith   /* Get the function pointers for the iterative method requested */
42*9a6cb015SBarry Smith   if (!MatSNESFDMFRegisterAllCalled) {ierr = MatSNESFDMFRegisterAll(PETSC_NULL); CHKERRQ(ierr);}
43*9a6cb015SBarry Smith 
44*9a6cb015SBarry Smith   ierr =  FListFind(ctx->comm, MatSNESFDMFList, ftype,(int (**)(void *)) &r );CHKERRQ(ierr);
45*9a6cb015SBarry Smith 
46*9a6cb015SBarry Smith   if (!r) SETERRQ(1,1,"Unknown MatSNESFDMF type given");
47*9a6cb015SBarry Smith 
48*9a6cb015SBarry Smith   ierr = (*r)(ctx); CHKERRQ(ierr);
49*9a6cb015SBarry Smith   ierr = PetscStrncpy(ctx->type_name,ftype,256);CHKERRQ(ierr);
50*9a6cb015SBarry Smith 
51*9a6cb015SBarry Smith   PetscFunctionReturn(0);
52*9a6cb015SBarry Smith }
53*9a6cb015SBarry Smith 
54*9a6cb015SBarry Smith /*MC
55*9a6cb015SBarry Smith    MatSNESFDMFRegister - Adds a method to the MatSNESFDMF registry
56*9a6cb015SBarry Smith 
57*9a6cb015SBarry Smith    Synopsis:
58*9a6cb015SBarry Smith    MatSNESFDMFRegister(char *name_solver,char *path,char *name_create,int (*routine_create)(MatSNESFDMF))
59*9a6cb015SBarry Smith 
60*9a6cb015SBarry Smith    Not Collective
61*9a6cb015SBarry Smith 
62*9a6cb015SBarry Smith    Input Parameters:
63*9a6cb015SBarry Smith +  name_solver - name of a new user-defined compute-h module
64*9a6cb015SBarry Smith .  path - path (either absolute or relative) the library containing this solver
65*9a6cb015SBarry Smith .  name_create - name of routine to create method context
66*9a6cb015SBarry Smith -  routine_create - routine to create method context
67*9a6cb015SBarry Smith 
68*9a6cb015SBarry Smith    Notes:
69*9a6cb015SBarry Smith    MatSNESFDMFRegister() may be called multiple times to add several user-defined solvers.
70*9a6cb015SBarry Smith 
71*9a6cb015SBarry Smith    If dynamic libraries are used, then the fourth input argument (routine_create)
72*9a6cb015SBarry Smith    is ignored.
73*9a6cb015SBarry Smith 
74*9a6cb015SBarry Smith    Sample usage:
75*9a6cb015SBarry Smith .vb
76*9a6cb015SBarry Smith    MatSNESFDMFRegister("my_h",/home/username/my_lib/lib/libO/solaris/mylib.a,
77*9a6cb015SBarry Smith                "MyHCreate",MyHCreate);
78*9a6cb015SBarry Smith .ve
79*9a6cb015SBarry Smith 
80*9a6cb015SBarry Smith    Then, your solver can be chosen with the procedural interface via
81*9a6cb015SBarry Smith $     MatSNESFDMFSetType(mfctx,"my_h")
82*9a6cb015SBarry Smith    or at runtime via the option
83*9a6cb015SBarry Smith $     -snes_mf_type my_h
84*9a6cb015SBarry Smith 
85*9a6cb015SBarry Smith .keywords: MatSNESFDMF, register
86*9a6cb015SBarry Smith 
87*9a6cb015SBarry Smith .seealso: MatSNESFDMFRegisterAll(), MatSNESFDMFRegisterDestroy()
88*9a6cb015SBarry Smith M*/
89*9a6cb015SBarry Smith 
90*9a6cb015SBarry Smith #undef __FUNC__
91*9a6cb015SBarry Smith #define __FUNC__ "MatSNESFDMFRegister_Private"
92*9a6cb015SBarry Smith int MatSNESFDMFRegister_Private(char *sname,char *path,char *name,int (*function)(MatSNESFDMFCtx))
93*9a6cb015SBarry Smith {
94*9a6cb015SBarry Smith   int ierr;
95*9a6cb015SBarry Smith   char fullname[256];
96*9a6cb015SBarry Smith 
97*9a6cb015SBarry Smith   PetscFunctionBegin;
98*9a6cb015SBarry Smith   PetscStrcpy(fullname,path); PetscStrcat(fullname,":");PetscStrcat(fullname,name);
99*9a6cb015SBarry Smith   ierr = FListAdd_Private(&MatSNESFDMFList,sname,fullname,(int (*)(void*))function);CHKERRQ(ierr);
100*9a6cb015SBarry Smith   PetscFunctionReturn(0);
101*9a6cb015SBarry Smith }
102*9a6cb015SBarry Smith 
103*9a6cb015SBarry Smith 
104*9a6cb015SBarry Smith #undef __FUNC__
105*9a6cb015SBarry Smith #define __FUNC__ "MatSNESFDMFRegisterDestroy"
106*9a6cb015SBarry Smith /*@C
107*9a6cb015SBarry Smith    MatSNESFDMFRegisterDestroy - Frees the list of MatSNESFDMF methods that were
108*9a6cb015SBarry Smith    registered by MatSNESFDMFRegister().
109*9a6cb015SBarry Smith 
110*9a6cb015SBarry Smith    Not Collective
111*9a6cb015SBarry Smith 
112*9a6cb015SBarry Smith .keywords: MatSNESFDMF, register, destroy
113*9a6cb015SBarry Smith 
114*9a6cb015SBarry Smith .seealso: MatSNESFDMFRegister(), MatSNESFDMFRegisterAll()
115*9a6cb015SBarry Smith @*/
116*9a6cb015SBarry Smith int MatSNESFDMFRegisterDestroy(void)
117*9a6cb015SBarry Smith {
118*9a6cb015SBarry Smith   int ierr;
119*9a6cb015SBarry Smith 
120*9a6cb015SBarry Smith   PetscFunctionBegin;
121*9a6cb015SBarry Smith   if (MatSNESFDMFList) {
122*9a6cb015SBarry Smith     ierr = FListDestroy( MatSNESFDMFList );CHKERRQ(ierr);
123*9a6cb015SBarry Smith     MatSNESFDMFList = 0;
124*9a6cb015SBarry Smith   }
125*9a6cb015SBarry Smith   MatSNESFDMFRegisterAllCalled = 0;
126*9a6cb015SBarry Smith   PetscFunctionReturn(0);
127*9a6cb015SBarry Smith }
128*9a6cb015SBarry Smith 
129*9a6cb015SBarry Smith /* ----------------------------------------------------------------------------------------*/
130a4d4d686SBarry Smith #undef __FUNC__
131a4d4d686SBarry Smith #define __FUNC__ "MatSNESFDMFDestroy_Private"
132a4d4d686SBarry Smith int MatSNESFDMFDestroy_Private(Mat mat)
133a4d4d686SBarry Smith {
134a4d4d686SBarry Smith   int            ierr;
135a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
136fae171e0SBarry Smith 
1373a40ed3dSBarry Smith   PetscFunctionBegin;
1380a25c783SBarry Smith   ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr);
139b9fa9cd0SBarry Smith   ierr = VecDestroy(ctx->w); CHKERRQ(ierr);
140*9a6cb015SBarry Smith   if (ctx->ops->destroy) {ierr = (*ctx->ops->destroy)(ctx);CHKERRQ(ierr);}
141b4fd4287SBarry Smith   if (ctx->sp) {ierr = PCNullSpaceDestroy(ctx->sp);CHKERRQ(ierr);}
142*9a6cb015SBarry Smith   PetscFree(ctx->ops);
143fae171e0SBarry Smith   PetscFree(ctx);
1443a40ed3dSBarry Smith   PetscFunctionReturn(0);
145b9fa9cd0SBarry Smith }
14650361f65SLois Curfman McInnes 
1475615d1e5SSatish Balay #undef __FUNC__
148a4d4d686SBarry Smith #define __FUNC__ "MatSNESFDMFView_Private"
14939e2f89bSBarry Smith /*
150a4d4d686SBarry Smith    MatSNESFDMFView_Private - Views matrix-free parameters.
1518f6e3e37SBarry Smith 
15239e2f89bSBarry Smith */
153a4d4d686SBarry Smith int MatSNESFDMFView_Private(Mat J,Viewer viewer)
154eb9086c3SLois Curfman McInnes {
155eb9086c3SLois Curfman McInnes   int            ierr;
156a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
157eb9086c3SLois Curfman McInnes   MPI_Comm       comm;
158eb9086c3SLois Curfman McInnes   FILE           *fd;
159eb9086c3SLois Curfman McInnes   ViewerType     vtype;
160eb9086c3SLois Curfman McInnes 
1613a40ed3dSBarry Smith   PetscFunctionBegin;
162eb9086c3SLois Curfman McInnes   PetscObjectGetComm((PetscObject)J,&comm);
163eb9086c3SLois Curfman McInnes   ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr);
164eb9086c3SLois Curfman McInnes   ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr);
165eb9086c3SLois Curfman McInnes   ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr);
166eb9086c3SLois Curfman McInnes   if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) {
167eb9086c3SLois Curfman McInnes      PetscFPrintf(comm,fd,"  SNES matrix-free approximation:\n");
168eb9086c3SLois Curfman McInnes      PetscFPrintf(comm,fd,"    err=%g (relative error in function evaluation)\n",ctx->error_rel);
169*9a6cb015SBarry Smith      PetscFPrintf(ctx->comm,fd,"    Using %s compute h routine \n",ctx->type_name);
170*9a6cb015SBarry Smith      if (ctx->ops->view) {
171*9a6cb015SBarry Smith        ierr = (*ctx->ops->view)(ctx,viewer);CHKERRQ(ierr);
172*9a6cb015SBarry Smith      }
1735cd90555SBarry Smith   } else {
1745cd90555SBarry Smith     SETERRQ(1,1,"Viewer type not supported for this object");
175eb9086c3SLois Curfman McInnes   }
1763a40ed3dSBarry Smith   PetscFunctionReturn(0);
177eb9086c3SLois Curfman McInnes }
178eb9086c3SLois Curfman McInnes 
179c481317fSBarry Smith 
1805615d1e5SSatish Balay #undef __FUNC__
181a4d4d686SBarry Smith #define __FUNC__ "MatSNESFDMFMult_Private"
182eb9086c3SLois Curfman McInnes /*
183a4d4d686SBarry Smith   MatSNESFDMFMult_Private - Default matrix-free form for Jacobian-vector
184eb9086c3SLois Curfman McInnes   product, y = F'(u)*a:
185a4d4d686SBarry Smith 
186*9a6cb015SBarry Smith         y ~= ( F(u + ha) - F(u) )/h,
187eb9086c3SLois Curfman McInnes   where F = nonlinear function, as set by SNESSetFunction()
188eb9086c3SLois Curfman McInnes         u = current iterate
189eb9086c3SLois Curfman McInnes         h = difference interval
190eb9086c3SLois Curfman McInnes */
191a4d4d686SBarry Smith int MatSNESFDMFMult_Private(Mat mat,Vec a,Vec y)
19239e2f89bSBarry Smith {
193a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
194fae171e0SBarry Smith   SNES           snes;
195a4d4d686SBarry Smith   Scalar         h, mone = -1.0;
196fae171e0SBarry Smith   Vec            w,U,F;
19783e56afdSLois Curfman McInnes   int            ierr, (*eval_fct)(SNES,Vec,Vec);
19839e2f89bSBarry Smith 
1993a40ed3dSBarry Smith   PetscFunctionBegin;
200*9a6cb015SBarry Smith   /* We log matrix-free matrix-vector products separately, so that we can
201*9a6cb015SBarry Smith      separate the performance monitoring from the cases that use conventional
202*9a6cb015SBarry Smith      storage.  We may eventually modify event logging to associate events
203*9a6cb015SBarry Smith      with particular objects, hence alleviating the more general problem. */
20456cd22aeSBarry Smith   PLogEventBegin(MAT_MatrixFreeMult,a,y,0,0);
20556cd22aeSBarry Smith 
2066e38a044SBarry Smith   ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr);
207fae171e0SBarry Smith   snes = ctx->snes;
208fae171e0SBarry Smith   w    = ctx->w;
209fae171e0SBarry Smith 
21078b31e54SBarry Smith   ierr = SNESGetSolution(snes,&U); CHKERRQ(ierr);
21183e56afdSLois Curfman McInnes   if (snes->method_class == SNES_NONLINEAR_EQUATIONS) {
21283e56afdSLois Curfman McInnes     eval_fct = SNESComputeFunction;
21378b31e54SBarry Smith     ierr = SNESGetFunction(snes,&F); CHKERRQ(ierr);
214a4d4d686SBarry Smith   } else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) {
21583e56afdSLois Curfman McInnes     eval_fct = SNESComputeGradient;
21683e56afdSLois Curfman McInnes     ierr = SNESGetGradient(snes,&F); CHKERRQ(ierr);
217a4d4d686SBarry Smith   } else SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Invalid method class");
21850361f65SLois Curfman McInnes 
219*9a6cb015SBarry Smith   if (!ctx->ops->compute) {
220*9a6cb015SBarry Smith     ierr = MatSNESFDMFSetType(mat,"default");CHKERRQ(ierr);
221*9a6cb015SBarry Smith     ierr = MatSNESFDMFSetFromOptions(mat);CHKERRQ(ierr);
222*9a6cb015SBarry Smith   }
223*9a6cb015SBarry Smith   ierr = (*ctx->ops->compute)(ctx,U,a,&h); CHKERRQ(ierr);
224a4d4d686SBarry Smith 
225a4d4d686SBarry Smith   /* keep a record of the current differencing parameter h */
226a4d4d686SBarry Smith   ctx->currenth = h;
227a4d4d686SBarry Smith #if defined(USE_PETSC_COMPLEX)
228a4d4d686SBarry Smith   PLogInfo(mat,"Current differencing parameter: %g + %g i\n",PetscReal(h),PetscImaginary(h));
229a4d4d686SBarry Smith #else
230a4d4d686SBarry Smith   PLogInfo(mat,"Current differencing parameter: %g\n",h);
231a4d4d686SBarry Smith #endif
232a4d4d686SBarry Smith   if (ctx->historyh && ctx->ncurrenth < ctx->maxcurrenth) {
233a4d4d686SBarry Smith     ctx->historyh[ctx->ncurrenth++] = h;
234a4d4d686SBarry Smith   }
235a4d4d686SBarry Smith 
236a4d4d686SBarry Smith   /* Evaluate function at F(u + ha) */
237a4d4d686SBarry Smith   ierr = VecWAXPY(&h,a,U,w); CHKERRQ(ierr);
238a4d4d686SBarry Smith   ierr = eval_fct(snes,w,y); CHKERRQ(ierr);
239a4d4d686SBarry Smith 
240a4d4d686SBarry Smith   ierr = VecAXPY(&mone,F,y); CHKERRQ(ierr);
241a4d4d686SBarry Smith   h = 1.0/h;
242a4d4d686SBarry Smith   ierr = VecScale(&h,y); CHKERRQ(ierr);
243a4d4d686SBarry Smith   if (ctx->sp) {ierr = PCNullSpaceRemove(ctx->sp,y); CHKERRQ(ierr);}
244a4d4d686SBarry Smith 
245a4d4d686SBarry Smith   PLogEventEnd(MAT_MatrixFreeMult,a,y,0,0);
246a4d4d686SBarry Smith   PetscFunctionReturn(0);
247a4d4d686SBarry Smith }
248a4d4d686SBarry Smith 
249a4d4d686SBarry Smith #undef __FUNC__
250a4d4d686SBarry Smith #define __FUNC__ "MatCreateSNESFDMF"
251a4d4d686SBarry Smith /*@C
252a4d4d686SBarry Smith    MatCreateSNESFDMF - Creates a matrix-free matrix
253a4d4d686SBarry Smith    context for use with a SNES solver.  This matrix can be used as
254a4d4d686SBarry Smith    the Jacobian argument for the routine SNESSetJacobian().
255a4d4d686SBarry Smith 
256a4d4d686SBarry Smith    Collective on SNES and Vec
257a4d4d686SBarry Smith 
258a4d4d686SBarry Smith    Input Parameters:
259a4d4d686SBarry Smith +  snes - the SNES context
260a4d4d686SBarry Smith -  x - vector where SNES solution is to be stored.
261a4d4d686SBarry Smith 
262a4d4d686SBarry Smith    Output Parameter:
263a4d4d686SBarry Smith .  J - the matrix-free matrix
264a4d4d686SBarry Smith 
265a4d4d686SBarry Smith    Notes:
266a4d4d686SBarry Smith    The matrix-free matrix context merely contains the function pointers
267a4d4d686SBarry Smith    and work space for performing finite difference approximations of
268*9a6cb015SBarry Smith    Jacobian-vector products, J(u)*a,
269*9a6cb015SBarry Smith 
270*9a6cb015SBarry Smith    The default code uses the following approach to compute h
271a4d4d686SBarry Smith 
272a4d4d686SBarry Smith .vb
273a4d4d686SBarry Smith      J(u)*a = [J(u+h*a) - J(u)]/h where
274a4d4d686SBarry Smith      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
275a4d4d686SBarry Smith        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   otherwise
276a4d4d686SBarry Smith  where
277a4d4d686SBarry Smith      error_rel = square root of relative error in function evaluation
278a4d4d686SBarry Smith      umin = minimum iterate parameter
279a4d4d686SBarry Smith .ve
280a4d4d686SBarry Smith 
281*9a6cb015SBarry Smith    The user can set the error_rel via MatSNESFDMFSetFunctionError() and
282*9a6cb015SBarry Smith    umin via MatSNESFDMFDefaultSetUmin()
283a4d4d686SBarry Smith    See the nonlinear solvers chapter of the users manual for details.
284a4d4d686SBarry Smith 
285a4d4d686SBarry Smith    The user should call MatDestroy() when finished with the matrix-free
286a4d4d686SBarry Smith    matrix context.
287a4d4d686SBarry Smith 
288a4d4d686SBarry Smith    Options Database Keys:
289a4d4d686SBarry Smith +  -snes_mf_err <error_rel> - Sets error_rel
290*9a6cb015SBarry Smith .  -snes_mf_unim <umin> - Sets umin (for default PETSc routine that computes h only)
291a4d4d686SBarry Smith -  -snes_mf_ksp_monitor - KSP monitor routine that prints differencing h
292a4d4d686SBarry Smith 
293a4d4d686SBarry Smith .keywords: SNES, default, matrix-free, create, matrix
294a4d4d686SBarry Smith 
295*9a6cb015SBarry Smith .seealso: MatDestroy(), MatSNESFDMFSetFunctionError(), MatSNESFDMFDefaultSetUmin()
296a4d4d686SBarry Smith           MatSNESFDMFSetHHistory(), MatSNESFDMFResetHHistory(),
297*9a6cb015SBarry Smith           MatSNESFDMFGetH(),MatSNESFDMFKSPMonitor(), MatSNESFDMFRegister()
298a4d4d686SBarry Smith 
299a4d4d686SBarry Smith @*/
300a4d4d686SBarry Smith int MatCreateSNESFDMF(SNES snes,Vec x, Mat *J)
301a4d4d686SBarry Smith {
302a4d4d686SBarry Smith   MPI_Comm       comm;
303a4d4d686SBarry Smith   MatSNESFDMFCtx mfctx;
304*9a6cb015SBarry Smith   int            n, nloc, ierr;
305a4d4d686SBarry Smith 
306a4d4d686SBarry Smith   PetscFunctionBegin;
307a4d4d686SBarry Smith   mfctx = (MatSNESFDMFCtx) PetscMalloc(sizeof(struct _p_MatSNESFDMFCtx)); CHKPTRQ(mfctx);
308a4d4d686SBarry Smith   PLogObjectMemory(snes,sizeof(MatSNESFDMFCtx));
309*9a6cb015SBarry Smith   mfctx->comm         = snes->comm;
310a4d4d686SBarry Smith   mfctx->sp           = 0;
311a4d4d686SBarry Smith   mfctx->snes         = snes;
312a4d4d686SBarry Smith   mfctx->error_rel    = 1.e-8; /* assumes double precision */
313a4d4d686SBarry Smith   mfctx->currenth     = 0.0;
314a4d4d686SBarry Smith   mfctx->historyh     = PETSC_NULL;
315a4d4d686SBarry Smith   mfctx->ncurrenth    = 0;
316a4d4d686SBarry Smith   mfctx->maxcurrenth  = 0;
317a4d4d686SBarry Smith 
318*9a6cb015SBarry Smith   /*
319*9a6cb015SBarry Smith      Create the empty data structure to contain compute-h routines.
320*9a6cb015SBarry Smith      These will be filled in below from the command line options or
321*9a6cb015SBarry Smith      a later call with MatSNESFDMFSetType() or if that is not called
322*9a6cb015SBarry Smith      then it will default in the first use of MatSNESFDMFMult_private()
323*9a6cb015SBarry Smith   */
324*9a6cb015SBarry Smith   mfctx->ops                 = (MFOps *)PetscMalloc(sizeof(MFOps)); CHKPTRQ(mfctx->ops);
325*9a6cb015SBarry Smith   mfctx->ops->compute        = 0;
326*9a6cb015SBarry Smith   mfctx->ops->destroy        = 0;
327*9a6cb015SBarry Smith   mfctx->ops->view           = 0;
328*9a6cb015SBarry Smith   mfctx->ops->printhelp      = 0;
329*9a6cb015SBarry Smith   mfctx->ops->setfromoptions = 0;
330*9a6cb015SBarry Smith   mfctx->hctx                = 0;
331*9a6cb015SBarry Smith 
332a4d4d686SBarry Smith   ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr);
333a4d4d686SBarry Smith   ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr);
334a4d4d686SBarry Smith   ierr = VecGetSize(x,&n); CHKERRQ(ierr);
335a4d4d686SBarry Smith   ierr = VecGetLocalSize(x,&nloc); CHKERRQ(ierr);
336a4d4d686SBarry Smith   ierr = MatCreateShell(comm,nloc,n,n,n,mfctx,J); CHKERRQ(ierr);
337a4d4d686SBarry Smith   ierr = MatShellSetOperation(*J,MATOP_MULT,(void*)MatSNESFDMFMult_Private);CHKERRQ(ierr);
338a4d4d686SBarry Smith   ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void *)MatSNESFDMFDestroy_Private);CHKERRQ(ierr);
339a4d4d686SBarry Smith   ierr = MatShellSetOperation(*J,MATOP_VIEW,(void *)MatSNESFDMFView_Private); CHKERRQ(ierr);
340a4d4d686SBarry Smith   PLogObjectParent(*J,mfctx->w);
341a4d4d686SBarry Smith   PLogObjectParent(snes,*J);
342*9a6cb015SBarry Smith 
343*9a6cb015SBarry Smith   mfctx->mat = *J;
344*9a6cb015SBarry Smith 
345*9a6cb015SBarry Smith 
346*9a6cb015SBarry Smith   PetscFunctionReturn(0);
347*9a6cb015SBarry Smith }
348*9a6cb015SBarry Smith 
349*9a6cb015SBarry Smith #undef __FUNC__
350*9a6cb015SBarry Smith #define __FUNC__ "MatSNESFDMFGetH"
351*9a6cb015SBarry Smith /*@
352*9a6cb015SBarry Smith    MatSNESFDMFSetFromOptions - Sets the MatSNESFDMF options from the command line
353*9a6cb015SBarry Smith      parameter.
354*9a6cb015SBarry Smith 
355*9a6cb015SBarry Smith    Collective on Mat
356*9a6cb015SBarry Smith 
357*9a6cb015SBarry Smith    Input Parameters:
358*9a6cb015SBarry Smith .   mat - the matrix obtained with MatCreateSNESFDMF()
359*9a6cb015SBarry Smith 
360*9a6cb015SBarry Smith .keywords: SNES, matrix-free, parameters
361*9a6cb015SBarry Smith 
362*9a6cb015SBarry Smith .seealso: MatCreateSNESFDMF(),MatSNESFDMFSetHHistory(),
363*9a6cb015SBarry Smith           MatSNESFDMFResetHHistory(),MatSNESFDMFKSPMonitor()
364*9a6cb015SBarry Smith @*/
365*9a6cb015SBarry Smith int MatSNESFDMFSetFromOptions(Mat mat)
366*9a6cb015SBarry Smith {
367*9a6cb015SBarry Smith   MatSNESFDMFCtx mfctx;
368*9a6cb015SBarry Smith   int            ierr,flg;
369*9a6cb015SBarry Smith   char           ftype[256],p[64];
370*9a6cb015SBarry Smith 
371*9a6cb015SBarry Smith   PetscFunctionBegin;
372*9a6cb015SBarry Smith   ierr = MatShellGetContext(mat,(void **)&mfctx); CHKERRQ(ierr);
373*9a6cb015SBarry Smith   if (mfctx) {
374*9a6cb015SBarry Smith     /* allow user to set the type */
375*9a6cb015SBarry Smith     ierr = OptionsGetString(mfctx->snes->prefix,"-snes_mf_type",ftype,256,&flg);CHKERRQ(ierr);
376*9a6cb015SBarry Smith     if (flg) {
377*9a6cb015SBarry Smith       ierr = MatSNESFDMFSetType(mat,ftype);CHKERRQ(ierr);
378*9a6cb015SBarry Smith     }
379*9a6cb015SBarry Smith 
380*9a6cb015SBarry Smith     ierr = OptionsGetDouble(mfctx->snes->prefix,"-snes_mf_err",&mfctx->error_rel,&flg);CHKERRQ(ierr);
381*9a6cb015SBarry Smith     if (mfctx->ops->setfromoptions) {
382*9a6cb015SBarry Smith       ierr = (*mfctx->ops->setfromoptions)(mfctx);CHKERRQ(ierr);
383*9a6cb015SBarry Smith     }
384*9a6cb015SBarry Smith 
385*9a6cb015SBarry Smith     ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr);
386*9a6cb015SBarry Smith     PetscStrcpy(p,"-");
387*9a6cb015SBarry Smith     if (mfctx->snes->prefix) PetscStrcat(p,mfctx->snes->prefix);
388*9a6cb015SBarry Smith     if (flg) {
389*9a6cb015SBarry Smith       (*PetscHelpPrintf)(mfctx->snes->comm,"   %ssnes_mf_err <err>: set sqrt rel error in function (default %g)\n",p,mfctx->error_rel);
390*9a6cb015SBarry Smith       if (mfctx->ops->printhelp) {
391*9a6cb015SBarry Smith         (*mfctx->ops->printhelp)(mfctx);CHKERRQ(ierr);
392*9a6cb015SBarry Smith       }
393*9a6cb015SBarry Smith     }
394*9a6cb015SBarry Smith   }
395a4d4d686SBarry Smith   PetscFunctionReturn(0);
396a4d4d686SBarry Smith }
397a4d4d686SBarry Smith 
398a4d4d686SBarry Smith #undef __FUNC__
399a4d4d686SBarry Smith #define __FUNC__ "MatSNESFDMFGetH"
400a4d4d686SBarry Smith /*@
401a4d4d686SBarry Smith    MatSNESFDMFGetH - Gets the last h that was used as the differencing
402a4d4d686SBarry Smith      parameter.
403a4d4d686SBarry Smith 
404a4d4d686SBarry Smith    Not Collective
405a4d4d686SBarry Smith 
406a4d4d686SBarry Smith    Input Parameters:
407a4d4d686SBarry Smith .   mat - the matrix obtained with MatCreateSNESFDMF()
408a4d4d686SBarry Smith 
409a4d4d686SBarry Smith    Output Paramter:
410a4d4d686SBarry Smith .  h - the differencing step size
411a4d4d686SBarry Smith 
412a4d4d686SBarry Smith .keywords: SNES, matrix-free, parameters
413a4d4d686SBarry Smith 
414a4d4d686SBarry Smith .seealso: MatCreateSNESFDMF(),MatSNESFDMFSetHHistory(),
415a4d4d686SBarry Smith           MatSNESFDMFResetHHistory(),MatSNESFDMFKSPMonitor()
416a4d4d686SBarry Smith @*/
417a4d4d686SBarry Smith int MatSNESFDMFGetH(Mat mat,Scalar *h)
418a4d4d686SBarry Smith {
419a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
420a4d4d686SBarry Smith   int            ierr;
421a4d4d686SBarry Smith 
422a4d4d686SBarry Smith   PetscFunctionBegin;
423a4d4d686SBarry Smith   ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr);
424a4d4d686SBarry Smith   if (ctx) {
425a4d4d686SBarry Smith     *h = ctx->currenth;
426a4d4d686SBarry Smith   }
427a4d4d686SBarry Smith   PetscFunctionReturn(0);
428a4d4d686SBarry Smith }
429a4d4d686SBarry Smith 
430a4d4d686SBarry Smith #undef __FUNC__
431a4d4d686SBarry Smith #define __FUNC__ "MatSNESFDMFKSPMonitor"
432a4d4d686SBarry Smith /*
433a4d4d686SBarry Smith    MatSNESFDMFKSPMonitor - A KSP monitor for use with the default PETSc
434a4d4d686SBarry Smith       SNES matrix free routines. Prints the h differencing parameter used at each
435a4d4d686SBarry Smith       timestep.
436a4d4d686SBarry Smith 
437a4d4d686SBarry Smith */
438a4d4d686SBarry Smith int MatSNESFDMFKSPMonitor(KSP ksp,int n,double rnorm,void *dummy)
439a4d4d686SBarry Smith {
440a4d4d686SBarry Smith   PC             pc;
441a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
442a4d4d686SBarry Smith   int            ierr;
443a4d4d686SBarry Smith   Mat            mat;
444a4d4d686SBarry Smith   MPI_Comm       comm;
445a4d4d686SBarry Smith   PetscTruth     nonzeroinitialguess;
446a4d4d686SBarry Smith 
447a4d4d686SBarry Smith   PetscFunctionBegin;
448a4d4d686SBarry Smith   ierr = PetscObjectGetComm((PetscObject)ksp,&comm);CHKERRQ(ierr);
449a4d4d686SBarry Smith   ierr = KSPGetPC(ksp,&pc); CHKERRQ(ierr);
450a4d4d686SBarry Smith   ierr = KSPGetInitialGuessNonzero(ksp,&nonzeroinitialguess);CHKERRQ(ierr);
451a4d4d686SBarry Smith   ierr = PCGetOperators(pc,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
452a4d4d686SBarry Smith   ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr);
453*9a6cb015SBarry Smith   if (!ctx) {
454*9a6cb015SBarry Smith     SETERRQ(1,1,"Matrix is not a matrix free shell matrix");
455*9a6cb015SBarry Smith   }
456a4d4d686SBarry Smith   if (n > 0 || nonzeroinitialguess) {
457a4d4d686SBarry Smith #if defined(USE_PETSC_COMPLEX)
458a4d4d686SBarry Smith     PetscPrintf(comm,"%d KSP Residual norm %14.12e h %g + %g i\n",n,rnorm,
459a4d4d686SBarry Smith                 PetscReal(ctx->currenth),PetscImaginary(ctx->currenth));
460a4d4d686SBarry Smith #else
461a4d4d686SBarry Smith     PetscPrintf(comm,"%d KSP Residual norm %14.12e h %g \n",n,rnorm,ctx->currenth);
462a4d4d686SBarry Smith #endif
463a4d4d686SBarry Smith   } else {
464a4d4d686SBarry Smith     PetscPrintf(comm,"%d KSP Residual norm %14.12e\n",n,rnorm);
465a4d4d686SBarry Smith   }
466a4d4d686SBarry Smith   PetscFunctionReturn(0);
467a4d4d686SBarry Smith }
468a4d4d686SBarry Smith 
469a4d4d686SBarry Smith #undef __FUNC__
470*9a6cb015SBarry Smith #define __FUNC__ "MatSNESFDMFSetFunctionError"
471a4d4d686SBarry Smith /*@
472*9a6cb015SBarry Smith    MatSNESFDMFSetFunctionError - Sets the error_rel for the approximation of
473a4d4d686SBarry Smith    matrix-vector products using finite differences.
474a4d4d686SBarry Smith 
475a4d4d686SBarry Smith    Collective on Mat
476a4d4d686SBarry Smith 
477a4d4d686SBarry Smith    Input Parameters:
478a4d4d686SBarry Smith +  mat - the matrix free matrix created via MatCreateSNESFDMF()
479*9a6cb015SBarry Smith -  error_rel - relative error (should be set to the square root of
480a4d4d686SBarry Smith                the relative error in the function evaluations)
481a4d4d686SBarry Smith 
482a4d4d686SBarry Smith    Notes:
483a4d4d686SBarry Smith    The default matrix-free matrix-vector product routine computes
484a4d4d686SBarry Smith .vb
485a4d4d686SBarry Smith      J(u)*a = [J(u+h*a) - J(u)]/h where
486a4d4d686SBarry Smith      h = error_rel*u'a/||a||^2                        if  |u'a| > umin*||a||_{1}
487a4d4d686SBarry Smith        = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2   else
488a4d4d686SBarry Smith .ve
489a4d4d686SBarry Smith 
490a4d4d686SBarry Smith    Options Database Keys:
491a4d4d686SBarry Smith +  -snes_mf_err <error_rel> - Sets error_rel
492a4d4d686SBarry Smith 
493a4d4d686SBarry Smith .keywords: SNES, matrix-free, parameters
494a4d4d686SBarry Smith 
495a4d4d686SBarry Smith .seealso: MatCreateSNESFDMF(),MatSNESFDMFGetH(),
496a4d4d686SBarry Smith           MatSNESFDMFSetHHistory(), MatSNESFDMFResetHHistory(),
497a4d4d686SBarry Smith           MatSNESFDMFKSPMonitor()
498a4d4d686SBarry Smith @*/
499*9a6cb015SBarry Smith int MatSNESFDMFSetFunctionError(Mat mat,double error)
500a4d4d686SBarry Smith {
501a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
502a4d4d686SBarry Smith   int            ierr;
503a4d4d686SBarry Smith 
504a4d4d686SBarry Smith   PetscFunctionBegin;
505a4d4d686SBarry Smith   ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr);
506a4d4d686SBarry Smith   if (ctx) {
507a4d4d686SBarry Smith     if (error != PETSC_DEFAULT) ctx->error_rel = error;
508a4d4d686SBarry Smith   }
509a4d4d686SBarry Smith   PetscFunctionReturn(0);
510a4d4d686SBarry Smith }
511a4d4d686SBarry Smith 
512a4d4d686SBarry Smith #undef __FUNC__
513a4d4d686SBarry Smith #define __FUNC__ "MatSNESFDMFAddNullSpace"
514a4d4d686SBarry Smith /*@
515a4d4d686SBarry Smith    MatSNESFDMFAddNullSpace - Provides a null space that
516a4d4d686SBarry Smith    an operator is supposed to have.  Since roundoff will create a
517a4d4d686SBarry Smith    small component in the null space, if you know the null space
518a4d4d686SBarry Smith    you may have it automatically removed.
519a4d4d686SBarry Smith 
520a4d4d686SBarry Smith    Collective on Mat
521a4d4d686SBarry Smith 
522a4d4d686SBarry Smith    Input Parameters:
523a4d4d686SBarry Smith +  J - the matrix-free matrix context
524a4d4d686SBarry Smith .  has_cnst - PETSC_TRUE or PETSC_FALSE, indicating if null space has constants
525a4d4d686SBarry Smith .  n - number of vectors (excluding constant vector) in null space
526a4d4d686SBarry Smith -  vecs - the vectors that span the null space (excluding the constant vector);
527a4d4d686SBarry Smith           these vectors must be orthonormal
528a4d4d686SBarry Smith 
529a4d4d686SBarry Smith .keywords: SNES, matrix-free, null space
530a4d4d686SBarry Smith 
531a4d4d686SBarry Smith .seealso: MatSNESFDMFGetH(), MatCreateSNESFDMF(),
532a4d4d686SBarry Smith           MatSNESFDMFSetHHistory(), MatSNESFDMFResetHHistory(),
533*9a6cb015SBarry Smith           MatSNESFDMFKSPMonitor(), MatSNESFDMFErrorRel()
534a4d4d686SBarry Smith 
535a4d4d686SBarry Smith @*/
536a4d4d686SBarry Smith int MatSNESFDMFAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs)
537a4d4d686SBarry Smith {
538a4d4d686SBarry Smith   int            ierr;
539a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
540a4d4d686SBarry Smith   MPI_Comm       comm;
541a4d4d686SBarry Smith 
542a4d4d686SBarry Smith   PetscFunctionBegin;
543a4d4d686SBarry Smith   PetscObjectGetComm((PetscObject)J,&comm);
544a4d4d686SBarry Smith 
545a4d4d686SBarry Smith   ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr);
546a4d4d686SBarry Smith   /* no context indicates that it is not the "matrix free" matrix type */
547a4d4d686SBarry Smith   if (!ctx) PetscFunctionReturn(0);
548a4d4d686SBarry Smith   ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr);
549a4d4d686SBarry Smith   PetscFunctionReturn(0);
550a4d4d686SBarry Smith }
551a4d4d686SBarry Smith 
552a4d4d686SBarry Smith #undef __FUNC__
553a4d4d686SBarry Smith #define __FUNC__ "MatSNESFDMFSetHHistory"
554a4d4d686SBarry Smith /*@
555a4d4d686SBarry Smith    MatSNESFDMFSetHHistory - Sets an array to collect a history
556a4d4d686SBarry Smith       of the differencing values h computed for the matrix free product
557a4d4d686SBarry Smith 
558a4d4d686SBarry Smith    Collective on Mat
559a4d4d686SBarry Smith 
560a4d4d686SBarry Smith    Input Parameters:
561a4d4d686SBarry Smith +  J - the matrix-free matrix context
562a4d4d686SBarry Smith .  histroy - space to hold the h history
563a4d4d686SBarry Smith -  nhistory - number of entries in history, if more h are generated than
564a4d4d686SBarry Smith               nhistory the later ones are discarded
565a4d4d686SBarry Smith 
566a4d4d686SBarry Smith    Notes:
567a4d4d686SBarry Smith     Use MatSNESFDMFResetHHistory() to reset the history counter
568a4d4d686SBarry Smith     and collect a new batch of h.
569a4d4d686SBarry Smith 
570a4d4d686SBarry Smith .keywords: SNES, matrix-free, h history, differencing history
571a4d4d686SBarry Smith 
572a4d4d686SBarry Smith .seealso: MatSNESFDMFGetH(), MatCreateSNESFDMF(),
573a4d4d686SBarry Smith           MatSNESFDMFResetHHistory(),
574*9a6cb015SBarry Smith           MatSNESFDMFKSPMonitor(), MatSNESFDMFSetFunctionError()
575a4d4d686SBarry Smith 
576a4d4d686SBarry Smith @*/
577a4d4d686SBarry Smith int MatSNESFDMFSetHHistory(Mat J,Scalar *history,int nhistory)
578a4d4d686SBarry Smith {
579a4d4d686SBarry Smith   int            ierr;
580a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
581a4d4d686SBarry Smith 
582a4d4d686SBarry Smith   PetscFunctionBegin;
583a4d4d686SBarry Smith 
584a4d4d686SBarry Smith   ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr);
585a4d4d686SBarry Smith   /* no context indicates that it is not the "matrix free" matrix type */
586a4d4d686SBarry Smith   if (!ctx) PetscFunctionReturn(0);
587a4d4d686SBarry Smith   ctx->historyh    = history;
588a4d4d686SBarry Smith   ctx->maxcurrenth = nhistory;
589a4d4d686SBarry Smith   ctx->currenth    = 0;
590a4d4d686SBarry Smith 
591a4d4d686SBarry Smith   PetscFunctionReturn(0);
592a4d4d686SBarry Smith }
593a4d4d686SBarry Smith 
594a4d4d686SBarry Smith #undef __FUNC__
595a4d4d686SBarry Smith #define __FUNC__ "MatSNESFDMFResetHHistory"
596a4d4d686SBarry Smith /*@
597a4d4d686SBarry Smith    MatSNESFDMFResetHHistory - Resets the counter to zero to begin
598a4d4d686SBarry Smith       collecting a new set of differencing histories.
599a4d4d686SBarry Smith 
600a4d4d686SBarry Smith    Collective on Mat
601a4d4d686SBarry Smith 
602a4d4d686SBarry Smith    Input Parameters:
603a4d4d686SBarry Smith .  J - the matrix-free matrix context
604a4d4d686SBarry Smith 
605a4d4d686SBarry Smith    Notes:
606a4d4d686SBarry Smith     Use MatSNESFDMFSetHHistory() to create the original history counter
607a4d4d686SBarry Smith 
608a4d4d686SBarry Smith .keywords: SNES, matrix-free, h history, differencing history
609a4d4d686SBarry Smith 
610a4d4d686SBarry Smith .seealso: MatSNESFDMFGetH(), MatCreateSNESFDMF(),
611a4d4d686SBarry Smith           MatSNESFDMFSetHHistory(),
612*9a6cb015SBarry Smith           MatSNESFDMFKSPMonitor(), MatSNESFDMFSetFunctionError()
613a4d4d686SBarry Smith 
614a4d4d686SBarry Smith @*/
615a4d4d686SBarry Smith int MatSNESFDMFResetHHistory(Mat J)
616a4d4d686SBarry Smith {
617a4d4d686SBarry Smith   int            ierr;
618a4d4d686SBarry Smith   MatSNESFDMFCtx ctx;
619a4d4d686SBarry Smith 
620a4d4d686SBarry Smith   PetscFunctionBegin;
621a4d4d686SBarry Smith 
622a4d4d686SBarry Smith   ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr);
623a4d4d686SBarry Smith   /* no context indicates that it is not the "matrix free" matrix type */
624a4d4d686SBarry Smith   if (!ctx) PetscFunctionReturn(0);
625a4d4d686SBarry Smith   ctx->currenth    = 0;
626a4d4d686SBarry Smith 
627a4d4d686SBarry Smith   PetscFunctionReturn(0);
628a4d4d686SBarry Smith }
629a4d4d686SBarry Smith 
630