1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.36 1996/09/28 16:24:48 curfman Exp curfman $"; 4 #endif 5 6 #include "draw.h" /*I "draw.h" I*/ 7 #include "src/snes/snesimpl.h" /*I "snes.h" I*/ 8 #include "pinclude/pviewer.h" 9 10 typedef struct { /* default context for matrix-free SNES */ 11 SNES snes; /* SNES context */ 12 Vec w; /* work vector */ 13 PCNullSpace sp; /* null space context */ 14 double error_rel; /* square root of relative error in computing function */ 15 double umin; /* minimum allowable u value */ 16 } MFCtx_Private; 17 18 int SNESMatrixFreeDestroy_Private(PetscObject obj) 19 { 20 int ierr; 21 Mat mat = (Mat) obj; 22 MFCtx_Private *ctx; 23 24 ierr = MatShellGetContext(mat,(void **)&ctx); 25 ierr = VecDestroy(ctx->w); CHKERRQ(ierr); 26 if (ctx->sp) {ierr = PCNullSpaceDestroy(ctx->sp); CHKERRQ(ierr);} 27 PetscFree(ctx); 28 return 0; 29 } 30 31 /* 32 SNESMatrixFreeView_Private - Views matrix-free parameters. 33 */ 34 int SNESMatrixFreeView_Private(Mat J,Viewer viewer) 35 { 36 int ierr; 37 MFCtx_Private *ctx; 38 MPI_Comm comm; 39 FILE *fd; 40 ViewerType vtype; 41 42 PetscObjectGetComm((PetscObject)J,&comm); 43 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 44 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 45 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 46 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 47 PetscFPrintf(comm,fd," SNES matrix-free approximation:\n"); 48 PetscFPrintf(comm,fd," err=%g (relative error in function evaluation)\n",ctx->error_rel); 49 PetscFPrintf(comm,fd," umin=%g (minimum iterate parameter)\n",ctx->umin); 50 } 51 return 0; 52 } 53 54 /* 55 SNESMatrixFreeMult_Private - Default matrix-free form for Jacobian-vector 56 product, y = F'(u)*a: 57 y = ( F(u + ha) - F(u) ) /h, 58 where F = nonlinear function, as set by SNESSetFunction() 59 u = current iterate 60 h = difference interval 61 */ 62 int SNESMatrixFreeMult_Private(Mat mat,Vec a,Vec y) 63 { 64 MFCtx_Private *ctx; 65 SNES snes; 66 double norm, sum, umin; 67 Scalar h, dot, mone = -1.0; 68 Vec w,U,F; 69 int ierr, (*eval_fct)(SNES,Vec,Vec); 70 71 MatShellGetContext(mat,(void **)&ctx); 72 snes = ctx->snes; 73 w = ctx->w; 74 umin = ctx->umin; 75 76 /* We log matrix-free matrix-vector products separately, so that we can 77 separate the performance monitoring from the cases that use conventional 78 storage. We may eventually modify event logging to associate events 79 with particular objects, hence alleviating the more general problem. */ 80 PLogEventBegin(MAT_MatrixFreeMult,a,y,0,0); 81 82 ierr = SNESGetSolution(snes,&U); CHKERRQ(ierr); 83 if (snes->method_class == SNES_NONLINEAR_EQUATIONS) { 84 eval_fct = SNESComputeFunction; 85 ierr = SNESGetFunction(snes,&F); CHKERRQ(ierr); 86 } 87 else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) { 88 eval_fct = SNESComputeGradient; 89 ierr = SNESGetGradient(snes,&F); CHKERRQ(ierr); 90 } 91 else SETERRQ(1,"SNESMatrixFreeMult_Private: Invalid method class"); 92 93 /* Determine a "good" step size, h */ 94 ierr = VecDot(U,a,&dot); CHKERRQ(ierr); 95 ierr = VecNorm(a,NORM_1,&sum); CHKERRQ(ierr); 96 ierr = VecNorm(a,NORM_2,&norm); CHKERRQ(ierr); 97 98 /* Safeguard for step sizes too small */ 99 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 100 #if defined(PETSC_COMPLEX) 101 else if (abs(dot) < umin*sum && real(dot) >= 0.0) dot = umin*sum; 102 else if (abs(dot) < 0.0 && real(dot) > -umin*sum) dot = -umin*sum; 103 #else 104 else if (dot < umin*sum && dot >= 0.0) dot = umin*sum; 105 else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum; 106 #endif 107 h = ctx->error_rel*dot/(norm*norm); 108 109 /* Evaluate function at F(u + ha) */ 110 ierr = VecWAXPY(&h,a,U,w); CHKERRQ(ierr); 111 ierr = eval_fct(snes,w,y); CHKERRQ(ierr); 112 ierr = VecAXPY(&mone,F,y); CHKERRQ(ierr); 113 h = 1.0/h; 114 ierr = VecScale(&h,y); CHKERRQ(ierr); 115 if (ctx->sp) {ierr = PCNullSpaceRemove(ctx->sp,y); CHKERRQ(ierr);} 116 117 PLogEventEnd(MAT_MatrixFreeMult,a,y,0,0); 118 return 0; 119 } 120 121 /*@C 122 SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 123 context for use with a SNES solver. This matrix can be used as 124 the Jacobian argument for the routine SNESSetJacobian(). 125 126 Input Parameters: 127 . snes - the SNES context 128 . x - vector where SNES solution is to be stored. 129 130 Output Parameter: 131 . J - the matrix-free matrix 132 133 Notes: 134 The matrix-free matrix context merely contains the function pointers 135 and work space for performing finite difference approximations of 136 matrix operations such as matrix-vector products. 137 138 The user should call MatDestroy() when finished with the matrix-free 139 matrix context. 140 141 .keywords: SNES, default, matrix-free, create, matrix 142 143 .seealso: MatDestroy() 144 @*/ 145 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 146 { 147 MPI_Comm comm; 148 MFCtx_Private *mfctx; 149 int n, nloc, ierr, flg; 150 151 mfctx = (MFCtx_Private *) PetscMalloc(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 152 PLogObjectMemory(snes,sizeof(MFCtx_Private)); 153 mfctx->sp = 0; 154 mfctx->snes = snes; 155 mfctx->error_rel = 1.e-8; /* assumes double precision */ 156 mfctx->umin = 1.e-8; 157 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_err",&mfctx->error_rel,&flg); CHKERRQ(ierr); 158 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_umin",&mfctx->umin,&flg); CHKERRQ(ierr); 159 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 160 ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr); 161 ierr = VecGetSize(x,&n); CHKERRQ(ierr); 162 ierr = VecGetLocalSize(x,&nloc); CHKERRQ(ierr); 163 ierr = MatCreateShell(comm,nloc,n,n,n,(void*)mfctx,J); CHKERRQ(ierr); 164 ierr = MatShellSetOperation(*J,MATOP_MULT,(void*)SNESMatrixFreeMult_Private); CHKERRQ(ierr); 165 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void *)SNESMatrixFreeDestroy_Private); CHKERRQ(ierr); 166 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void *)SNESMatrixFreeView_Private); CHKERRQ(ierr); 167 PLogObjectParent(*J,mfctx->w); 168 PLogObjectParent(snes,*J); 169 return 0; 170 } 171 172 /*@ 173 SNESSetMatrixFreeParameters - Sets the parameters for the approximation of 174 matrix-vector products using finite differences. 175 176 Input Parameters: 177 . snes - the SNES context 178 . error_rel - relative error 179 . umin - minimum allowable u-value 180 181 .keywords: SNES, matrix-free, parameters 182 @*/ 183 int SNESSetMatrixFreeParameters(SNES snes,double error,double umin) 184 { 185 MFCtx_Private *ctx; 186 int ierr; 187 Mat mat; 188 189 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 190 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 191 if (ctx) { 192 if (error != PETSC_DEFAULT) ctx->error_rel = error; 193 if (umin != PETSC_DEFAULT) ctx->umin = umin; 194 } 195 return 0; 196 } 197 198 /*@ 199 SNESDefaultMatrixFreeMatAddNullSpace - Provides a null space that 200 an operator is supposed to have. Since roundoff will create a 201 small component in the null space, if you know the null space 202 you may have it automatically removed. 203 204 Input Parameters: 205 . J - the matrix-free matrix context 206 . has_cnst - PETSC_TRUE or PETSC_FALSE, indicating if null space has constants 207 . n - number of vectors (excluding constant vector) in null space 208 . vecs - the vectors that span the null space (excluding the constant vector); 209 . these vectors must be orthonormal 210 211 .keywords: SNES, matrix-free, null space 212 @*/ 213 int SNESDefaultMatrixFreeMatAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs) 214 { 215 int ierr; 216 MFCtx_Private *ctx; 217 MPI_Comm comm; 218 219 PetscObjectGetComm((PetscObject)J,&comm); 220 221 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 222 /* no context indicates that it is not the "matrix free" matrix type */ 223 if (!ctx) return 0; 224 ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr); 225 return 0; 226 } 227 228 229 230 231