181e6777dSBarry Smith 281e6777dSBarry Smith #ifndef lint 3*39e2f89bSBarry Smith static char vcid[] = "$Id: snesmfj.c,v 1.1 1995/05/04 03:17:06 bsmith Exp bsmith $"; 481e6777dSBarry Smith #endif 581e6777dSBarry Smith 681e6777dSBarry Smith #include "draw.h" 781e6777dSBarry Smith #include "snesimpl.h" 881e6777dSBarry Smith #include "options.h" 981e6777dSBarry Smith 10*39e2f89bSBarry Smith typedef struct { 11*39e2f89bSBarry Smith SNES snes; 12*39e2f89bSBarry Smith Vec w; 13*39e2f89bSBarry Smith } MFCtx_Private; 14*39e2f89bSBarry Smith 15*39e2f89bSBarry Smith /* 16*39e2f89bSBarry Smith SNESMatrixFreeMult_Private - Default matrix free form of A*u. 17*39e2f89bSBarry Smith 18*39e2f89bSBarry Smith */ 19*39e2f89bSBarry Smith int SNESMatrixFreeMult_Private(void *ptr,Vec dx,Vec y) 20*39e2f89bSBarry Smith { 21*39e2f89bSBarry Smith MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 22*39e2f89bSBarry Smith SNES snes = ctx->snes; 23*39e2f89bSBarry Smith double norm,epsilon = 1.e-8; /* assumes double precision */ 24*39e2f89bSBarry Smith Scalar h,dot,sum; 25*39e2f89bSBarry Smith Scalar mone = -1.0; 26*39e2f89bSBarry Smith Vec w = ctx->w,U,F; 27*39e2f89bSBarry Smith int ierr; 28*39e2f89bSBarry Smith 29*39e2f89bSBarry Smith SNESGetSolution(snes,&U); SNESGetFunction(snes,&F); 30*39e2f89bSBarry Smith /* determine a "good" step size */ 31*39e2f89bSBarry Smith VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 32*39e2f89bSBarry Smith if (dot < 1.e-16*sum && dot >= 0.0) dot = 1.e-16*sum; 33*39e2f89bSBarry Smith else if (dot < 0.0 && dot > 1.e-16*sum) dot = -1.e-16*sum; 34*39e2f89bSBarry Smith h = epsilon*dot/(norm*norm); 35*39e2f89bSBarry Smith 36*39e2f89bSBarry Smith /* evaluate function at F(x + dx) */ 37*39e2f89bSBarry Smith VecWAXPY(&h,dx,U,w); 38*39e2f89bSBarry Smith ierr = SNESComputeFunction(snes,w,y); 39*39e2f89bSBarry Smith VecAXPY(&mone,F,y); 40*39e2f89bSBarry Smith h = -1.0/h; 41*39e2f89bSBarry Smith VecScale(&h,y); 42*39e2f89bSBarry Smith return 0; 43*39e2f89bSBarry Smith } 4481e6777dSBarry Smith /*@ 45*39e2f89bSBarry Smith SNESDefaultMatrixFreeComputeJacobian - Computes Jacobian using finite 46*39e2f89bSBarry Smith differences, matrix free style. 4781e6777dSBarry Smith 4881e6777dSBarry Smith Input Parameters: 4981e6777dSBarry Smith . x - compute Jacobian at this point 5081e6777dSBarry Smith . ctx - applications Function context 5181e6777dSBarry Smith 5281e6777dSBarry Smith Output Parameters: 5381e6777dSBarry Smith . J - Jacobian 5481e6777dSBarry Smith . B - preconditioner, same as Jacobian 5581e6777dSBarry Smith 5681e6777dSBarry Smith .keywords: finite differences, Jacobian 5781e6777dSBarry Smith 5881e6777dSBarry Smith .seealso: SNESSetJacobian, SNESTestJacobian 5981e6777dSBarry Smith @*/ 60*39e2f89bSBarry Smith int SNESDefaultMatrixFreeComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B, 61*39e2f89bSBarry Smith int *flag,void *ctx) 6281e6777dSBarry Smith { 63*39e2f89bSBarry Smith int n,ierr; 64*39e2f89bSBarry Smith MPI_Comm comm; 65*39e2f89bSBarry Smith PetscObject obj = (PetscObject) *J; 66*39e2f89bSBarry Smith VecGetSize(x1,&n); 67*39e2f89bSBarry Smith comm = ((PetscObject)x1)->comm; 68*39e2f89bSBarry Smith if (!*J || obj->type != MATSHELL) { 69*39e2f89bSBarry Smith MFCtx_Private *mfctx; 70*39e2f89bSBarry Smith /* first time in, therefore build datastructures */ 71*39e2f89bSBarry Smith mfctx = NEW(MFCtx_Private); CHKPTR(mfctx); 72*39e2f89bSBarry Smith mfctx->snes = snes; 73*39e2f89bSBarry Smith ierr = VecDuplicate(x1,&mfctx->w); CHKERR(ierr); 74*39e2f89bSBarry Smith ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERR(ierr); 75*39e2f89bSBarry Smith MatShellSetMult(*J,SNESMatrixFreeMult_Private); 76*39e2f89bSBarry Smith *B = *J; 7781e6777dSBarry Smith } 7881e6777dSBarry Smith return 0; 7981e6777dSBarry Smith } 8081e6777dSBarry Smith 81