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