181e6777dSBarry Smith 281e6777dSBarry Smith #ifndef lint 3*edd2f0e1SBarry Smith static char vcid[] = "$Id: snesmfj.c,v 1.6 1995/05/12 18:18:46 curfman Exp bsmith $"; 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*edd2f0e1SBarry Smith ierr = SNESGetSolution(snes,&U); CHKERR(ierr); 30*edd2f0e1SBarry Smith ierr = SNESGetFunction(snes,&F); CHKERR(ierr); 3139e2f89bSBarry Smith /* determine a "good" step size */ 3239e2f89bSBarry Smith VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 33*edd2f0e1SBarry Smith if (sum == 0.0) {dot = 1.0; norm = 1.0;} 34*edd2f0e1SBarry Smith else if (dot < 1.e-16*sum && dot >= 0.0) dot = 1.e-16*sum; 3539e2f89bSBarry Smith else if (dot < 0.0 && dot > 1.e-16*sum) dot = -1.e-16*sum; 3639e2f89bSBarry Smith h = epsilon*dot/(norm*norm); 3739e2f89bSBarry Smith 3839e2f89bSBarry Smith /* evaluate function at F(x + dx) */ 3939e2f89bSBarry Smith VecWAXPY(&h,dx,U,w); 40ad960d00SLois Curfman McInnes ierr = SNESComputeFunction(snes,w,y); CHKERR(ierr); 4139e2f89bSBarry Smith VecAXPY(&mone,F,y); 4239e2f89bSBarry Smith h = -1.0/h; 4339e2f89bSBarry Smith VecScale(&h,y); 4439e2f89bSBarry Smith return 0; 4539e2f89bSBarry Smith } 4681e6777dSBarry Smith /*@ 4739e2f89bSBarry Smith SNESDefaultMatrixFreeComputeJacobian - Computes Jacobian using finite 4871dce806SLois Curfman McInnes differences, matrix-free style. 4981e6777dSBarry Smith 5081e6777dSBarry Smith Input Parameters: 5181e6777dSBarry Smith . x - compute Jacobian at this point 52ad960d00SLois Curfman McInnes . ctx - application's function context, as set with SNESSetFunction() 5381e6777dSBarry Smith 5481e6777dSBarry Smith Output Parameters: 5581e6777dSBarry Smith . J - Jacobian 5681e6777dSBarry Smith . B - preconditioner, same as Jacobian 57ad960d00SLois Curfman McInnes . flag - matrix flag 58ad960d00SLois Curfman McInnes 59ad960d00SLois Curfman McInnes Options Database Key: 60ad960d00SLois Curfman McInnes $ -snes_mf 6181e6777dSBarry Smith 6271dce806SLois Curfman McInnes .keywords: SNES, finite differences, Jacobian 6381e6777dSBarry Smith 6471dce806SLois Curfman McInnes .seealso: SNESSetJacobian(), SNESTestJacobian() 6581e6777dSBarry Smith @*/ 6639e2f89bSBarry Smith int SNESDefaultMatrixFreeComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B, 67*edd2f0e1SBarry Smith MatStructure *flag,void *ctx) 6881e6777dSBarry Smith { 6939e2f89bSBarry Smith int n,ierr; 707f223b93SBarry Smith MatType type; 717f223b93SBarry Smith 7239e2f89bSBarry Smith VecGetSize(x1,&n); 737f223b93SBarry Smith if (*J) MatGetType(*J,&type); 747f223b93SBarry Smith if (!*J || type != MATSHELL) { 757f223b93SBarry Smith MPI_Comm comm; 7639e2f89bSBarry Smith MFCtx_Private *mfctx; 7739e2f89bSBarry Smith /* first time in, therefore build datastructures */ 7839e2f89bSBarry Smith mfctx = NEW(MFCtx_Private); CHKPTR(mfctx); 7939e2f89bSBarry Smith mfctx->snes = snes; 8039e2f89bSBarry Smith ierr = VecDuplicate(x1,&mfctx->w); CHKERR(ierr); 817f223b93SBarry Smith PetscObjectGetComm((PetscObject)x1,&comm); 8239e2f89bSBarry Smith ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERR(ierr); 8339e2f89bSBarry Smith MatShellSetMult(*J,SNESMatrixFreeMult_Private); 8439e2f89bSBarry Smith *B = *J; 8581e6777dSBarry Smith } 8681e6777dSBarry Smith return 0; 8781e6777dSBarry Smith } 8881e6777dSBarry Smith 89