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