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