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