1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.7 1995/05/14 16:35:00 bsmith Exp curfman $"; 4 #endif 5 6 #include "draw.h" 7 #include "snes.h" 8 9 typedef struct { 10 SNES snes; 11 Vec w; 12 } MFCtx_Private; 13 14 /* 15 SNESMatrixFreeMult_Private - Default matrix free form of A*u. 16 17 */ 18 int SNESMatrixFreeMult_Private(void *ptr,Vec dx,Vec y) 19 { 20 MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 21 SNES snes = ctx->snes; 22 double norm,epsilon = 1.e-8; /* assumes double precision */ 23 Scalar h,dot,sum; 24 Scalar mone = -1.0; 25 Vec w = ctx->w,U,F; 26 int ierr; 27 28 ierr = SNESGetSolution(snes,&U); CHKERR(ierr); 29 ierr = SNESGetFunction(snes,&F); CHKERR(ierr); 30 /* determine a "good" step size */ 31 VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 32 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 33 else 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 MatStructure *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