1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.12 1995/06/29 23:54:14 bsmith Exp bsmith $"; 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 for use with SNES solver. You may use this matrix as 64 Jacobian argument for the routine SNESSetJacobian. This is 65 most useful when you are using finite differences for a 66 matrix free Newton method but explictly are forming a 67 preconditioner matrix. 68 69 Input Parameters: 70 . x - vector where SNES solution is to be stored. 71 72 Output Parameters: 73 . J - the matrix-free matrix 74 75 @*/ 76 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 77 { 78 MPI_Comm comm; 79 MFCtx_Private *mfctx; 80 int n,ierr; 81 82 mfctx = (MFCtx_Private *) PETSCMALLOC(sizeof(MFCtx_Private));CHKPTRQ(mfctx); 83 mfctx->snes = snes; 84 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 85 PetscObjectGetComm((PetscObject)x,&comm); 86 VecGetSize(x,&n); 87 ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERRQ(ierr); 88 MatShellSetMult(*J,SNESMatrixFreeMult_Private); 89 MatShellSetDestroy(*J,SNESMatrixFreeDestroy_Private); 90 PLogObjectParent(*J,mfctx->w); 91 return 0; 92 } 93 94