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