181e6777dSBarry Smith 281e6777dSBarry Smith #ifndef lint 3*7f223b93SBarry Smith static char vcid[] = "$Id: snesmfj.c,v 1.2 1995/05/05 03:51:29 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 1039e2f89bSBarry Smith typedef struct { 1139e2f89bSBarry Smith SNES snes; 1239e2f89bSBarry Smith Vec w; 1339e2f89bSBarry Smith } MFCtx_Private; 1439e2f89bSBarry Smith 1539e2f89bSBarry Smith /* 1639e2f89bSBarry Smith SNESMatrixFreeMult_Private - Default matrix free form of A*u. 1739e2f89bSBarry Smith 1839e2f89bSBarry Smith */ 1939e2f89bSBarry Smith int SNESMatrixFreeMult_Private(void *ptr,Vec dx,Vec y) 2039e2f89bSBarry Smith { 2139e2f89bSBarry Smith MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 2239e2f89bSBarry Smith SNES snes = ctx->snes; 2339e2f89bSBarry Smith double norm,epsilon = 1.e-8; /* assumes double precision */ 2439e2f89bSBarry Smith Scalar h,dot,sum; 2539e2f89bSBarry Smith Scalar mone = -1.0; 2639e2f89bSBarry Smith Vec w = ctx->w,U,F; 2739e2f89bSBarry Smith int ierr; 2839e2f89bSBarry Smith 2939e2f89bSBarry Smith SNESGetSolution(snes,&U); SNESGetFunction(snes,&F); 3039e2f89bSBarry Smith /* determine a "good" step size */ 3139e2f89bSBarry Smith VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 3239e2f89bSBarry Smith if (dot < 1.e-16*sum && dot >= 0.0) dot = 1.e-16*sum; 3339e2f89bSBarry Smith else if (dot < 0.0 && dot > 1.e-16*sum) dot = -1.e-16*sum; 3439e2f89bSBarry Smith h = epsilon*dot/(norm*norm); 3539e2f89bSBarry Smith 3639e2f89bSBarry Smith /* evaluate function at F(x + dx) */ 3739e2f89bSBarry Smith VecWAXPY(&h,dx,U,w); 3839e2f89bSBarry Smith ierr = SNESComputeFunction(snes,w,y); 3939e2f89bSBarry Smith VecAXPY(&mone,F,y); 4039e2f89bSBarry Smith h = -1.0/h; 4139e2f89bSBarry Smith VecScale(&h,y); 4239e2f89bSBarry Smith return 0; 4339e2f89bSBarry Smith } 4481e6777dSBarry Smith /*@ 4539e2f89bSBarry Smith SNESDefaultMatrixFreeComputeJacobian - Computes Jacobian using finite 4639e2f89bSBarry 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 @*/ 6039e2f89bSBarry Smith int SNESDefaultMatrixFreeComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B, 6139e2f89bSBarry Smith int *flag,void *ctx) 6281e6777dSBarry Smith { 6339e2f89bSBarry Smith int n,ierr; 64*7f223b93SBarry Smith MatType type; 65*7f223b93SBarry Smith 6639e2f89bSBarry Smith VecGetSize(x1,&n); 67*7f223b93SBarry Smith if (*J) MatGetType(*J,&type); 68*7f223b93SBarry Smith if (!*J || type != MATSHELL) { 69*7f223b93SBarry Smith MPI_Comm comm; 7039e2f89bSBarry Smith MFCtx_Private *mfctx; 7139e2f89bSBarry Smith /* first time in, therefore build datastructures */ 7239e2f89bSBarry Smith mfctx = NEW(MFCtx_Private); CHKPTR(mfctx); 7339e2f89bSBarry Smith mfctx->snes = snes; 7439e2f89bSBarry Smith ierr = VecDuplicate(x1,&mfctx->w); CHKERR(ierr); 75*7f223b93SBarry Smith PetscObjectGetComm((PetscObject)x1,&comm); 7639e2f89bSBarry Smith ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERR(ierr); 7739e2f89bSBarry Smith MatShellSetMult(*J,SNESMatrixFreeMult_Private); 7839e2f89bSBarry Smith *B = *J; 7981e6777dSBarry Smith } 8081e6777dSBarry Smith return 0; 8181e6777dSBarry Smith } 8281e6777dSBarry Smith 83