1 #ifndef lint 2 static char vcid[] = "$Id: snesmfj2.c,v 1.3 1997/07/06 16:37:32 curfman Exp curfman $"; 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 Scalar 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(PetscObject obj) 33 { 34 int ierr; 35 Mat mat = (Mat) obj; 36 MFCtx_Private *ctx; 37 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 return 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 PetscObjectGetComm((PetscObject)J,&comm); 60 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 61 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 62 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 63 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 64 PetscFPrintf(comm,fd," SNES matrix-free approximation:\n"); 65 if (ctx->jorge) 66 PetscFPrintf(comm,fd," using Jorge's method of determining differencing parameter\n"); 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 return 0; 73 } 74 75 extern int VecDot_Seq(Vec,Vec,Scalar *); 76 extern int VecNorm_Seq(Vec,NormType,double *); 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 norm, sum, umin, noise, ovalues[3],values[3]; 93 Scalar h, dot, mone = -1.0; 94 Vec w,U,F; 95 int ierr, iter, (*eval_fct)(SNES,Vec,Vec); 96 MPI_Comm comm; 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 PetscObjectGetComm((PetscObject)mat,&comm); 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); CHKERRQ(ierr); 114 } 115 else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) { 116 eval_fct = SNESComputeGradient; 117 ierr = SNESGetGradient(snes,&F); 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 /* 145 ierr = VecDot(U,a,&dot); CHKERRQ(ierr); 146 ierr = VecNorm(a,NORM_1,&sum); CHKERRQ(ierr); 147 ierr = VecNorm(a,NORM_2,&norm); CHKERRQ(ierr); 148 */ 149 150 /* 151 Call the Seq Vector routines and then do a single reduction 152 to reduce the number of communications required 153 */ 154 155 #if !defined(PETSC_COMPLEX) 156 PLogEventBegin(VEC_Dot,U,a,0,0); 157 ierr = VecDot_Seq(U,a,ovalues); CHKERRQ(ierr); 158 PLogEventEnd(VEC_Dot,U,a,0,0); 159 PLogEventBegin(VEC_Norm,a,0,0,0); 160 ierr = VecNorm_Seq(a,NORM_1,ovalues+1); CHKERRQ(ierr); 161 ierr = VecNorm_Seq(a,NORM_2,ovalues+2); CHKERRQ(ierr); 162 ovalues[2] = ovalues[2]*ovalues[2]; 163 MPI_Allreduce(ovalues,values,3,MPI_DOUBLE,MPI_SUM,comm ); 164 dot = values[0]; sum = values[1]; norm = sqrt(values[2]); 165 PLogEventEnd(VEC_Norm,a,0,0,0); 166 #else 167 { 168 Scalar cvalues[3],covalues[3]; 169 170 PLogEventBegin(VEC_Dot,U,a,0,0); 171 ierr = VecDot_Seq(U,a,covalues); CHKERRQ(ierr); 172 PLogEventEnd(VEC_Dot,U,a,0,0); 173 PLogEventBegin(VEC_Norm,a,0,0,0); 174 ierr = VecNorm_Seq(a,NORM_1,ovalues+1); CHKERRQ(ierr); 175 ierr = VecNorm_Seq(a,NORM_2,ovalues+2); CHKERRQ(ierr); 176 covalues[1] = ovalues[1]; 177 covalues[2] = ovalues[2]*ovalues[2]; 178 MPI_Allreduce(covalues,cvalues,6,MPI_DOUBLE,MPI_SUM,comm ); 179 dot = cvalues[0]; sum = PetscReal(cvalues[1]); norm = sqrt(PetscReal(cvalues[2])); 180 PLogEventEnd(VEC_Norm,a,0,0,0); 181 } 182 #endif 183 184 /* Safeguard for step sizes too small */ 185 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 186 #if defined(PETSC_COMPLEX) 187 else if (abs(dot) < umin*sum && real(dot) >= 0.0) dot = umin*sum; 188 else if (abs(dot) < 0.0 && real(dot) > -umin*sum) dot = -umin*sum; 189 #else 190 else if (dot < umin*sum && dot >= 0.0) dot = umin*sum; 191 else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum; 192 #endif 193 h = ctx->error_rel*dot/(norm*norm); 194 } 195 } else { 196 h = ctx->h; 197 } 198 if (!ctx->jorge || !ctx->need_h) PLogInfo(snes,"SNESMatrixFreeMult2_Private: h = %g\n",h); 199 200 /* Evaluate function at F(u + ha) */ 201 ierr = VecWAXPY(&h,a,U,w); CHKERRQ(ierr); 202 ierr = eval_fct(snes,w,y); CHKERRQ(ierr); 203 ierr = VecAXPY(&mone,F,y); CHKERRQ(ierr); 204 h = 1.0/h; 205 ierr = VecScale(&h,y); CHKERRQ(ierr); 206 if (ctx->sp) {ierr = PCNullSpaceRemove(ctx->sp,y); CHKERRQ(ierr);} 207 208 PLogEventEnd(MAT_MatrixFreeMult,a,y,0,0); 209 return 0; 210 } 211 212 #undef __FUNC__ 213 #define __FUNC__ "SNESDefaultMatrixFreeMatCreate2" 214 /*@C 215 SNESMatrixFreeMatCreate2 - Creates a matrix-free matrix 216 context for use with a SNES solver. This matrix can be used as 217 the Jacobian argument for the routine SNESSetJacobian(). 218 219 Input Parameters: 220 . snes - the SNES context 221 . x - vector where SNES solution is to be stored. 222 223 Output Parameter: 224 . J - the matrix-free matrix 225 226 Notes: 227 The matrix-free matrix context merely contains the function pointers 228 and work space for performing finite difference approximations of 229 Jacobian-vector products, J(u)*a, via 230 231 $ J(u)*a = [J(u+h*a) - J(u)]/h, 232 $ where by default 233 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 234 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 235 $ where 236 $ error_rel = square root of relative error in 237 $ function evaluation 238 $ umin = minimum iterate parameter 239 $ Alternatively, the differencing parameter, h, can be set using 240 $ Jorge's nifty new strategy if one specifies the option 241 $ -snes_mf_jorge 242 243 The user can set these parameters via SNESSetMatrixFreeParameters(). 244 See the nonlinear solvers chapter of the users manual for details. 245 246 The user should call MatDestroy() when finished with the matrix-free 247 matrix context. 248 249 Options Database Keys: 250 $ -snes_mf_err <error_rel> 251 $ -snes_mf_unim <umin> 252 $ -snes_mf_compute_err 253 $ -snes_mf_freq_err <freq> 254 $ -snes_mf_jorge 255 256 .keywords: SNES, default, matrix-free, create, matrix 257 258 .seealso: MatDestroy(), SNESSetMatrixFreeParameters() 259 @*/ 260 int SNESMatrixFreeMatCreate2(SNES snes,Vec x, Mat *J) 261 { 262 MPI_Comm comm; 263 MFCtx_Private *mfctx; 264 int n, nloc, ierr, flg; 265 char p[64]; 266 267 mfctx = (MFCtx_Private *) PetscMalloc(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 268 PLogObjectMemory(snes,sizeof(MFCtx_Private)); 269 mfctx->sp = 0; 270 mfctx->snes = snes; 271 mfctx->error_rel = 1.e-8; /* assumes double precision */ 272 mfctx->umin = 1.e-6; 273 mfctx->h = 0.0; 274 mfctx->need_h = 1; 275 mfctx->need_err = 0; 276 mfctx->compute_err = 0; 277 mfctx->compute_err_freq = 0; 278 mfctx->compute_err_iter = -1; 279 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_err",&mfctx->error_rel,&flg); CHKERRQ(ierr); 280 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_umin",&mfctx->umin,&flg); CHKERRQ(ierr); 281 ierr = OptionsHasName(snes->prefix,"-snes_mf_jorge",&mfctx->jorge); CHKERRQ(ierr); 282 ierr = OptionsHasName(snes->prefix,"-snes_mf_compute_err",&mfctx->compute_err); CHKERRQ(ierr); 283 ierr = OptionsGetInt(snes->prefix,"-snes_mf_freq_err",&mfctx->compute_err_freq,&flg); CHKERRQ(ierr); 284 if (flg) { 285 if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0; 286 mfctx->compute_err = 1; 287 } 288 if (mfctx->compute_err == 1) mfctx->need_err = 1; 289 if (mfctx->jorge || mfctx->compute_err) { 290 ierr = DiffParameterCreate_More(snes,x,&mfctx->data); CHKERRQ(ierr); 291 } else mfctx->data = 0; 292 293 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 294 PetscStrcpy(p,"-"); 295 if (snes->prefix) PetscStrcat(p,snes->prefix); 296 if (flg) { 297 PetscPrintf(snes->comm," Matrix-free Options (via SNES):\n"); 298 PetscPrintf(snes->comm," %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n",p,mfctx->error_rel); 299 PetscPrintf(snes->comm," %ssnes_mf_umin <umin>: see users manual (default %g)\n",p,mfctx->umin); 300 PetscPrintf(snes->comm," %ssnes_mf_jorge: use Jorge More's method\n",p); 301 PetscPrintf(snes->comm," %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p); 302 PetscPrintf(snes->comm," %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p); 303 PetscPrintf(snes->comm," %ssnes_mf_noise_file <file>: set file for printing noise info\n",p); 304 } 305 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 306 ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr); 307 ierr = VecGetSize(x,&n); CHKERRQ(ierr); 308 ierr = VecGetLocalSize(x,&nloc); CHKERRQ(ierr); 309 ierr = MatCreateShell(comm,nloc,n,n,n,mfctx,J); CHKERRQ(ierr); 310 ierr = MatShellSetOperation(*J,MATOP_MULT,(void*)SNESMatrixFreeMult2_Private);CHKERRQ(ierr); 311 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void *)SNESMatrixFreeDestroy2_Private);CHKERRQ(ierr); 312 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void *)SNESMatrixFreeView2_Private); CHKERRQ(ierr); 313 314 PLogObjectParent(*J,mfctx->w); 315 PLogObjectParent(snes,*J); 316 return 0; 317 } 318 319 #undef __FUNC__ 320 #define __FUNC__ "SNESSetMatrixFreeParameters2" 321 /*@ 322 SNESSetMatrixFreeParameters2 - Sets the parameters for the approximation of 323 matrix-vector products using finite differences. 324 325 $ J(u)*a = [J(u+h*a) - J(u)]/h where 326 327 either the user sets h directly here, or this parameter is computed via 328 329 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 330 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 331 $ 332 333 Input Parameters: 334 . snes - the SNES context 335 . error_rel - relative error (should be set to the square root of 336 the relative error in the function evaluations) 337 . umin - minimum allowable u-value 338 . h - differencing parameter 339 340 Notes: 341 If the user sets the parameter h directly, then this value will be used 342 instead of the default computation indicated above. 343 344 .keywords: SNES, matrix-free, parameters 345 346 .seealso: SNESDefaultMatrixFreeMatCreate() 347 @*/ 348 int SNESSetMatrixFreeParameters2(SNES snes,double error,double umin,double h) 349 { 350 MFCtx_Private *ctx; 351 int ierr; 352 Mat mat; 353 354 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 355 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 356 if (ctx) { 357 if (error != PETSC_DEFAULT) ctx->error_rel = error; 358 if (umin != PETSC_DEFAULT) ctx->umin = umin; 359 if (h != PETSC_DEFAULT) { 360 ctx->h = h; 361 ctx->need_h = 0; 362 } 363 } 364 return 0; 365 } 366 367 int SNESUnSetMatrixFreeParameter(SNES snes) 368 { 369 MFCtx_Private *ctx; 370 int ierr; 371 Mat mat; 372 373 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 374 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 375 if (ctx) ctx->need_h = 1; 376 return 0; 377 } 378 379