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,(void**)&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,(void**)&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,(void**)&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 137 /* Safeguard for step sizes too small */ 138 if (sum == 0.0) { 139 dot = 1.0; 140 norm = 1.0; 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 h = PetscRealPart(ctx->error_rel*dot/(norm*norm)); 144 } 145 } else h = ctx->h; 146 147 if (!ctx->jorge || !ctx->need_h) {ierr = PetscInfo1(snes,"h = %g\n",(double)h);CHKERRQ(ierr);} 148 149 /* Evaluate function at F(u + ha) */ 150 hs = h; 151 ierr = VecWAXPY(w,hs,a,U);CHKERRQ(ierr); 152 ierr = eval_fct(snes,w,y);CHKERRQ(ierr); 153 ierr = VecAXPY(y,-1.0,F);CHKERRQ(ierr); 154 ierr = VecScale(y,1.0/hs);CHKERRQ(ierr); 155 if (mat->nullsp) {ierr = MatNullSpaceRemove(mat->nullsp,y);CHKERRQ(ierr);} 156 157 ierr = PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);CHKERRQ(ierr); 158 PetscFunctionReturn(0); 159 } 160 161 /*@C 162 SNESMatrixFreeCreate2 - Creates a matrix-free matrix 163 context for use with a SNES solver. This matrix can be used as 164 the Jacobian argument for the routine SNESSetJacobian(). 165 166 Input Parameters: 167 . snes - the SNES context 168 . x - vector where SNES solution is to be stored. 169 170 Output Parameter: 171 . J - the matrix-free matrix 172 173 Level: advanced 174 175 Notes: 176 The matrix-free matrix context merely contains the function pointers 177 and work space for performing finite difference approximations of 178 Jacobian-vector products, J(u)*a, via 179 180 $ J(u)*a = [J(u+h*a) - J(u)]/h, 181 $ where by default 182 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 183 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 184 $ where 185 $ error_rel = square root of relative error in 186 $ function evaluation 187 $ umin = minimum iterate parameter 188 $ Alternatively, the differencing parameter, h, can be set using 189 $ Jorge's nifty new strategy if one specifies the option 190 $ -snes_mf_jorge 191 192 The user can set these parameters via MatMFFDSetFunctionError(). 193 See Users-Manual: ch_snes for details 194 195 The user should call MatDestroy() when finished with the matrix-free 196 matrix context. 197 198 Options Database Keys: 199 $ -snes_mf_err <error_rel> 200 $ -snes_mf_unim <umin> 201 $ -snes_mf_compute_err 202 $ -snes_mf_freq_err <freq> 203 $ -snes_mf_jorge 204 205 .keywords: SNES, default, matrix-free, create, matrix 206 207 .seealso: MatDestroy(), MatMFFDSetFunctionError() 208 @*/ 209 PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J) 210 { 211 MPI_Comm comm; 212 MFCtx_Private *mfctx; 213 PetscErrorCode ierr; 214 PetscInt n,nloc; 215 PetscBool flg; 216 char p[64]; 217 218 PetscFunctionBegin; 219 ierr = PetscNewLog(snes,&mfctx);CHKERRQ(ierr); 220 mfctx->sp = 0; 221 mfctx->snes = snes; 222 mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON; 223 mfctx->umin = 1.e-6; 224 mfctx->h = 0.0; 225 mfctx->need_h = PETSC_TRUE; 226 mfctx->need_err = PETSC_FALSE; 227 mfctx->compute_err = PETSC_FALSE; 228 mfctx->compute_err_freq = 0; 229 mfctx->compute_err_iter = -1; 230 mfctx->compute_err = PETSC_FALSE; 231 mfctx->jorge = PETSC_FALSE; 232 233 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_err",&mfctx->error_rel,NULL);CHKERRQ(ierr); 234 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_umin",&mfctx->umin,NULL);CHKERRQ(ierr); 235 ierr = PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_jorge",&mfctx->jorge,NULL);CHKERRQ(ierr); 236 ierr = PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_compute_err",&mfctx->compute_err,NULL);CHKERRQ(ierr); 237 ierr = PetscOptionsGetInt(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_freq_err",&mfctx->compute_err_freq,&flg);CHKERRQ(ierr); 238 if (flg) { 239 if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0; 240 mfctx->compute_err = PETSC_TRUE; 241 } 242 if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE; 243 if (mfctx->jorge || mfctx->compute_err) { 244 ierr = SNESDiffParameterCreate_More(snes,x,&mfctx->data);CHKERRQ(ierr); 245 } else mfctx->data = 0; 246 247 ierr = PetscOptionsHasHelp(((PetscObject)snes)->options,&flg);CHKERRQ(ierr); 248 ierr = PetscStrncpy(p,"-",sizeof(p));CHKERRQ(ierr); 249 if (((PetscObject)snes)->prefix) {ierr = PetscStrlcat(p,((PetscObject)snes)->prefix,sizeof(p));CHKERRQ(ierr);} 250 if (flg) { 251 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," Matrix-free Options (via SNES):\n");CHKERRQ(ierr); 252 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); 253 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_umin <umin>: see users manual (default %g)\n",p,(double)mfctx->umin);CHKERRQ(ierr); 254 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_jorge: use Jorge More's method\n",p);CHKERRQ(ierr); 255 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);CHKERRQ(ierr); 256 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);CHKERRQ(ierr); 257 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_noise_file <file>: set file for printing noise info\n",p);CHKERRQ(ierr); 258 } 259 ierr = VecDuplicate(x,&mfctx->w);CHKERRQ(ierr); 260 ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr); 261 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 262 ierr = VecGetLocalSize(x,&nloc);CHKERRQ(ierr); 263 ierr = MatCreate(comm,J);CHKERRQ(ierr); 264 ierr = MatSetSizes(*J,nloc,n,n,n);CHKERRQ(ierr); 265 ierr = MatSetType(*J,MATSHELL);CHKERRQ(ierr); 266 ierr = MatShellSetContext(*J,mfctx);CHKERRQ(ierr); 267 ierr = MatShellSetOperation(*J,MATOP_MULT,(void (*)(void))SNESMatrixFreeMult2_Private);CHKERRQ(ierr); 268 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void (*)(void))SNESMatrixFreeDestroy2_Private);CHKERRQ(ierr); 269 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void (*)(void))SNESMatrixFreeView2_Private);CHKERRQ(ierr); 270 ierr = MatSetUp(*J);CHKERRQ(ierr); 271 272 ierr = PetscLogObjectParent((PetscObject)*J,(PetscObject)mfctx->w);CHKERRQ(ierr); 273 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)*J);CHKERRQ(ierr); 274 PetscFunctionReturn(0); 275 } 276 277 /*@C 278 SNESDefaultMatrixFreeSetParameters2 - Sets the parameters for the approximation of 279 matrix-vector products using finite differences. 280 281 $ J(u)*a = [J(u+h*a) - J(u)]/h where 282 283 either the user sets h directly here, or this parameter is computed via 284 285 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 286 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 287 $ 288 289 Input Parameters: 290 + mat - the matrix 291 . error_rel - relative error (should be set to the square root of 292 the relative error in the function evaluations) 293 . umin - minimum allowable u-value 294 - h - differencing parameter 295 296 Level: advanced 297 298 Notes: 299 If the user sets the parameter h directly, then this value will be used 300 instead of the default computation indicated above. 301 302 .keywords: SNES, matrix-free, parameters 303 304 .seealso: MatCreateSNESMF() 305 @*/ 306 PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat mat,PetscReal error,PetscReal umin,PetscReal h) 307 { 308 MFCtx_Private *ctx; 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 ierr = MatShellGetContext(mat,(void**)&ctx);CHKERRQ(ierr); 313 if (ctx) { 314 if (error != PETSC_DEFAULT) ctx->error_rel = error; 315 if (umin != PETSC_DEFAULT) ctx->umin = umin; 316 if (h != PETSC_DEFAULT) { 317 ctx->h = h; 318 ctx->need_h = PETSC_FALSE; 319 } 320 } 321 PetscFunctionReturn(0); 322 } 323 324 PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes) 325 { 326 MFCtx_Private *ctx; 327 PetscErrorCode ierr; 328 Mat mat; 329 330 PetscFunctionBegin; 331 ierr = SNESGetJacobian(snes,&mat,NULL,NULL,NULL);CHKERRQ(ierr); 332 ierr = MatShellGetContext(mat,(void**)&ctx);CHKERRQ(ierr); 333 if (ctx) ctx->need_h = PETSC_TRUE; 334 PetscFunctionReturn(0); 335 } 336 337