1 2 #include "src/snes/snesimpl.h" /*I "petscsnes.h" I*/ 3 4 EXTERN PetscErrorCode DiffParameterCreate_More(SNES,Vec,void**); 5 EXTERN PetscErrorCode DiffParameterCompute_More(SNES,void*,Vec,Vec,PetscReal*,PetscReal*); 6 EXTERN PetscErrorCode 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 PetscTruth jorge; /* flag indicating use of Jorge's method for determining 15 the differencing parameter */ 16 PetscReal h; /* differencing parameter */ 17 PetscTruth need_h; /* flag indicating whether we must compute h */ 18 PetscTruth need_err; /* flag indicating whether we must currently compute error_rel */ 19 PetscTruth compute_err; /* flag indicating whether we must ever compute error_rel */ 20 PetscInt compute_err_iter; /* last iter where we've computer error_rel */ 21 PetscInt 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 PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat) 28 { 29 PetscErrorCode 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 PetscErrorCode SNESMatrixFreeView2_Private(Mat J,PetscViewer viewer) 47 { 48 PetscErrorCode 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(PETSC_ERR_SUP,"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 PetscErrorCode 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 PetscErrorCode ierr,(*eval_fct)(SNES,Vec,Vec); 89 MPI_Comm comm; 90 PetscInt iter; 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(MATSNESMF_Mult,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 eval_fct = SNESComputeFunction; 108 ierr = SNESGetFunction(snes,&F,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 109 110 /* Determine a "good" step size, h */ 111 if (ctx->need_h) { 112 113 /* Use Jorge's method to compute h */ 114 if (ctx->jorge) { 115 ierr = DiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);CHKERRQ(ierr); 116 117 /* Use the Brown/Saad method to compute h */ 118 } else { 119 /* Compute error if desired */ 120 ierr = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr); 121 if ((ctx->need_err) || ((ctx->compute_err_freq) && (ctx->compute_err_iter != iter) && (!((iter-1)%ctx->compute_err_freq)))) { 122 /* Use Jorge's method to compute noise */ 123 ierr = DiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);CHKERRQ(ierr); 124 ctx->error_rel = sqrt(noise); 125 ierr = PetscLogInfo((snes,"SNESMatrixFreeMult2_Private: Using Jorge's noise: noise=%g, sqrt(noise)=%g, h_more=%g\n",noise,ctx->error_rel,h));CHKERRQ(ierr); 126 ctx->compute_err_iter = iter; 127 ctx->need_err = PETSC_FALSE; 128 } 129 130 ierr = VecDotBegin(U,a,&dot);CHKERRQ(ierr); 131 ierr = VecNormBegin(a,NORM_1,&sum);CHKERRQ(ierr); 132 ierr = VecNormBegin(a,NORM_2,&norm);CHKERRQ(ierr); 133 ierr = VecDotEnd(U,a,&dot);CHKERRQ(ierr); 134 ierr = VecNormEnd(a,NORM_1,&sum);CHKERRQ(ierr); 135 ierr = VecNormEnd(a,NORM_2,&norm);CHKERRQ(ierr); 136 137 138 /* Safeguard for step sizes too small */ 139 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 140 #if defined(PETSC_USE_COMPLEX) 141 else if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum; 142 else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum; 143 #else 144 else if (dot < umin*sum && dot >= 0.0) dot = umin*sum; 145 else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum; 146 #endif 147 h = PetscRealPart(ctx->error_rel*dot/(norm*norm)); 148 } 149 } else { 150 h = ctx->h; 151 } 152 if (!ctx->jorge || !ctx->need_h) {ierr = PetscLogInfo((snes,"SNESMatrixFreeMult2_Private: h = %g\n",h));CHKERRQ(ierr);} 153 154 /* Evaluate function at F(u + ha) */ 155 hs = h; 156 ierr = VecWAXPY(&hs,a,U,w);CHKERRQ(ierr); 157 ierr = eval_fct(snes,w,y);CHKERRQ(ierr); 158 ierr = VecAXPY(&mone,F,y);CHKERRQ(ierr); 159 hs = 1.0/hs; 160 ierr = VecScale(&hs,y);CHKERRQ(ierr); 161 if (ctx->sp) {ierr = MatNullSpaceRemove(ctx->sp,y,PETSC_NULL);CHKERRQ(ierr);} 162 163 ierr = PetscLogEventEnd(MATSNESMF_Mult,a,y,0,0);CHKERRQ(ierr); 164 PetscFunctionReturn(0); 165 } 166 167 #undef __FUNCT__ 168 #define __FUNCT__ "SNESMatrixFreeMatCreate2" 169 /*@C 170 SNESMatrixFreeMatCreate2 - Creates a matrix-free matrix 171 context for use with a SNES solver. This matrix can be used as 172 the Jacobian argument for the routine SNESSetJacobian(). 173 174 Input Parameters: 175 . snes - the SNES context 176 . x - vector where SNES solution is to be stored. 177 178 Output Parameter: 179 . J - the matrix-free matrix 180 181 Level: advanced 182 183 Notes: 184 The matrix-free matrix context merely contains the function pointers 185 and work space for performing finite difference approximations of 186 Jacobian-vector products, J(u)*a, via 187 188 $ J(u)*a = [J(u+h*a) - J(u)]/h, 189 $ where by default 190 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 191 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 192 $ where 193 $ error_rel = square root of relative error in 194 $ function evaluation 195 $ umin = minimum iterate parameter 196 $ Alternatively, the differencing parameter, h, can be set using 197 $ Jorge's nifty new strategy if one specifies the option 198 $ -snes_mf_jorge 199 200 The user can set these parameters via MatSNESMFSetFunctionError(). 201 See the nonlinear solvers chapter of the users manual for details. 202 203 The user should call MatDestroy() when finished with the matrix-free 204 matrix context. 205 206 Options Database Keys: 207 $ -snes_mf_err <error_rel> 208 $ -snes_mf_unim <umin> 209 $ -snes_mf_compute_err 210 $ -snes_mf_freq_err <freq> 211 $ -snes_mf_jorge 212 213 .keywords: SNES, default, matrix-free, create, matrix 214 215 .seealso: MatDestroy(), MatSNESMFSetFunctionError() 216 @*/ 217 PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J) 218 { 219 MPI_Comm comm; 220 MFCtx_Private *mfctx; 221 PetscErrorCode ierr; 222 PetscInt n,nloc; 223 PetscTruth flg; 224 char p[64]; 225 226 PetscFunctionBegin; 227 ierr = PetscNew(MFCtx_Private,&mfctx);CHKERRQ(ierr); 228 ierr = PetscLogObjectMemory(snes,sizeof(MFCtx_Private));CHKERRQ(ierr); 229 mfctx->sp = 0; 230 mfctx->snes = snes; 231 mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON; 232 mfctx->umin = 1.e-6; 233 mfctx->h = 0.0; 234 mfctx->need_h = PETSC_TRUE; 235 mfctx->need_err = PETSC_FALSE; 236 mfctx->compute_err = PETSC_FALSE; 237 mfctx->compute_err_freq = 0; 238 mfctx->compute_err_iter = -1; 239 ierr = PetscOptionsGetReal(snes->prefix,"-snes_mf_err",&mfctx->error_rel,PETSC_NULL);CHKERRQ(ierr); 240 ierr = PetscOptionsGetReal(snes->prefix,"-snes_mf_umin",&mfctx->umin,PETSC_NULL);CHKERRQ(ierr); 241 ierr = PetscOptionsHasName(snes->prefix,"-snes_mf_jorge",&mfctx->jorge);CHKERRQ(ierr); 242 ierr = PetscOptionsHasName(snes->prefix,"-snes_mf_compute_err",&mfctx->compute_err);CHKERRQ(ierr); 243 ierr = PetscOptionsGetInt(snes->prefix,"-snes_mf_freq_err",&mfctx->compute_err_freq,&flg);CHKERRQ(ierr); 244 if (flg) { 245 if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0; 246 mfctx->compute_err = PETSC_TRUE; 247 } 248 if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE; 249 if (mfctx->jorge || mfctx->compute_err) { 250 ierr = DiffParameterCreate_More(snes,x,&mfctx->data);CHKERRQ(ierr); 251 } else mfctx->data = 0; 252 253 ierr = PetscOptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 254 ierr = PetscStrcpy(p,"-");CHKERRQ(ierr); 255 if (snes->prefix) PetscStrcat(p,snes->prefix); 256 if (flg) { 257 ierr = PetscPrintf(snes->comm," Matrix-free Options (via SNES):\n");CHKERRQ(ierr); 258 ierr = PetscPrintf(snes->comm," %ssnes_mf_err <err>: set sqrt of relative error in function (default %g)\n",p,mfctx->error_rel);CHKERRQ(ierr); 259 ierr = PetscPrintf(snes->comm," %ssnes_mf_umin <umin>: see users manual (default %g)\n",p,mfctx->umin);CHKERRQ(ierr); 260 ierr = PetscPrintf(snes->comm," %ssnes_mf_jorge: use Jorge More's method\n",p);CHKERRQ(ierr); 261 ierr = PetscPrintf(snes->comm," %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);CHKERRQ(ierr); 262 ierr = PetscPrintf(snes->comm," %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);CHKERRQ(ierr); 263 ierr = PetscPrintf(snes->comm," %ssnes_mf_noise_file <file>: set file for printing noise info\n",p);CHKERRQ(ierr); 264 } 265 ierr = VecDuplicate(x,&mfctx->w);CHKERRQ(ierr); 266 ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr); 267 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 268 ierr = VecGetLocalSize(x,&nloc);CHKERRQ(ierr); 269 ierr = MatCreate(comm,nloc,n,n,n,J);CHKERRQ(ierr); 270 ierr = MatSetType(*J,MATSHELL);CHKERRQ(ierr); 271 ierr = MatShellSetContext(*J,mfctx);CHKERRQ(ierr); 272 ierr = MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))SNESMatrixFreeMult2_Private);CHKERRQ(ierr); 273 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))SNESMatrixFreeDestroy2_Private);CHKERRQ(ierr); 274 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))SNESMatrixFreeView2_Private);CHKERRQ(ierr); 275 276 ierr = PetscLogObjectParent(*J,mfctx->w);CHKERRQ(ierr); 277 ierr = PetscLogObjectParent(snes,*J);CHKERRQ(ierr); 278 PetscFunctionReturn(0); 279 } 280 281 #undef __FUNCT__ 282 #define __FUNCT__ "SNESDefaultMatrixFreeSetParameters2" 283 /*@C 284 SNESDefaultMatrixFreeSetParameters2 - Sets the parameters for the approximation of 285 matrix-vector products using finite differences. 286 287 $ J(u)*a = [J(u+h*a) - J(u)]/h where 288 289 either the user sets h directly here, or this parameter is computed via 290 291 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 292 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 293 $ 294 295 Input Parameters: 296 + mat - the matrix 297 . error_rel - relative error (should be set to the square root of 298 the relative error in the function evaluations) 299 . umin - minimum allowable u-value 300 - h - differencing parameter 301 302 Level: advanced 303 304 Notes: 305 If the user sets the parameter h directly, then this value will be used 306 instead of the default computation indicated above. 307 308 .keywords: SNES, matrix-free, parameters 309 310 .seealso: MatCreateSNESMF() 311 @*/ 312 PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat mat,PetscReal error,PetscReal umin,PetscReal h) 313 { 314 MFCtx_Private *ctx; 315 PetscErrorCode ierr; 316 317 PetscFunctionBegin; 318 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 319 if (ctx) { 320 if (error != PETSC_DEFAULT) ctx->error_rel = error; 321 if (umin != PETSC_DEFAULT) ctx->umin = umin; 322 if (h != PETSC_DEFAULT) { 323 ctx->h = h; 324 ctx->need_h = PETSC_FALSE; 325 } 326 } 327 PetscFunctionReturn(0); 328 } 329 330 PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes) 331 { 332 MFCtx_Private *ctx; 333 PetscErrorCode ierr; 334 Mat mat; 335 336 PetscFunctionBegin; 337 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 338 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 339 if (ctx) ctx->need_h = PETSC_TRUE; 340 PetscFunctionReturn(0); 341 } 342 343