1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.4 1995/05/05 11:48:16 bsmith 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); SNESGetFunction(snes,&F); 30 /* determine a "good" step size */ 31 VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 32 if (dot < 1.e-16*sum && dot >= 0.0) dot = 1.e-16*sum; 33 else if (dot < 0.0 && dot > 1.e-16*sum) dot = -1.e-16*sum; 34 h = epsilon*dot/(norm*norm); 35 36 /* evaluate function at F(x + dx) */ 37 VecWAXPY(&h,dx,U,w); 38 ierr = SNESComputeFunction(snes,w,y); 39 VecAXPY(&mone,F,y); 40 h = -1.0/h; 41 VecScale(&h,y); 42 return 0; 43 } 44 /*@ 45 SNESDefaultMatrixFreeComputeJacobian - Computes Jacobian using finite 46 differences, matrix-free style. 47 48 Input Parameters: 49 . x - compute Jacobian at this point 50 . ctx - applications Function context 51 52 Output Parameters: 53 . J - Jacobian 54 . B - preconditioner, same as Jacobian 55 56 .keywords: SNES, finite differences, Jacobian 57 58 .seealso: SNESSetJacobian(), SNESTestJacobian() 59 @*/ 60 int SNESDefaultMatrixFreeComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B, 61 int *flag,void *ctx) 62 { 63 int n,ierr; 64 MatType type; 65 66 VecGetSize(x1,&n); 67 if (*J) MatGetType(*J,&type); 68 if (!*J || type != MATSHELL) { 69 MPI_Comm comm; 70 MFCtx_Private *mfctx; 71 /* first time in, therefore build datastructures */ 72 mfctx = NEW(MFCtx_Private); CHKPTR(mfctx); 73 mfctx->snes = snes; 74 ierr = VecDuplicate(x1,&mfctx->w); CHKERR(ierr); 75 PetscObjectGetComm((PetscObject)x1,&comm); 76 ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERR(ierr); 77 MatShellSetMult(*J,SNESMatrixFreeMult_Private); 78 *B = *J; 79 } 80 return 0; 81 } 82 83