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