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 Notes: 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. 179 180 Input Parameters: 181 + A - the matrix created with MatCreateSNESMF() 182 - umin - the parameter 183 184 Level: advanced 185 186 Notes: 187 See the manual page for MatCreateSNESMF() for a complete description of the 188 algorithm used to compute h. 189 190 .seealso: `MatMFFDSetFunctionError()`, `MatCreateSNESMF()` 191 192 @*/ 193 PetscErrorCode MatMFFDDSSetUmin(Mat A, PetscReal umin) { 194 PetscFunctionBegin; 195 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 196 PetscTryMethod(A, "MatMFFDDSSetUmin_C", (Mat, PetscReal), (A, umin)); 197 PetscFunctionReturn(0); 198 } 199 200 /*MC 201 MATMFFD_DS - the code for compute the "h" used in the finite difference 202 matrix-free matrix vector product. This code 203 implements the strategy in Dennis and Schnabel, "Numerical Methods for Unconstrained 204 Optimization and Nonlinear Equations". 205 206 Options Database Keys: 207 . -mat_mffd_umin <umin> - see MatMFFDDSSetUmin() 208 209 Level: intermediate 210 211 Notes: 212 Requires 2 norms and 1 inner product, but they are computed together 213 so only one parallel collective operation is needed. See MATMFFD_WP for a method 214 (with GMRES) that requires NO collective operations. 215 216 Formula used: 217 F'(u)*a = [F(u+h*a) - F(u)]/h where 218 h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 219 = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 220 where 221 error_rel = square root of relative error in function evaluation 222 umin = minimum iterate parameter 223 224 .seealso: `MATMFFD`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_WP`, `MatMFFDDSSetUmin()` 225 226 M*/ 227 PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx) { 228 MatMFFD_DS *hctx; 229 230 PetscFunctionBegin; 231 /* allocate my own private data structure */ 232 PetscCall(PetscNewLog(ctx, &hctx)); 233 ctx->hctx = (void *)hctx; 234 /* set a default for my parameter */ 235 hctx->umin = 1.e-6; 236 237 /* set the functions I am providing */ 238 ctx->ops->compute = MatMFFDCompute_DS; 239 ctx->ops->destroy = MatMFFDDestroy_DS; 240 ctx->ops->view = MatMFFDView_DS; 241 ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS; 242 243 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDDSSetUmin_C", MatMFFDDSSetUmin_DS)); 244 PetscFunctionReturn(0); 245 } 246