1 2 /* 3 Implements the DS PETSc approach for computing the h 4 parameter used with the finite difference based matrix-free 5 Jacobian-vector products. 6 7 To make your own: clone this file and modify for your needs. 8 9 Mandatory functions: 10 ------------------- 11 MatMFFDCompute_ - for a given point and direction computes h 12 13 MatCreateMFFD _ - fills in the MatMFFD data structure 14 for this particular implementation 15 16 Optional functions: 17 ------------------- 18 MatMFFDView_ - prints information about the parameters being used. 19 This is called when SNESView() or -snes_view is used. 20 21 MatMFFDSetFromOptions_ - checks the options database for options that 22 apply to this method. 23 24 MatMFFDDestroy_ - frees any space allocated by the routines above 25 26 */ 27 28 /* 29 This include file defines the data structure MatMFFD that 30 includes information about the computation of h. It is shared by 31 all implementations that people provide 32 */ 33 #include <petsc/private/matimpl.h> 34 #include <../src/mat/impls/mffd/mffdimpl.h> /*I "petscmat.h" I*/ 35 36 /* 37 The method has one parameter that is used to 38 "cutoff" very small values. This is stored in a data structure 39 that is only visible to this file. If your method has no parameters 40 it can omit this, if it has several simply reorganize the data structure. 41 The data structure is "hung-off" the MatMFFD data structure in 42 the void *hctx; field. 43 */ 44 typedef struct { 45 PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */ 46 } MatMFFD_DS; 47 48 /* 49 MatMFFDCompute_DS - Standard PETSc code for computing the 50 differencing parameter (h) for use with matrix-free finite differences. 51 52 Input Parameters: 53 + ctx - the matrix free context 54 . U - the location at which you want the Jacobian 55 - a - the direction you want the derivative 56 57 Output Parameter: 58 . h - the scale computed 59 60 */ 61 static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx, Vec U, Vec a, PetscScalar *h, PetscBool *zeroa) { 62 MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx; 63 PetscReal nrm, sum, umin = hctx->umin; 64 PetscScalar dot; 65 66 PetscFunctionBegin; 67 if (!(ctx->count % ctx->recomputeperiod)) { 68 /* 69 This algorithm requires 2 norms and 1 inner product. Rather than 70 use directly the VecNorm() and VecDot() routines (and thus have 71 three separate collective operations, we use the VecxxxBegin/End() routines 72 */ 73 PetscCall(VecDotBegin(U, a, &dot)); 74 PetscCall(VecNormBegin(a, NORM_1, &sum)); 75 PetscCall(VecNormBegin(a, NORM_2, &nrm)); 76 PetscCall(VecDotEnd(U, a, &dot)); 77 PetscCall(VecNormEnd(a, NORM_1, &sum)); 78 PetscCall(VecNormEnd(a, NORM_2, &nrm)); 79 80 if (nrm == 0.0) { 81 *zeroa = PETSC_TRUE; 82 PetscFunctionReturn(0); 83 } 84 *zeroa = PETSC_FALSE; 85 86 /* 87 Safeguard for step sizes that are "too small" 88 */ 89 if (PetscAbsScalar(dot) < umin * sum && PetscRealPart(dot) >= 0.0) dot = umin * sum; 90 else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin * sum) dot = -umin * sum; 91 *h = ctx->error_rel * dot / (nrm * nrm); 92 PetscCheck(!PetscIsInfOrNanScalar(*h), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Differencing parameter is not a number sum = %g dot = %g norm = %g", (double)sum, (double)PetscRealPart(dot), (double)nrm); 93 } else { 94 *h = ctx->currenth; 95 } 96 ctx->count++; 97 PetscFunctionReturn(0); 98 } 99 100 /* 101 MatMFFDView_DS - Prints information about this particular 102 method for computing h. Note that this does not print the general 103 information about the matrix-free method, as such info is printed 104 by the calling routine. 105 106 Input Parameters: 107 + ctx - the matrix free context 108 - viewer - the PETSc viewer 109 */ 110 static PetscErrorCode MatMFFDView_DS(MatMFFD ctx, PetscViewer viewer) { 111 MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx; 112 PetscBool iascii; 113 114 PetscFunctionBegin; 115 /* 116 Currently this only handles the ascii file viewers, others 117 could be added, but for this type of object other viewers 118 make less sense 119 */ 120 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 121 if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " umin=%g (minimum iterate parameter)\n", (double)hctx->umin)); 122 PetscFunctionReturn(0); 123 } 124 125 /* 126 MatMFFDSetFromOptions_DS - Looks in the options database for 127 any options appropriate for this method. 128 129 Input Parameter: 130 . ctx - the matrix free context 131 132 */ 133 static PetscErrorCode MatMFFDSetFromOptions_DS(MatMFFD ctx, PetscOptionItems *PetscOptionsObject) { 134 MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx; 135 136 PetscFunctionBegin; 137 PetscOptionsHeadBegin(PetscOptionsObject, "Finite difference matrix free parameters"); 138 PetscCall(PetscOptionsReal("-mat_mffd_umin", "umin", "MatMFFDDSSetUmin", hctx->umin, &hctx->umin, NULL)); 139 PetscOptionsHeadEnd(); 140 PetscFunctionReturn(0); 141 } 142 143 /* 144 MatMFFDDestroy_DS - Frees the space allocated by 145 MatCreateMFFD_DS(). 146 147 Input Parameter: 148 . ctx - the matrix free context 149 150 Note: 151 Does not free the ctx, that is handled by the calling routine 152 */ 153 static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx) { 154 PetscFunctionBegin; 155 PetscCall(PetscFree(ctx->hctx)); 156 PetscFunctionReturn(0); 157 } 158 159 /* 160 The following two routines use the PetscObjectCompose() and PetscObjectQuery() 161 mechanism to allow the user to change the Umin parameter used in this method. 162 */ 163 PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat, PetscReal umin) { 164 MatMFFD ctx = NULL; 165 MatMFFD_DS *hctx; 166 167 PetscFunctionBegin; 168 PetscCall(MatShellGetContext(mat, &ctx)); 169 PetscCheck(ctx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatMFFDDSSetUmin() attached to non-shell matrix"); 170 hctx = (MatMFFD_DS *)ctx->hctx; 171 hctx->umin = umin; 172 PetscFunctionReturn(0); 173 } 174 175 /*@ 176 MatMFFDDSSetUmin - Sets the "umin" parameter used by the 177 PETSc routine for computing the differencing parameter, h, which is used 178 for matrix-free Jacobian-vector products for a `MATMFFD` matrix. 179 180 Input Parameters: 181 + A - the `MATMFFD` matrix 182 - umin - the parameter 183 184 Level: advanced 185 186 Note: 187 See the manual page for `MatCreateSNESMF()` for a complete description of the 188 algorithm used to compute h. 189 190 .seealso: `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()` 191 @*/ 192 PetscErrorCode MatMFFDDSSetUmin(Mat A, PetscReal umin) { 193 PetscFunctionBegin; 194 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 195 PetscTryMethod(A, "MatMFFDDSSetUmin_C", (Mat, PetscReal), (A, umin)); 196 PetscFunctionReturn(0); 197 } 198 199 /*MC 200 MATMFFD_DS - algorithm for compute the "h" used in the finite difference matrix-free matrix vector product, `MatMult()`. 201 202 Options Database Keys: 203 . -mat_mffd_umin <umin> - see `MatMFFDDSSetUmin()` 204 205 Level: intermediate 206 207 Notes: 208 Requires 2 norms and 1 inner product, but they are computed together 209 so only one parallel collective operation is needed. See `MATMFFD_WP` for a method 210 (with `KSPGMRES`) that requires NO collective operations. 211 212 Formula used: 213 F'(u)*a = [F(u+h*a) - F(u)]/h where 214 h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 215 = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 216 where 217 error_rel = square root of relative error in function evaluation 218 umin = minimum iterate parameter 219 220 References: 221 . * - Dennis and Schnabel, "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" 222 223 .seealso: `MATMFFD`, `MATMFFD_WP`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_WP`, `MatMFFDDSSetUmin()` 224 M*/ 225 PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx) { 226 MatMFFD_DS *hctx; 227 228 PetscFunctionBegin; 229 /* allocate my own private data structure */ 230 PetscCall(PetscNew(&hctx)); 231 ctx->hctx = (void *)hctx; 232 /* set a default for my parameter */ 233 hctx->umin = 1.e-6; 234 235 /* set the functions I am providing */ 236 ctx->ops->compute = MatMFFDCompute_DS; 237 ctx->ops->destroy = MatMFFDDestroy_DS; 238 ctx->ops->view = MatMFFDView_DS; 239 ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS; 240 241 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDDSSetUmin_C", MatMFFDDSSetUmin_DS)); 242 PetscFunctionReturn(0); 243 } 244