1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.13 1995/08/02 04:18:56 bsmith Exp curfman $"; 4 #endif 5 6 #include "draw.h" /*I "draw.h" I*/ 7 #include "snes.h" /*I "snes.h" I*/ 8 #include "ptscimpl.h" 9 #include "plog.h" 10 11 typedef struct { 12 SNES snes; 13 Vec w; 14 } MFCtx_Private; 15 16 int SNESMatrixFreeDestroy_Private(void *ptr) 17 { 18 int ierr; 19 MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 20 ierr = VecDestroy(ctx->w); CHKERRQ(ierr); 21 PETSCFREE(ptr); 22 return 0; 23 } 24 /* 25 SNESMatrixFreeMult_Private - Default matrix free form of A*u. 26 27 */ 28 int SNESMatrixFreeMult_Private(void *ptr,Vec dx,Vec y) 29 { 30 MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 31 SNES snes = ctx->snes; 32 double norm,epsilon = 1.e-8; /* assumes double precision */ 33 Scalar h,dot; 34 double sum; 35 Scalar mone = -1.0; 36 Vec w = ctx->w,U,F; 37 int ierr; 38 39 ierr = SNESGetSolution(snes,&U); CHKERRQ(ierr); 40 ierr = SNESGetFunction(snes,&F); CHKERRQ(ierr); 41 /* determine a "good" step size */ 42 VecDot(U,dx,&dot); VecASum(dx,&sum); VecNorm(dx,&norm); 43 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 44 #if defined(PETSC_COMPLEX) 45 else if (abs(dot) < 1.e-16*sum && real(dot) >= 0.0) dot = 1.e-16*sum; 46 else if (abs(dot) < 0.0 && real(dot) > 1.e-16*sum) dot = -1.e-16*sum; 47 #else 48 else if (dot < 1.e-16*sum && dot >= 0.0) dot = 1.e-16*sum; 49 else if (dot < 0.0 && dot > 1.e-16*sum) dot = -1.e-16*sum; 50 #endif 51 h = epsilon*dot/(norm*norm); 52 53 /* evaluate function at F(x + dx) */ 54 VecWAXPY(&h,dx,U,w); 55 ierr = SNESComputeFunction(snes,w,y); CHKERRQ(ierr); 56 VecAXPY(&mone,F,y); 57 h = -1.0/h; 58 VecScale(&h,y); 59 return 0; 60 } 61 /*@ 62 SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 63 context for use with a SNES solver. You can use this matrix as the 64 Jacobian argument for the routine SNESSetJacobian(). 65 66 Input Parameters: 67 . x - vector where SNES solution is to be stored. 68 69 Output Parameters: 70 . J - the matrix-free matrix 71 72 Notes: 73 The matrix-free matrix context merely contains the function pointers 74 and work space for performing finite difference approximations of 75 matrix operations such as matrix-vector products. 76 77 The user should call MatDestroy() when finished with the matrix-free 78 matrix context. 79 80 .keywords: SNES, default, matrix-free, create, matrix 81 82 .seealso: MatDestroy() 83 @*/ 84 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 85 { 86 MPI_Comm comm; 87 MFCtx_Private *mfctx; 88 int n,ierr; 89 90 mfctx = (MFCtx_Private *) PETSCMALLOC(sizeof(MFCtx_Private));CHKPTRQ(mfctx); 91 mfctx->snes = snes; 92 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 93 PetscObjectGetComm((PetscObject)x,&comm); 94 VecGetSize(x,&n); 95 ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERRQ(ierr); 96 MatShellSetMult(*J,SNESMatrixFreeMult_Private); 97 MatShellSetDestroy(*J,SNESMatrixFreeDestroy_Private); 98 PLogObjectParent(*J,mfctx->w); 99 return 0; 100 } 101 102