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 .seealso: MatDestroy(), MatMFFDSetFunctionError() 206 @*/ 207 PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J) 208 { 209 MPI_Comm comm; 210 MFCtx_Private *mfctx; 211 PetscErrorCode ierr; 212 PetscInt n,nloc; 213 PetscBool flg; 214 char p[64]; 215 216 PetscFunctionBegin; 217 ierr = PetscNewLog(snes,&mfctx);CHKERRQ(ierr); 218 mfctx->sp = 0; 219 mfctx->snes = snes; 220 mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON; 221 mfctx->umin = 1.e-6; 222 mfctx->h = 0.0; 223 mfctx->need_h = PETSC_TRUE; 224 mfctx->need_err = PETSC_FALSE; 225 mfctx->compute_err = PETSC_FALSE; 226 mfctx->compute_err_freq = 0; 227 mfctx->compute_err_iter = -1; 228 mfctx->compute_err = PETSC_FALSE; 229 mfctx->jorge = PETSC_FALSE; 230 231 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_err",&mfctx->error_rel,NULL);CHKERRQ(ierr); 232 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_umin",&mfctx->umin,NULL);CHKERRQ(ierr); 233 ierr = PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_jorge",&mfctx->jorge,NULL);CHKERRQ(ierr); 234 ierr = PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_compute_err",&mfctx->compute_err,NULL);CHKERRQ(ierr); 235 ierr = PetscOptionsGetInt(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_freq_err",&mfctx->compute_err_freq,&flg);CHKERRQ(ierr); 236 if (flg) { 237 if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0; 238 mfctx->compute_err = PETSC_TRUE; 239 } 240 if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE; 241 if (mfctx->jorge || mfctx->compute_err) { 242 ierr = SNESDiffParameterCreate_More(snes,x,&mfctx->data);CHKERRQ(ierr); 243 } else mfctx->data = 0; 244 245 ierr = PetscOptionsHasHelp(((PetscObject)snes)->options,&flg);CHKERRQ(ierr); 246 ierr = PetscStrncpy(p,"-",sizeof(p));CHKERRQ(ierr); 247 if (((PetscObject)snes)->prefix) {ierr = PetscStrlcat(p,((PetscObject)snes)->prefix,sizeof(p));CHKERRQ(ierr);} 248 if (flg) { 249 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," Matrix-free Options (via SNES):\n");CHKERRQ(ierr); 250 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); 251 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_umin <umin>: see users manual (default %g)\n",p,(double)mfctx->umin);CHKERRQ(ierr); 252 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_jorge: use Jorge More's method\n",p);CHKERRQ(ierr); 253 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);CHKERRQ(ierr); 254 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);CHKERRQ(ierr); 255 ierr = PetscPrintf(PetscObjectComm((PetscObject)snes)," %ssnes_mf_noise_file <file>: set file for printing noise info\n",p);CHKERRQ(ierr); 256 } 257 ierr = VecDuplicate(x,&mfctx->w);CHKERRQ(ierr); 258 ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr); 259 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 260 ierr = VecGetLocalSize(x,&nloc);CHKERRQ(ierr); 261 ierr = MatCreate(comm,J);CHKERRQ(ierr); 262 ierr = MatSetSizes(*J,nloc,n,n,n);CHKERRQ(ierr); 263 ierr = MatSetType(*J,MATSHELL);CHKERRQ(ierr); 264 ierr = MatShellSetContext(*J,mfctx);CHKERRQ(ierr); 265 ierr = MatShellSetOperation(*J,MATOP_MULT,(void (*)(void))SNESMatrixFreeMult2_Private);CHKERRQ(ierr); 266 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void (*)(void))SNESMatrixFreeDestroy2_Private);CHKERRQ(ierr); 267 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void (*)(void))SNESMatrixFreeView2_Private);CHKERRQ(ierr); 268 ierr = MatSetUp(*J);CHKERRQ(ierr); 269 270 ierr = PetscLogObjectParent((PetscObject)*J,(PetscObject)mfctx->w);CHKERRQ(ierr); 271 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)*J);CHKERRQ(ierr); 272 PetscFunctionReturn(0); 273 } 274 275 /*@C 276 SNESDefaultMatrixFreeSetParameters2 - Sets the parameters for the approximation of 277 matrix-vector products using finite differences. 278 279 $ J(u)*a = [J(u+h*a) - J(u)]/h where 280 281 either the user sets h directly here, or this parameter is computed via 282 283 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 284 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 285 $ 286 287 Input Parameters: 288 + mat - the matrix 289 . error_rel - relative error (should be set to the square root of 290 the relative error in the function evaluations) 291 . umin - minimum allowable u-value 292 - h - differencing parameter 293 294 Level: advanced 295 296 Notes: 297 If the user sets the parameter h directly, then this value will be used 298 instead of the default computation indicated above. 299 300 .seealso: MatCreateSNESMF() 301 @*/ 302 PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat mat,PetscReal error,PetscReal umin,PetscReal h) 303 { 304 MFCtx_Private *ctx; 305 PetscErrorCode ierr; 306 307 PetscFunctionBegin; 308 ierr = MatShellGetContext(mat,(void**)&ctx);CHKERRQ(ierr); 309 if (ctx) { 310 if (error != PETSC_DEFAULT) ctx->error_rel = error; 311 if (umin != PETSC_DEFAULT) ctx->umin = umin; 312 if (h != PETSC_DEFAULT) { 313 ctx->h = h; 314 ctx->need_h = PETSC_FALSE; 315 } 316 } 317 PetscFunctionReturn(0); 318 } 319 320 PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes) 321 { 322 MFCtx_Private *ctx; 323 PetscErrorCode ierr; 324 Mat mat; 325 326 PetscFunctionBegin; 327 ierr = SNESGetJacobian(snes,&mat,NULL,NULL,NULL);CHKERRQ(ierr); 328 ierr = MatShellGetContext(mat,(void**)&ctx);CHKERRQ(ierr); 329 if (ctx) ctx->need_h = PETSC_TRUE; 330 PetscFunctionReturn(0); 331 } 332 333