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