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