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