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