181e6777dSBarry Smith 281e6777dSBarry Smith #ifndef lint 3*78b31e54SBarry Smith static char vcid[] = "$Id: snesmfj.c,v 1.8 1995/05/16 00:37:02 curfman Exp bsmith $"; 481e6777dSBarry Smith #endif 581e6777dSBarry Smith 681e6777dSBarry Smith #include "draw.h" 7f3b845d0SBarry Smith #include "snes.h" 881e6777dSBarry Smith 939e2f89bSBarry Smith typedef struct { 1039e2f89bSBarry Smith SNES snes; 1139e2f89bSBarry Smith Vec w; 1239e2f89bSBarry Smith } MFCtx_Private; 1339e2f89bSBarry Smith 1439e2f89bSBarry Smith /* 1539e2f89bSBarry Smith SNESMatrixFreeMult_Private - Default matrix free form of A*u. 1639e2f89bSBarry Smith 1739e2f89bSBarry Smith */ 1839e2f89bSBarry Smith int SNESMatrixFreeMult_Private(void *ptr,Vec dx,Vec y) 1939e2f89bSBarry Smith { 2039e2f89bSBarry Smith MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 2139e2f89bSBarry Smith SNES snes = ctx->snes; 2239e2f89bSBarry Smith double norm,epsilon = 1.e-8; /* assumes double precision */ 2339e2f89bSBarry Smith Scalar h,dot,sum; 2439e2f89bSBarry Smith Scalar mone = -1.0; 2539e2f89bSBarry Smith Vec w = ctx->w,U,F; 2639e2f89bSBarry Smith int ierr; 2739e2f89bSBarry Smith 28*78b31e54SBarry Smith ierr = SNESGetSolution(snes,&U); CHKERRQ(ierr); 29*78b31e54SBarry Smith ierr = SNESGetFunction(snes,&F); CHKERRQ(ierr); 3039e2f89bSBarry Smith /* determine a "good" step size */ 3139e2f89bSBarry Smith VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 32edd2f0e1SBarry Smith if (sum == 0.0) {dot = 1.0; norm = 1.0;} 33edd2f0e1SBarry Smith else 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*78b31e54SBarry Smith ierr = SNESComputeFunction(snes,w,y); CHKERRQ(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 51ad960d00SLois 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 56ad960d00SLois Curfman McInnes . flag - matrix flag 57ad960d00SLois Curfman McInnes 58ad960d00SLois Curfman McInnes Options Database Key: 59ad960d00SLois 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, 66edd2f0e1SBarry Smith MatStructure *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 */ 77*78b31e54SBarry Smith mfctx = (MFCtx_Private *) PETSCMALLOC(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 7839e2f89bSBarry Smith mfctx->snes = snes; 79*78b31e54SBarry Smith ierr = VecDuplicate(x1,&mfctx->w); CHKERRQ(ierr); 807f223b93SBarry Smith PetscObjectGetComm((PetscObject)x1,&comm); 81*78b31e54SBarry Smith ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERRQ(ierr); 8239e2f89bSBarry Smith MatShellSetMult(*J,SNESMatrixFreeMult_Private); 8339e2f89bSBarry Smith *B = *J; 8481e6777dSBarry Smith } 8581e6777dSBarry Smith return 0; 8681e6777dSBarry Smith } 8781e6777dSBarry Smith 88