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