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