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