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); 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 } else { 68 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by SNES matrix free Jorge",((PetscObject)viewer)->type_name); 69 } 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "SNESMatrixFreeMult2_Private" 75 /* 76 SNESMatrixFreeMult2_Private - Default matrix-free form for Jacobian-vector 77 product, y = F'(u)*a: 78 y = (F(u + ha) - F(u)) /h, 79 where F = nonlinear function, as set by SNESSetFunction() 80 u = current iterate 81 h = difference interval 82 */ 83 PetscErrorCode SNESMatrixFreeMult2_Private(Mat mat,Vec a,Vec y) 84 { 85 MFCtx_Private *ctx; 86 SNES snes; 87 PetscReal h,norm,sum,umin,noise; 88 PetscScalar hs,dot; 89 Vec w,U,F; 90 PetscErrorCode ierr,(*eval_fct)(SNES,Vec,Vec); 91 MPI_Comm comm; 92 PetscInt iter; 93 94 PetscFunctionBegin; 95 96 /* We log matrix-free matrix-vector products separately, so that we can 97 separate the performance monitoring from the cases that use conventional 98 storage. We may eventually modify event logging to associate events 99 with particular objects, hence alleviating the more general problem. */ 100 ierr = PetscLogEventBegin(MATMFFD_Mult,a,y,0,0);CHKERRQ(ierr); 101 102 ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr); 103 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 104 snes = ctx->snes; 105 w = ctx->w; 106 umin = ctx->umin; 107 108 ierr = SNESGetSolution(snes,&U);CHKERRQ(ierr); 109 eval_fct = SNESComputeFunction; 110 ierr = SNESGetFunction(snes,&F,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 111 112 /* Determine a "good" step size, h */ 113 if (ctx->need_h) { 114 115 /* Use Jorge's method to compute h */ 116 if (ctx->jorge) { 117 ierr = SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);CHKERRQ(ierr); 118 119 /* Use the Brown/Saad method to compute h */ 120 } else { 121 /* Compute error if desired */ 122 ierr = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr); 123 if ((ctx->need_err) || ((ctx->compute_err_freq) && (ctx->compute_err_iter != iter) && (!((iter-1)%ctx->compute_err_freq)))) { 124 /* Use Jorge's method to compute noise */ 125 ierr = SNESDiffParameterCompute_More(snes,ctx->data,U,a,&noise,&h);CHKERRQ(ierr); 126 ctx->error_rel = sqrt(noise); 127 ierr = PetscInfo3(snes,"Using Jorge's noise: noise=%G, sqrt(noise)=%G, h_more=%G\n",noise,ctx->error_rel,h);CHKERRQ(ierr); 128 ctx->compute_err_iter = iter; 129 ctx->need_err = PETSC_FALSE; 130 } 131 132 ierr = VecDotBegin(U,a,&dot);CHKERRQ(ierr); 133 ierr = VecNormBegin(a,NORM_1,&sum);CHKERRQ(ierr); 134 ierr = VecNormBegin(a,NORM_2,&norm);CHKERRQ(ierr); 135 ierr = VecDotEnd(U,a,&dot);CHKERRQ(ierr); 136 ierr = VecNormEnd(a,NORM_1,&sum);CHKERRQ(ierr); 137 ierr = VecNormEnd(a,NORM_2,&norm);CHKERRQ(ierr); 138 139 140 /* Safeguard for step sizes too small */ 141 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 142 else if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum; 143 else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum; 144 h = PetscRealPart(ctx->error_rel*dot/(norm*norm)); 145 } 146 } else { 147 h = ctx->h; 148 } 149 if (!ctx->jorge || !ctx->need_h) {ierr = PetscInfo1(snes,"h = %G\n",h);CHKERRQ(ierr);} 150 151 /* Evaluate function at F(u + ha) */ 152 hs = h; 153 ierr = VecWAXPY(w,hs,a,U);CHKERRQ(ierr); 154 ierr = eval_fct(snes,w,y);CHKERRQ(ierr); 155 ierr = VecAXPY(y,-1.0,F);CHKERRQ(ierr); 156 ierr = VecScale(y,1.0/hs);CHKERRQ(ierr); 157 if (ctx->sp) {ierr = MatNullSpaceRemove(ctx->sp,y,PETSC_NULL);CHKERRQ(ierr);} 158 159 ierr = PetscLogEventEnd(MATMFFD_Mult,a,y,0,0);CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "SNESDefaultMatrixFreeCreate2" 165 /*@C 166 SNESMatrixFreeCreate2 - Creates a matrix-free matrix 167 context for use with a SNES solver. This matrix can be used as 168 the Jacobian argument for the routine SNESSetJacobian(). 169 170 Input Parameters: 171 . snes - the SNES context 172 . x - vector where SNES solution is to be stored. 173 174 Output Parameter: 175 . J - the matrix-free matrix 176 177 Level: advanced 178 179 Notes: 180 The matrix-free matrix context merely contains the function pointers 181 and work space for performing finite difference approximations of 182 Jacobian-vector products, J(u)*a, via 183 184 $ J(u)*a = [J(u+h*a) - J(u)]/h, 185 $ where by default 186 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 187 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 188 $ where 189 $ error_rel = square root of relative error in 190 $ function evaluation 191 $ umin = minimum iterate parameter 192 $ Alternatively, the differencing parameter, h, can be set using 193 $ Jorge's nifty new strategy if one specifies the option 194 $ -snes_mf_jorge 195 196 The user can set these parameters via MatMFFDSetFunctionError(). 197 See the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>. 198 199 The user should call MatDestroy() when finished with the matrix-free 200 matrix context. 201 202 Options Database Keys: 203 $ -snes_mf_err <error_rel> 204 $ -snes_mf_unim <umin> 205 $ -snes_mf_compute_err 206 $ -snes_mf_freq_err <freq> 207 $ -snes_mf_jorge 208 209 .keywords: SNES, default, matrix-free, create, matrix 210 211 .seealso: MatDestroy(), MatMFFDSetFunctionError() 212 @*/ 213 PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES snes,Vec x,Mat *J) 214 { 215 MPI_Comm comm; 216 MFCtx_Private *mfctx; 217 PetscErrorCode ierr; 218 PetscInt n,nloc; 219 PetscBool flg; 220 char p[64]; 221 222 PetscFunctionBegin; 223 ierr = PetscNewLog(snes,MFCtx_Private,&mfctx);CHKERRQ(ierr); 224 mfctx->sp = 0; 225 mfctx->snes = snes; 226 mfctx->error_rel = PETSC_SQRT_MACHINE_EPSILON; 227 mfctx->umin = 1.e-6; 228 mfctx->h = 0.0; 229 mfctx->need_h = PETSC_TRUE; 230 mfctx->need_err = PETSC_FALSE; 231 mfctx->compute_err = PETSC_FALSE; 232 mfctx->compute_err_freq = 0; 233 mfctx->compute_err_iter = -1; 234 mfctx->compute_err = PETSC_FALSE; 235 mfctx->jorge = PETSC_FALSE; 236 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_mf_err",&mfctx->error_rel,PETSC_NULL);CHKERRQ(ierr); 237 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_mf_umin",&mfctx->umin,PETSC_NULL);CHKERRQ(ierr); 238 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_mf_jorge",&mfctx->jorge,PETSC_NULL);CHKERRQ(ierr); 239 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_mf_compute_err",&mfctx->compute_err,PETSC_NULL);CHKERRQ(ierr); 240 ierr = PetscOptionsGetInt(((PetscObject)snes)->prefix,"-snes_mf_freq_err",&mfctx->compute_err_freq,&flg);CHKERRQ(ierr); 241 if (flg) { 242 if (mfctx->compute_err_freq < 0) mfctx->compute_err_freq = 0; 243 mfctx->compute_err = PETSC_TRUE; 244 } 245 if (mfctx->compute_err) mfctx->need_err = PETSC_TRUE; 246 if (mfctx->jorge || mfctx->compute_err) { 247 ierr = SNESDiffParameterCreate_More(snes,x,&mfctx->data);CHKERRQ(ierr); 248 } else mfctx->data = 0; 249 250 ierr = PetscOptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 251 ierr = PetscStrcpy(p,"-");CHKERRQ(ierr); 252 if (((PetscObject)snes)->prefix) PetscStrcat(p,((PetscObject)snes)->prefix); 253 if (flg) { 254 ierr = PetscPrintf(((PetscObject)snes)->comm," Matrix-free Options (via SNES):\n");CHKERRQ(ierr); 255 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); 256 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_umin <umin>: see users manual (default %G)\n",p,mfctx->umin);CHKERRQ(ierr); 257 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_jorge: use Jorge More's method\n",p);CHKERRQ(ierr); 258 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_compute_err: compute sqrt or relative error in function\n",p);CHKERRQ(ierr); 259 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_freq_err <freq>: frequency to recompute this (default only once)\n",p);CHKERRQ(ierr); 260 ierr = PetscPrintf(((PetscObject)snes)->comm," %ssnes_mf_noise_file <file>: set file for printing noise info\n",p);CHKERRQ(ierr); 261 } 262 ierr = VecDuplicate(x,&mfctx->w);CHKERRQ(ierr); 263 ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr); 264 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 265 ierr = VecGetLocalSize(x,&nloc);CHKERRQ(ierr); 266 ierr = MatCreate(comm,J);CHKERRQ(ierr); 267 ierr = MatSetSizes(*J,nloc,n,n,n);CHKERRQ(ierr); 268 ierr = MatSetType(*J,MATSHELL);CHKERRQ(ierr); 269 ierr = MatShellSetContext(*J,mfctx);CHKERRQ(ierr); 270 ierr = MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))SNESMatrixFreeMult2_Private);CHKERRQ(ierr); 271 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))SNESMatrixFreeDestroy2_Private);CHKERRQ(ierr); 272 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))SNESMatrixFreeView2_Private);CHKERRQ(ierr); 273 274 ierr = PetscLogObjectParent(*J,mfctx->w);CHKERRQ(ierr); 275 ierr = PetscLogObjectParent(snes,*J);CHKERRQ(ierr); 276 PetscFunctionReturn(0); 277 } 278 279 #undef __FUNCT__ 280 #define __FUNCT__ "SNESDefaultMatrixFreeSetParameters2" 281 /*@C 282 SNESDefaultMatrixFreeSetParameters2 - Sets the parameters for the approximation of 283 matrix-vector products using finite differences. 284 285 $ J(u)*a = [J(u+h*a) - J(u)]/h where 286 287 either the user sets h directly here, or this parameter is computed via 288 289 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 290 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 291 $ 292 293 Input Parameters: 294 + mat - the matrix 295 . error_rel - relative error (should be set to the square root of 296 the relative error in the function evaluations) 297 . umin - minimum allowable u-value 298 - h - differencing parameter 299 300 Level: advanced 301 302 Notes: 303 If the user sets the parameter h directly, then this value will be used 304 instead of the default computation indicated above. 305 306 .keywords: SNES, matrix-free, parameters 307 308 .seealso: MatCreateSNESMF() 309 @*/ 310 PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat mat,PetscReal error,PetscReal umin,PetscReal h) 311 { 312 MFCtx_Private *ctx; 313 PetscErrorCode ierr; 314 315 PetscFunctionBegin; 316 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 317 if (ctx) { 318 if (error != PETSC_DEFAULT) ctx->error_rel = error; 319 if (umin != PETSC_DEFAULT) ctx->umin = umin; 320 if (h != PETSC_DEFAULT) { 321 ctx->h = h; 322 ctx->need_h = PETSC_FALSE; 323 } 324 } 325 PetscFunctionReturn(0); 326 } 327 328 PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes) 329 { 330 MFCtx_Private *ctx; 331 PetscErrorCode ierr; 332 Mat mat; 333 334 PetscFunctionBegin; 335 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 336 ierr = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(ierr); 337 if (ctx) ctx->need_h = PETSC_TRUE; 338 PetscFunctionReturn(0); 339 } 340 341