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