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