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