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