1 2 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 3 /* matimpl.h is needed only for logging of matrix operation */ 4 #include <petsc-private/matimpl.h> 5 6 extern PetscErrorCode SNESDiffParameterCreate_More(SNES,Vec,void**); 7 extern PetscErrorCode SNESDiffParameterCompute_More(SNES,void*,Vec,Vec,PetscReal*,PetscReal*); 8 extern PetscErrorCode SNESDiffParameterDestroy_More(void*); 9 10 typedef struct { /* default context for matrix-free SNES */ 11 SNES snes; /* SNES context */ 12 Vec w; /* work vector */ 13 MatNullSpace sp; /* null space context */ 14 PetscReal error_rel; /* square root of relative error in computing function */ 15 PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */ 16 PetscBool jorge; /* flag indicating use of Jorge's method for determining 17 the differencing parameter */ 18 PetscReal h; /* differencing parameter */ 19 PetscBool need_h; /* flag indicating whether we must compute h */ 20 PetscBool need_err; /* flag indicating whether we must currently compute error_rel */ 21 PetscBool compute_err; /* flag indicating whether we must ever compute error_rel */ 22 PetscInt compute_err_iter; /* last iter where we've computer error_rel */ 23 PetscInt compute_err_freq; /* frequency of computing error_rel */ 24 void *data; /* implementation-specific data */ 25 } MFCtx_Private; 26 27 #undef __FUNCT__ 28 #define __FUNCT__ "SNESMatrixFreeDestroy2_Private" 29 PetscErrorCode SNESMatrixFreeDestroy2_Private(Mat mat) 30 { 31 PetscErrorCode ierr; 32 MFCtx_Private *ctx; 33 34 PetscFunctionBegin; 35 ierr = MatShellGetContext(mat,(void**)&ctx);CHKERRQ(ierr); 36 ierr = VecDestroy(&ctx->w);CHKERRQ(ierr); 37 ierr = MatNullSpaceDestroy(&ctx->sp);CHKERRQ(ierr); 38 if (ctx->jorge || ctx->compute_err) {ierr = SNESDiffParameterDestroy_More(ctx->data);CHKERRQ(ierr);} 39 ierr = PetscFree(ctx);CHKERRQ(ierr); 40 PetscFunctionReturn(0); 41 } 42 43 #undef __FUNCT__ 44 #define __FUNCT__ "SNESMatrixFreeView2_Private" 45 /* 46 SNESMatrixFreeView2_Private - Views matrix-free parameters. 47 */ 48 PetscErrorCode SNESMatrixFreeView2_Private(Mat J,PetscViewer viewer) 49 { 50 PetscErrorCode ierr; 51 MFCtx_Private *ctx; 52 PetscBool iascii; 53 54 PetscFunctionBegin; 55 ierr = MatShellGetContext(J,(void**)&ctx);CHKERRQ(ierr); 56 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 57 if (iascii) { 58 ierr = PetscViewerASCIIPrintf(viewer," SNES matrix-free approximation:\n");CHKERRQ(ierr); 59 if (ctx->jorge) { 60 ierr = PetscViewerASCIIPrintf(viewer," using Jorge's method of determining differencing parameter\n");CHKERRQ(ierr); 61 } 62 ierr = PetscViewerASCIIPrintf(viewer," err=%G (relative error in function evaluation)\n",ctx->error_rel);CHKERRQ(ierr); 63 ierr = PetscViewerASCIIPrintf(viewer," umin=%G (minimum iterate parameter)\n",ctx->umin);CHKERRQ(ierr); 64 if (ctx->compute_err) { 65 ierr = PetscViewerASCIIPrintf(viewer," freq_err=%D (frequency for computing err)\n",ctx->compute_err_freq);CHKERRQ(ierr); 66 } 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; 87 Vec w,U,F; 88 PetscErrorCode ierr,(*eval_fct)(SNES,Vec,Vec); 89 MPI_Comm comm; 90 PetscInt iter; 91 92 PetscFunctionBegin; 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(MATMFFD_Mult,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 = SNESDiffParameterCompute_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 = SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);CHKERRQ(ierr); 123 124 ctx->error_rel = sqrt(noise); 125 126 ierr = PetscInfo3(snes,"Using Jorge's noise: noise=%G, sqrt(noise)=%G, h_more=%G\n",noise,ctx->error_rel,h);CHKERRQ(ierr); 127 128 ctx->compute_err_iter = iter; 129 ctx->need_err = PETSC_FALSE; 130 } 131 132 ierr = VecDotBegin(U,a,&dot);CHKERRQ(ierr); 133 ierr = VecNormBegin(a,NORM_1,&sum);CHKERRQ(ierr); 134 ierr = VecNormBegin(a,NORM_2,&norm);CHKERRQ(ierr); 135 ierr = VecDotEnd(U,a,&dot);CHKERRQ(ierr); 136 ierr = VecNormEnd(a,NORM_1,&sum);CHKERRQ(ierr); 137 ierr = VecNormEnd(a,NORM_2,&norm);CHKERRQ(ierr); 138 139 140 /* Safeguard for step sizes too small */ 141 if (sum == 0.0) { 142 dot = 1.0; 143 norm = 1.0; 144 } else if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum; 145 else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum; 146 h = PetscRealPart(ctx->error_rel*dot/(norm*norm)); 147 } 148 } else h = ctx->h; 149 150 if (!ctx->jorge || !ctx->need_h) {ierr = PetscInfo1(snes,"h = %G\n",h);CHKERRQ(ierr);} 151 152 /* Evaluate function at F(u + ha) */ 153 hs = h; 154 ierr = VecWAXPY(w,hs,a,U);CHKERRQ(ierr); 155 ierr = eval_fct(snes,w,y);CHKERRQ(ierr); 156 ierr = VecAXPY(y,-1.0,F);CHKERRQ(ierr); 157 ierr = VecScale(y,1.0/hs);CHKERRQ(ierr); 158 if (ctx->sp) {ierr = MatNullSpaceRemove(ctx->sp,y,PETSC_NULL);CHKERRQ(ierr);} 159 160 ierr = PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);CHKERRQ(ierr); 161 PetscFunctionReturn(0); 162 } 163 164 #undef __FUNCT__ 165 #define __FUNCT__ "SNESDefaultMatrixFreeCreate2" 166 /*@C 167 SNESMatrixFreeCreate2 - Creates a matrix-free matrix 168 context for use with a SNES solver. This matrix can be used as 169 the Jacobian argument for the routine SNESSetJacobian(). 170 171 Input Parameters: 172 . snes - the SNES context 173 . x - vector where SNES solution is to be stored. 174 175 Output Parameter: 176 . J - the matrix-free matrix 177 178 Level: advanced 179 180 Notes: 181 The matrix-free matrix context merely contains the function pointers 182 and work space for performing finite difference approximations of 183 Jacobian-vector products, J(u)*a, via 184 185 $ J(u)*a = [J(u+h*a) - J(u)]/h, 186 $ where by default 187 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 188 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 189 $ where 190 $ error_rel = square root of relative error in 191 $ function evaluation 192 $ umin = minimum iterate parameter 193 $ Alternatively, the differencing parameter, h, can be set using 194 $ Jorge's nifty new strategy if one specifies the option 195 $ -snes_mf_jorge 196 197 The user can set these parameters via MatMFFDSetFunctionError(). 198 See the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>. 199 200 The user should call MatDestroy() when finished with the matrix-free 201 matrix context. 202 203 Options Database Keys: 204 $ -snes_mf_err <error_rel> 205 $ -snes_mf_unim <umin> 206 $ -snes_mf_compute_err 207 $ -snes_mf_freq_err <freq> 208 $ -snes_mf_jorge 209 210 .keywords: SNES, default, matrix-free, create, matrix 211 212 .seealso: MatDestroy(), MatMFFDSetFunctionError() 213 @*/ 214 PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J) 215 { 216 MPI_Comm comm; 217 MFCtx_Private *mfctx; 218 PetscErrorCode ierr; 219 PetscInt n,nloc; 220 PetscBool flg; 221 char p[64]; 222 223 PetscFunctionBegin; 224 ierr = PetscNewLog(snes,MFCtx_Private,&mfctx);CHKERRQ(ierr); 225 mfctx->sp = 0; 226 mfctx->snes = snes; 227 mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON; 228 mfctx->umin = 1.e-6; 229 mfctx->h = 0.0; 230 mfctx->need_h = PETSC_TRUE; 231 mfctx->need_err = PETSC_FALSE; 232 mfctx->compute_err = PETSC_FALSE; 233 mfctx->compute_err_freq = 0; 234 mfctx->compute_err_iter = -1; 235 mfctx->compute_err = PETSC_FALSE; 236 mfctx->jorge = PETSC_FALSE; 237 238 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_mf_err",&mfctx->error_rel,PETSC_NULL);CHKERRQ(ierr); 239 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_mf_umin",&mfctx->umin,PETSC_NULL);CHKERRQ(ierr); 240 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_mf_jorge",&mfctx->jorge,PETSC_NULL);CHKERRQ(ierr); 241 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_mf_compute_err",&mfctx->compute_err,PETSC_NULL);CHKERRQ(ierr); 242 ierr = PetscOptionsGetInt(((PetscObject)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 = PETSC_TRUE; 246 } 247 if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE; 248 if (mfctx->jorge || mfctx->compute_err) { 249 ierr = SNESDiffParameterCreate_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 (((PetscObject)snes)->prefix) PetscStrcat(p,((PetscObject)snes)->prefix); 255 if (flg) { 256 ierr = PetscPrintf(((PetscObject)snes)->comm," Matrix-free Options (via SNES):\n");CHKERRQ(ierr); 257 ierr = PetscPrintf(((PetscObject)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(((PetscObject)snes)->comm," %ssnes_mf_umin <umin>: see users manual (default %G)\n",p,mfctx->umin);CHKERRQ(ierr); 259 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_jorge: use Jorge More's method\n",p);CHKERRQ(ierr); 260 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);CHKERRQ(ierr); 261 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);CHKERRQ(ierr); 262 ierr = PetscPrintf(((PetscObject)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,J);CHKERRQ(ierr); 269 ierr = MatSetSizes(*J,nloc,n,n,n);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