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