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 { 63 MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx; 64 PetscReal nrm, sum, umin = hctx->umin; 65 PetscScalar dot; 66 67 PetscFunctionBegin; 68 if (!(ctx->count % ctx->recomputeperiod)) { 69 /* 70 This algorithm requires 2 norms and 1 inner product. Rather than 71 use directly the VecNorm() and VecDot() routines (and thus have 72 three separate collective operations, we use the VecxxxBegin/End() routines 73 */ 74 PetscCall(VecDotBegin(U, a, &dot)); 75 PetscCall(VecNormBegin(a, NORM_1, &sum)); 76 PetscCall(VecNormBegin(a, NORM_2, &nrm)); 77 PetscCall(VecDotEnd(U, a, &dot)); 78 PetscCall(VecNormEnd(a, NORM_1, &sum)); 79 PetscCall(VecNormEnd(a, NORM_2, &nrm)); 80 81 if (nrm == 0.0) { 82 *zeroa = PETSC_TRUE; 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 *zeroa = PETSC_FALSE; 86 87 /* 88 Safeguard for step sizes that are "too small" 89 */ 90 if (PetscAbsScalar(dot) < umin * sum && PetscRealPart(dot) >= 0.0) dot = umin * sum; 91 else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin * sum) dot = -umin * sum; 92 *h = ctx->error_rel * dot / (nrm * nrm); 93 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); 94 } else { 95 *h = ctx->currenth; 96 } 97 ctx->count++; 98 PetscFunctionReturn(PETSC_SUCCESS); 99 } 100 101 /* 102 MatMFFDView_DS - Prints information about this particular 103 method for computing h. Note that this does not print the general 104 information about the matrix-free method, as such info is printed 105 by the calling routine. 106 107 Input Parameters: 108 + ctx - the matrix free context 109 - viewer - the PETSc viewer 110 */ 111 static PetscErrorCode MatMFFDView_DS(MatMFFD ctx, PetscViewer viewer) 112 { 113 MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx; 114 PetscBool iascii; 115 116 PetscFunctionBegin; 117 /* 118 Currently this only handles the ascii file viewers, others 119 could be added, but for this type of object other viewers 120 make less sense 121 */ 122 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 123 if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " umin=%g (minimum iterate parameter)\n", (double)hctx->umin)); 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 /* 128 MatMFFDSetFromOptions_DS - Looks in the options database for 129 any options appropriate for this method. 130 131 Input Parameter: 132 . ctx - the matrix free context 133 134 */ 135 static PetscErrorCode MatMFFDSetFromOptions_DS(MatMFFD ctx, PetscOptionItems *PetscOptionsObject) 136 { 137 MatMFFD_DS *hctx = (MatMFFD_DS *)ctx->hctx; 138 139 PetscFunctionBegin; 140 PetscOptionsHeadBegin(PetscOptionsObject, "Finite difference matrix free parameters"); 141 PetscCall(PetscOptionsReal("-mat_mffd_umin", "umin", "MatMFFDDSSetUmin", hctx->umin, &hctx->umin, NULL)); 142 PetscOptionsHeadEnd(); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 /* 147 MatMFFDDestroy_DS - Frees the space allocated by 148 MatCreateMFFD_DS(). 149 150 Input Parameter: 151 . ctx - the matrix free context 152 153 Note: 154 Does not free the ctx, that is handled by the calling routine 155 */ 156 static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx) 157 { 158 PetscFunctionBegin; 159 PetscCall(PetscFree(ctx->hctx)); 160 PetscFunctionReturn(PETSC_SUCCESS); 161 } 162 163 /* 164 The following two routines use the PetscObjectCompose() and PetscObjectQuery() 165 mechanism to allow the user to change the Umin parameter used in this method. 166 */ 167 PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat, PetscReal umin) 168 { 169 MatMFFD ctx = NULL; 170 MatMFFD_DS *hctx; 171 172 PetscFunctionBegin; 173 PetscCall(MatShellGetContext(mat, &ctx)); 174 PetscCheck(ctx, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatMFFDDSSetUmin() attached to non-shell matrix"); 175 hctx = (MatMFFD_DS *)ctx->hctx; 176 hctx->umin = umin; 177 PetscFunctionReturn(PETSC_SUCCESS); 178 } 179 180 /*@ 181 MatMFFDDSSetUmin - Sets the "umin" parameter used by the 182 PETSc routine for computing the differencing parameter, h, which is used 183 for matrix-free Jacobian-vector products for a `MATMFFD` matrix. 184 185 Input Parameters: 186 + A - the `MATMFFD` matrix 187 - umin - the parameter 188 189 Level: advanced 190 191 Note: 192 See the manual page for `MatCreateSNESMF()` for a complete description of the 193 algorithm used to compute h. 194 195 .seealso: `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()` 196 @*/ 197 PetscErrorCode MatMFFDDSSetUmin(Mat A, PetscReal umin) 198 { 199 PetscFunctionBegin; 200 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 201 PetscTryMethod(A, "MatMFFDDSSetUmin_C", (Mat, PetscReal), (A, umin)); 202 PetscFunctionReturn(PETSC_SUCCESS); 203 } 204 205 /*MC 206 MATMFFD_DS - algorithm for compute the "h" used in the finite difference matrix-free matrix vector product, `MatMult()`. 207 208 Options Database Keys: 209 . -mat_mffd_umin <umin> - see `MatMFFDDSSetUmin()` 210 211 Level: intermediate 212 213 Notes: 214 Requires 2 norms and 1 inner product, but they are computed together 215 so only one parallel collective operation is needed. See `MATMFFD_WP` for a method 216 (with `KSPGMRES`) that requires NO collective operations. 217 218 Formula used: 219 F'(u)*a = [F(u+h*a) - F(u)]/h where 220 h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 221 = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 222 where 223 error_rel = square root of relative error in function evaluation 224 umin = minimum iterate parameter 225 226 References: 227 . * - Dennis and Schnabel, "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" 228 229 .seealso: `MATMFFD`, `MATMFFD_WP`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_WP`, `MatMFFDDSSetUmin()` 230 M*/ 231 PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx) 232 { 233 MatMFFD_DS *hctx; 234 235 PetscFunctionBegin; 236 /* allocate my own private data structure */ 237 PetscCall(PetscNew(&hctx)); 238 ctx->hctx = (void *)hctx; 239 /* set a default for my parameter */ 240 hctx->umin = 1.e-6; 241 242 /* set the functions I am providing */ 243 ctx->ops->compute = MatMFFDCompute_DS; 244 ctx->ops->destroy = MatMFFDDestroy_DS; 245 ctx->ops->view = MatMFFDView_DS; 246 ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS; 247 248 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDDSSetUmin_C", MatMFFDDSSetUmin_DS)); 249 PetscFunctionReturn(PETSC_SUCCESS); 250 } 251