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