181e6777dSBarry Smith 281e6777dSBarry Smith #ifndef lint 3*ad960d00SLois Curfman McInnes static char vcid[] = "$Id: snesmfj.c,v 1.5 1995/05/11 17:20:55 curfman Exp curfman $"; 481e6777dSBarry Smith #endif 581e6777dSBarry Smith 681e6777dSBarry Smith #include "draw.h" 7f3b845d0SBarry Smith #include "snes.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 29*ad960d00SLois Curfman McInnes SNESGetSolution(snes,&U); CHKERR(ierr); 30*ad960d00SLois Curfman McInnes SNESGetFunction(snes,&F); CHKERR(ierr); 3139e2f89bSBarry Smith /* determine a "good" step size */ 3239e2f89bSBarry Smith VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 3339e2f89bSBarry Smith if (dot < 1.e-16*sum && dot >= 0.0) dot = 1.e-16*sum; 3439e2f89bSBarry Smith else if (dot < 0.0 && dot > 1.e-16*sum) dot = -1.e-16*sum; 3539e2f89bSBarry Smith h = epsilon*dot/(norm*norm); 3639e2f89bSBarry Smith 3739e2f89bSBarry Smith /* evaluate function at F(x + dx) */ 3839e2f89bSBarry Smith VecWAXPY(&h,dx,U,w); 39*ad960d00SLois Curfman McInnes ierr = SNESComputeFunction(snes,w,y); CHKERR(ierr); 4039e2f89bSBarry Smith VecAXPY(&mone,F,y); 4139e2f89bSBarry Smith h = -1.0/h; 4239e2f89bSBarry Smith VecScale(&h,y); 4339e2f89bSBarry Smith return 0; 4439e2f89bSBarry Smith } 4581e6777dSBarry Smith /*@ 4639e2f89bSBarry Smith SNESDefaultMatrixFreeComputeJacobian - Computes Jacobian using finite 4771dce806SLois Curfman McInnes differences, matrix-free style. 4881e6777dSBarry Smith 4981e6777dSBarry Smith Input Parameters: 5081e6777dSBarry Smith . x - compute Jacobian at this point 51*ad960d00SLois Curfman McInnes . ctx - application's function context, as set with SNESSetFunction() 5281e6777dSBarry Smith 5381e6777dSBarry Smith Output Parameters: 5481e6777dSBarry Smith . J - Jacobian 5581e6777dSBarry Smith . B - preconditioner, same as Jacobian 56*ad960d00SLois Curfman McInnes . flag - matrix flag 57*ad960d00SLois Curfman McInnes 58*ad960d00SLois Curfman McInnes Options Database Key: 59*ad960d00SLois Curfman McInnes $ -snes_mf 6081e6777dSBarry Smith 6171dce806SLois Curfman McInnes .keywords: SNES, finite differences, Jacobian 6281e6777dSBarry Smith 6371dce806SLois Curfman McInnes .seealso: SNESSetJacobian(), SNESTestJacobian() 6481e6777dSBarry Smith @*/ 6539e2f89bSBarry Smith int SNESDefaultMatrixFreeComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B, 6639e2f89bSBarry Smith int *flag,void *ctx) 6781e6777dSBarry Smith { 6839e2f89bSBarry Smith int n,ierr; 697f223b93SBarry Smith MatType type; 707f223b93SBarry Smith 7139e2f89bSBarry Smith VecGetSize(x1,&n); 727f223b93SBarry Smith if (*J) MatGetType(*J,&type); 737f223b93SBarry Smith if (!*J || type != MATSHELL) { 747f223b93SBarry Smith MPI_Comm comm; 7539e2f89bSBarry Smith MFCtx_Private *mfctx; 7639e2f89bSBarry Smith /* first time in, therefore build datastructures */ 7739e2f89bSBarry Smith mfctx = NEW(MFCtx_Private); CHKPTR(mfctx); 7839e2f89bSBarry Smith mfctx->snes = snes; 7939e2f89bSBarry Smith ierr = VecDuplicate(x1,&mfctx->w); CHKERR(ierr); 807f223b93SBarry Smith PetscObjectGetComm((PetscObject)x1,&comm); 8139e2f89bSBarry Smith ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERR(ierr); 8239e2f89bSBarry Smith MatShellSetMult(*J,SNESMatrixFreeMult_Private); 8339e2f89bSBarry Smith *B = *J; 8481e6777dSBarry Smith } 8581e6777dSBarry Smith return 0; 8681e6777dSBarry Smith } 8781e6777dSBarry Smith 88