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