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 ctx->error_rel = sqrt(noise); 124 ierr = PetscInfo3(snes,"Using Jorge's noise: noise=%G, sqrt(noise)=%G, h_more=%G\n",noise,ctx->error_rel,h);CHKERRQ(ierr); 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) {dot = 1.0; norm = 1.0;} 139 else if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum; 140 else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum; 141 h = PetscRealPart(ctx->error_rel*dot/(norm*norm)); 142 } 143 } else { 144 h = ctx->h; 145 } 146 if (!ctx->jorge || !ctx->need_h) {ierr = PetscInfo1(snes,"h = %G\n",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 (ctx->sp) {ierr = MatNullSpaceRemove(ctx->sp,y,PETSC_NULL);CHKERRQ(ierr);} 155 156 ierr = PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);CHKERRQ(ierr); 157 PetscFunctionReturn(0); 158 } 159 160 #undef __FUNCT__ 161 #define __FUNCT__ "SNESDefaultMatrixFreeCreate2" 162 /*@C 163 SNESMatrixFreeCreate2 - Creates a matrix-free matrix 164 context for use with a SNES solver. This matrix can be used as 165 the Jacobian argument for the routine SNESSetJacobian(). 166 167 Input Parameters: 168 . snes - the SNES context 169 . x - vector where SNES solution is to be stored. 170 171 Output Parameter: 172 . J - the matrix-free matrix 173 174 Level: advanced 175 176 Notes: 177 The matrix-free matrix context merely contains the function pointers 178 and work space for performing finite difference approximations of 179 Jacobian-vector products, J(u)*a, via 180 181 $ J(u)*a = [J(u+h*a) - J(u)]/h, 182 $ where by default 183 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 184 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 185 $ where 186 $ error_rel = square root of relative error in 187 $ function evaluation 188 $ umin = minimum iterate parameter 189 $ Alternatively, the differencing parameter, h, can be set using 190 $ Jorge's nifty new strategy if one specifies the option 191 $ -snes_mf_jorge 192 193 The user can set these parameters via MatMFFDSetFunctionError(). 194 See the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>. 195 196 The user should call MatDestroy() when finished with the matrix-free 197 matrix context. 198 199 Options Database Keys: 200 $ -snes_mf_err <error_rel> 201 $ -snes_mf_unim <umin> 202 $ -snes_mf_compute_err 203 $ -snes_mf_freq_err <freq> 204 $ -snes_mf_jorge 205 206 .keywords: SNES, default, matrix-free, create, matrix 207 208 .seealso: MatDestroy(), MatMFFDSetFunctionError() 209 @*/ 210 PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J) 211 { 212 MPI_Comm comm; 213 MFCtx_Private *mfctx; 214 PetscErrorCode ierr; 215 PetscInt n,nloc; 216 PetscBool flg; 217 char p[64]; 218 219 PetscFunctionBegin; 220 ierr = PetscNewLog(snes,MFCtx_Private,&mfctx);CHKERRQ(ierr); 221 mfctx->sp = 0; 222 mfctx->snes = snes; 223 mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON; 224 mfctx->umin = 1.e-6; 225 mfctx->h = 0.0; 226 mfctx->need_h = PETSC_TRUE; 227 mfctx->need_err = PETSC_FALSE; 228 mfctx->compute_err = PETSC_FALSE; 229 mfctx->compute_err_freq = 0; 230 mfctx->compute_err_iter = -1; 231 mfctx->compute_err = PETSC_FALSE; 232 mfctx->jorge = PETSC_FALSE; 233 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_mf_err",&mfctx->error_rel,PETSC_NULL);CHKERRQ(ierr); 234 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_mf_umin",&mfctx->umin,PETSC_NULL);CHKERRQ(ierr); 235 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_mf_jorge",&mfctx->jorge,PETSC_NULL);CHKERRQ(ierr); 236 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_mf_compute_err",&mfctx->compute_err,PETSC_NULL);CHKERRQ(ierr); 237 ierr = PetscOptionsGetInt(((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 = PetscOptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 248 ierr = PetscStrcpy(p,"-");CHKERRQ(ierr); 249 if (((PetscObject)snes)->prefix) PetscStrcat(p,((PetscObject)snes)->prefix); 250 if (flg) { 251 ierr = PetscPrintf(((PetscObject)snes)->comm," Matrix-free Options (via SNES):\n");CHKERRQ(ierr); 252 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); 253 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_umin <umin>: see users manual (default %G)\n",p,mfctx->umin);CHKERRQ(ierr); 254 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_jorge: use Jorge More's method\n",p);CHKERRQ(ierr); 255 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);CHKERRQ(ierr); 256 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);CHKERRQ(ierr); 257 ierr = PetscPrintf(((PetscObject)snes)->comm," %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 271 ierr = PetscLogObjectParent(*J,mfctx->w);CHKERRQ(ierr); 272 ierr = PetscLogObjectParent(snes,*J);CHKERRQ(ierr); 273 PetscFunctionReturn(0); 274 } 275 276 #undef __FUNCT__ 277 #define __FUNCT__ "SNESDefaultMatrixFreeSetParameters2" 278 /*@C 279 SNESDefaultMatrixFreeSetParameters2 - Sets the parameters for the approximation of 280 matrix-vector products using finite differences. 281 282 $ J(u)*a = [J(u+h*a) - J(u)]/h where 283 284 either the user sets h directly here, or this parameter is computed via 285 286 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 287 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 288 $ 289 290 Input Parameters: 291 + mat - the matrix 292 . error_rel - relative error (should be set to the square root of 293 the relative error in the function evaluations) 294 . umin - minimum allowable u-value 295 - h - differencing parameter 296 297 Level: advanced 298 299 Notes: 300 If the user sets the parameter h directly, then this value will be used 301 instead of the default computation indicated above. 302 303 .keywords: SNES, matrix-free, parameters 304 305 .seealso: MatCreateSNESMF() 306 @*/ 307 PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat mat,PetscReal error,PetscReal umin,PetscReal h) 308 { 309 MFCtx_Private *ctx; 310 PetscErrorCode ierr; 311 312 PetscFunctionBegin; 313 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 314 if (ctx) { 315 if (error != PETSC_DEFAULT) ctx->error_rel = error; 316 if (umin != PETSC_DEFAULT) ctx->umin = umin; 317 if (h != PETSC_DEFAULT) { 318 ctx->h = h; 319 ctx->need_h = PETSC_FALSE; 320 } 321 } 322 PetscFunctionReturn(0); 323 } 324 325 PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes) 326 { 327 MFCtx_Private *ctx; 328 PetscErrorCode ierr; 329 Mat mat; 330 331 PetscFunctionBegin; 332 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 333 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 334 if (ctx) ctx->need_h = PETSC_TRUE; 335 PetscFunctionReturn(0); 336 } 337 338