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