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