1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.22 1995/11/01 23:20:55 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,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 42 /* Determine a "good" step size */ 43 ierr = VecDot(U,dx,&dot); CHKERRQ(ierr); 44 ierr = VecNorm(dx,NORM_1,&sum); CHKERRQ(ierr); 45 ierr = VecNorm(dx,NORM_2,&norm); CHKERRQ(ierr); 46 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 47 #if defined(PETSC_COMPLEX) 48 else if (abs(dot) < 1.e-16*sum && real(dot) >= 0.0) dot = 1.e-16*sum; 49 else if (abs(dot) < 0.0 && real(dot) > 1.e-16*sum) dot = -1.e-16*sum; 50 #else 51 else if (dot < 1.e-16*sum && dot >= 0.0) dot = 1.e-16*sum; 52 else if (dot < 0.0 && dot > 1.e-16*sum) dot = -1.e-16*sum; 53 #endif 54 h = epsilon*dot/(norm*norm); 55 56 /* Evaluate function at F(x + dx) */ 57 ierr = VecWAXPY(&h,dx,U,w); CHKERRQ(ierr); 58 ierr = SNESComputeFunction(snes,w,y); CHKERRQ(ierr); 59 ierr = VecAXPY(&mone,F,y); CHKERRQ(ierr); 60 h = -1.0/h; 61 ierr = VecScale(&h,y); CHKERRQ(ierr); 62 if (ctx->sp) { ierr = PCNullSpaceRemove(ctx->sp,y); CHKERRQ(ierr);} 63 return 0; 64 } 65 /*@C 66 SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 67 context for use with a SNES solver. This matrix can be used as 68 the Jacobian argument for the routine SNESSetJacobian(). 69 70 Input Parameters: 71 . x - vector where SNES solution is to be stored. 72 73 Output Parameters: 74 . J - the matrix-free matrix 75 76 Notes: 77 The matrix-free matrix context merely contains the function pointers 78 and work space for performing finite difference approximations of 79 matrix operations such as matrix-vector products. 80 81 The user should call MatDestroy() when finished with the matrix-free 82 matrix context. 83 84 .keywords: SNES, default, matrix-free, create, matrix 85 86 .seealso: MatDestroy() 87 @*/ 88 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 89 { 90 MPI_Comm comm; 91 MFCtx_Private *mfctx; 92 int n,ierr; 93 94 mfctx = (MFCtx_Private *) PetscMalloc(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 95 PLogObjectMemory(snes,sizeof(MFCtx_Private)); 96 mfctx->sp = 0; 97 mfctx->snes = snes; 98 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 99 ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr); 100 ierr = VecGetSize(x,&n); CHKERRQ(ierr); 101 ierr = MatShellCreate(comm,n,n,(void*)mfctx,J); CHKERRQ(ierr); 102 ierr = MatShellSetMult(*J,SNESMatrixFreeMult_Private); CHKERRQ(ierr); 103 ierr = MatShellSetDestroy(*J,SNESMatrixFreeDestroy_Private); CHKERRQ(ierr); 104 PLogObjectParent(*J,mfctx->w); 105 PLogObjectParent(snes,*J); 106 return 0; 107 } 108 109 /*@ 110 SNESDefaultMatrixFreeMatAddNullSpace - Provide a null space that 111 an operator is suppose to have. Since round off will create a 112 small component in the null space, if you know the null space 113 you may have it automatically removed. 114 115 Input Parameters: 116 . J - the matrix free matrix 117 . has_cnst - PETSC_TRUE or PETSC_FALSE indicating if null space has constants 118 . n - number of vectors (excluding constant vector) in null space 119 . vecs - the vectors that span the null space (excluding the constant vector) 120 . these vectors must be orthonormal 121 122 .keywords: SNES, matrix-free, null space 123 @*/ 124 int SNESDefaultMatrixFreeMatAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs) 125 { 126 int ierr; 127 MFCtx_Private *ctx; 128 MPI_Comm comm; 129 130 PetscObjectGetComm((PetscObject)J,&comm); 131 132 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 133 /* no context indicates that it is not the "matrix free" matrix type */ 134 if (!ctx) return 0; 135 ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr); 136 return 0; 137 } 138 139