1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.5 1995/05/11 17:20:55 curfman Exp curfman $"; 4 #endif 5 6 #include "draw.h" 7 #include "snes.h" 8 #include "options.h" 9 10 typedef struct { 11 SNES snes; 12 Vec w; 13 } MFCtx_Private; 14 15 /* 16 SNESMatrixFreeMult_Private - Default matrix free form of A*u. 17 18 */ 19 int SNESMatrixFreeMult_Private(void *ptr,Vec dx,Vec y) 20 { 21 MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 22 SNES snes = ctx->snes; 23 double norm,epsilon = 1.e-8; /* assumes double precision */ 24 Scalar h,dot,sum; 25 Scalar mone = -1.0; 26 Vec w = ctx->w,U,F; 27 int ierr; 28 29 SNESGetSolution(snes,&U); CHKERR(ierr); 30 SNESGetFunction(snes,&F); CHKERR(ierr); 31 /* determine a "good" step size */ 32 VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 33 if (dot < 1.e-16*sum && dot >= 0.0) dot = 1.e-16*sum; 34 else if (dot < 0.0 && dot > 1.e-16*sum) dot = -1.e-16*sum; 35 h = epsilon*dot/(norm*norm); 36 37 /* evaluate function at F(x + dx) */ 38 VecWAXPY(&h,dx,U,w); 39 ierr = SNESComputeFunction(snes,w,y); CHKERR(ierr); 40 VecAXPY(&mone,F,y); 41 h = -1.0/h; 42 VecScale(&h,y); 43 return 0; 44 } 45 /*@ 46 SNESDefaultMatrixFreeComputeJacobian - Computes Jacobian using finite 47 differences, matrix-free style. 48 49 Input Parameters: 50 . x - compute Jacobian at this point 51 . ctx - application's function context, as set with SNESSetFunction() 52 53 Output Parameters: 54 . J - Jacobian 55 . B - preconditioner, same as Jacobian 56 . flag - matrix flag 57 58 Options Database Key: 59 $ -snes_mf 60 61 .keywords: SNES, finite differences, Jacobian 62 63 .seealso: SNESSetJacobian(), SNESTestJacobian() 64 @*/ 65 int SNESDefaultMatrixFreeComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B, 66 int *flag,void *ctx) 67 { 68 int n,ierr; 69 MatType type; 70 71 VecGetSize(x1,&n); 72 if (*J) MatGetType(*J,&type); 73 if (!*J || type != MATSHELL) { 74 MPI_Comm comm; 75 MFCtx_Private *mfctx; 76 /* first time in, therefore build datastructures */ 77 mfctx = NEW(MFCtx_Private); CHKPTR(mfctx); 78 mfctx->snes = snes; 79 ierr = VecDuplicate(x1,&mfctx->w); CHKERR(ierr); 80 PetscObjectGetComm((PetscObject)x1,&comm); 81 ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERR(ierr); 82 MatShellSetMult(*J,SNESMatrixFreeMult_Private); 83 *B = *J; 84 } 85 return 0; 86 } 87 88