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