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