1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.9 1995/06/08 03:11:42 bsmith Exp bsmith $"; 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); CHKERRQ(ierr); 29 ierr = SNESGetFunction(snes,&F); CHKERRQ(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); CHKERRQ(ierr); 40 VecAXPY(&mone,F,y); 41 h = -1.0/h; 42 VecScale(&h,y); 43 return 0; 44 } 45 /*@ 46 SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 47 for use with SNES solver. You may use this matrix as 48 Jacobian argument for the routine SNESSetJacobian. This is 49 most useful when you are using finite differences for a 50 matrix free Newton method but explictly are forming a 51 preconditioner matrix. 52 53 Input Parameters: 54 . x - vector where SNES solution is to be stored. 55 56 Output Parameters: 57 . J - the matrix-free matrix 58 59 @*/ 60 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 61 { 62 MPI_Comm comm; 63 MFCtx_Private *mfctx; 64 int n,ierr; 65 66 mfctx = (MFCtx_Private *) PETSCMALLOC(sizeof(MFCtx_Private));CHKPTRQ(mfctx); 67 mfctx->snes = snes; 68 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 69 PetscObjectGetComm((PetscObject)x,&comm); 70 VecGetSize(x,&n); 71 ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERRQ(ierr); 72 MatShellSetMult(*J,SNESMatrixFreeMult_Private); 73 return 0; 74 } 75 /*@ 76 SNESDefaultMatrixFreeComputeJacobian - Computes Jacobian using finite 77 differences, matrix-free style. 78 79 Input Parameters: 80 . x - compute Jacobian at this point 81 . ctx - application's function context, as set with SNESSetFunction() 82 83 Output Parameters: 84 . J - Jacobian 85 . B - preconditioner, same as Jacobian 86 . flag - matrix flag 87 88 Options Database Key: 89 $ -snes_mf 90 91 .keywords: SNES, finite differences, Jacobian 92 93 .seealso: SNESSetJacobian(), SNESTestJacobian() 94 @*/ 95 int SNESDefaultMatrixFreeComputeJacobian(SNES snes, Vec x1,Mat *J,Mat *B, 96 MatStructure *flag,void *ctx) 97 { 98 int n,ierr; 99 MatType type; 100 101 VecGetSize(x1,&n); 102 if (*J) MatGetType(*J,&type); 103 if (!*J || type != MATSHELL) { 104 MPI_Comm comm; 105 MFCtx_Private *mfctx; 106 /* first time in, therefore build datastructures */ 107 mfctx = (MFCtx_Private *) PETSCMALLOC(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 108 mfctx->snes = snes; 109 ierr = VecDuplicate(x1,&mfctx->w); CHKERRQ(ierr); 110 PetscObjectGetComm((PetscObject)x1,&comm); 111 ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERRQ(ierr); 112 MatShellSetMult(*J,SNESMatrixFreeMult_Private); 113 *B = *J; 114 } 115 return 0; 116 } 117 118