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