1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.23 1995/12/21 18:33:54 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 { /* default context for matrix-free SNES */ 10 SNES snes; 11 Vec w; 12 PCNullSpace sp; 13 } MFCtx_Private; 14 15 int SNESMatrixFreeDestroy_Private(void *ptr) 16 { 17 int ierr; 18 MFCtx_Private *ctx = (MFCtx_Private* ) ptr; 19 ierr = VecDestroy(ctx->w); CHKERRQ(ierr); 20 if (ctx->sp) {ierr = PCNullSpaceDestroy(ctx->sp); CHKERRQ(ierr);} 21 PetscFree(ptr); 22 return 0; 23 } 24 25 /* 26 SNESMatrixFreeMult_Private - Default matrix free form of A*u. 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,sum,epsilon = 1.e-8; /* assumes double precision */ 33 Scalar h,dot,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 = VecNorm(dx,NORM_1,&sum); CHKERRQ(ierr); 43 ierr = VecNorm(dx,NORM_2,&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 if (ctx->sp) { ierr = PCNullSpaceRemove(ctx->sp,y); CHKERRQ(ierr);} 61 return 0; 62 } 63 /*@C 64 SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 65 context for use with a SNES solver. This matrix can be used as 66 the Jacobian argument for the routine SNESSetJacobian(). 67 68 Input Parameters: 69 . x - vector where SNES solution is to be stored. 70 71 Output Parameters: 72 . J - the matrix-free matrix 73 74 Notes: 75 The matrix-free matrix context merely contains the function pointers 76 and work space for performing finite difference approximations of 77 matrix operations such as matrix-vector products. 78 79 The user should call MatDestroy() when finished with the matrix-free 80 matrix context. 81 82 .keywords: SNES, default, matrix-free, create, matrix 83 84 .seealso: MatDestroy() 85 @*/ 86 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 87 { 88 MPI_Comm comm; 89 MFCtx_Private *mfctx; 90 int n,ierr; 91 92 mfctx = (MFCtx_Private *) PetscMalloc(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 93 PLogObjectMemory(snes,sizeof(MFCtx_Private)); 94 mfctx->sp = 0; 95 mfctx->snes = snes; 96 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 97 ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr); 98 ierr = VecGetSize(x,&n); CHKERRQ(ierr); 99 ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERRQ(ierr); 100 ierr = MatShellSetMult(*J,SNESMatrixFreeMult_Private); CHKERRQ(ierr); 101 ierr = MatShellSetDestroy(*J,SNESMatrixFreeDestroy_Private); CHKERRQ(ierr); 102 PLogObjectParent(*J,mfctx->w); 103 PLogObjectParent(snes,*J); 104 return 0; 105 } 106 107 /*@ 108 SNESDefaultMatrixFreeMatAddNullSpace - Provide a null space that 109 an operator is suppose to have. Since round off will create a 110 small component in the null space, if you know the null space 111 you may have it automatically removed. 112 113 Input Parameters: 114 . J - the matrix free matrix 115 . has_cnst - PETSC_TRUE or PETSC_FALSE indicating if null space has constants 116 . n - number of vectors (excluding constant vector) in null space 117 . vecs - the vectors that span the null space (excluding the constant vector) 118 . these vectors must be orthonormal 119 120 .keywords: SNES, matrix-free, null space 121 @*/ 122 int SNESDefaultMatrixFreeMatAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs) 123 { 124 int ierr; 125 MFCtx_Private *ctx; 126 MPI_Comm comm; 127 128 PetscObjectGetComm((PetscObject)J,&comm); 129 130 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 131 /* no context indicates that it is not the "matrix free" matrix type */ 132 if (!ctx) return 0; 133 ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr); 134 return 0; 135 } 136 137 138 139 140