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